Question
Let $$ \pi_{rm_c}(x) = \sum_{ \substack{ {n\leq x}\\{(n+a,P(\sqrt{n}))=1}}} 1-1, $$ where $P(x)$ is the product of all primes less or equal to $x$ and $a$ is a random integer constrained to those values such that $(n+a,P(\sqrt{n}))\leq n$ for all $n\leq x$.
What is the asymptotic behaviour of the expected value of $\pi_{rm_c}(x)$?
Below follows some background for the question, a possible partial approach, and numerical results, the latter suggesting that the answer I seek is $\sim x/\log x$. I'm interested in any advancements towards a solution.
##Background##
The prime counting function $\pi(x)$ can be written on the form $$ \pi(x) = \sum_{ \substack{ {n\leq x}\\{(n,P(\sqrt{n}))=1}}} 1-1 \sim \frac{x}{\log x}. $$ We construct a random model with the same multiplicative structure as the primes in terms of $$ \pi_{rm}(x) = \sum_{ \substack{ {n\leq x}\\{(n+a,P(\sqrt{n}))=1}}} 1-1, $$ where $a$ initially is any random integer. Thus, $(n+a,P(\sqrt{n}))$ can take any value from $1$ to $P(\sqrt{n})$. The expected value of $\pi_{rm}(x)$ in this case is simply $$ \mathbf{E}\left[ \pi_{rm}(x) \right] = \sum_{n\leq x} W(\sqrt{n}) \sim 2 \operatorname{e}^{-\gamma} \frac{x}{\log x}, $$ where $W(x)= \prod_{p\leq x} \left(1-1/p\right)$, and the asymptotic equality follows from Merten's product theorem. The variance in this case satisfies $\operatorname{Var}(\pi_{rm}(x)) < \sum_{n\leq x} W(\sqrt{n})(1-W(\sqrt{n}))$ for $x\geq 2$.
Consider now the fact that the constraint $(n, P(\sqrt{n}))\leq n$ is satisfied for the primes. The prime counting function $\pi(x)$ therefore lies in a subspace of $\pi_{rm}(x)$ corresponding to those values of $a$ such that $(n+a,P(\sqrt{n}))\leq n$ for all $n\leq x$. This gives us the random model $\pi_{rm_c}(x)$ in the question.
Legendre sieve perspective
Can the Legendre sieve be a possible approach? Let $A(a) = \left\{ m: 1+a \leq m \leq x + a \right\}$ where $a$ is an integer such that $(n+a,P(\sqrt{n}))\leq n$ for all $n\leq x$. Also, let $A_d(a)$ be the set of integers in $A(a)$ divisible by $d$. In general, when $d\leq x$, $|A_d(a)|$ take either of the values $\lfloor x/d \rfloor=x/d -\{x/d\}$ or $\lfloor x/d \rfloor + 1 = x/d -\{x/d\}+1$. When $d>x$, $|A_d(a)| = \lfloor x/d \rfloor = 0$, meaning we only need to consider $d\leq x$. We then obtain the Legendre identity:
\begin{align} S(A(a),P(\sqrt{x})) &= \sum_{\substack{ {d\mid P(\sqrt{x})}\\ {d\leq x}}}\mu(d)|A_d(a)|\\ &= x \sum_{\substack{ {d\mid P(\sqrt{x})}\\ {d\leq x} } } \frac{\mu(d) }{d} - \sum_{\substack{ {d\mid P(\sqrt{x})}\\ {d\leq x} } } \mu(d) \left\{ \frac{x}{d}\right\} + \sum_{\substack{ {d\mid P(\sqrt{x}) }\\ {d\leq x} \\ {|\mathcal{A}_d(a)| = \lfloor x/d \rfloor+1 } }} \mu(d). \end{align}
In the case of the primes, $a=0$ and $|A_d(0)| = \lfloor x/d \rfloor$. The last term in the previous equation becomes zero and we obtain
\begin{align} \pi(x) - \pi(\sqrt{x}) + 1 &= S(A(0),P(\sqrt{x})) \\ &= x \sum_{\substack{ {d\mid P(\sqrt{x})}\\ {d\leq x} } } \frac{\mu(d) }{d} - \sum_{\substack{ {d\mid P(\sqrt{x})}\\ {d\leq x} } } \mu(d) \left\{ \frac{x}{d}\right\}. \end{align}
From the prime number theorem the left hand side of this equation is $\sim x/\log x$. Also, as shown by @Lucia in the MO post Asymptotic limit of truncated Legendre sieve, we have that
$$ x \sum_{\substack{ {d\mid P(\sqrt{x})}\\ {d\leq x} } } \frac{\mu(d) }{d} \sim \frac{x}{\log x}. $$
It therefore follows that
$$ \sum_{\substack{ {d\mid P(\sqrt{x})}\\ {d\leq x} } } \mu(d) \left\{ \frac{x}{d}\right\} = o\left(\frac{x}{\log x}\right). $$
To obtain an asymptotic estimate of $\mathbf{E}[S(A(a),P(\sqrt{x}))]$ one therefore needs to evaluate or bound the expected value of
$$ \sum_{\substack{ {d\mid P(\sqrt{x}) }\\ {d\leq x} \\ {|\mathcal{A}_d(a)| = \lfloor x/d \rfloor+1 } }} \mu(d). $$
Numerical results
Consider the random models with and without the constraint on $a$ for $x=p_{41}^2-1$. For this value of $x$ the sample space of the unconstrained model contains more than $1.6 \times 10^{68}$ elements, while the sample space of the constrained model contains only 88 elements.
In Figure A we see 88 realisations of each of the two error terms $\pi_{rm}(x) - \mathbf{E}[\pi_{rm}(x)]$ (dark gray) and $\pi_{rm_c}(x) - \mathbf{E}[\pi_{rm}(x)]$ (light gray). The black line shows $\operatorname{li}(x)-\mathbf{E}[\pi_{rm}(x)]$.
In Figure B we see the 88 realisations of the error term $\pi_{rm_c}(x) - \operatorname{li}(x)$ (light gray). The mean of the 88 realisations is displayed as dark gray. Also shown are $\pi(x) - \operatorname{li}(x)$ (black) and $R(x) - \operatorname{li}(x)$ (white), where $R(x)$ is the Riemann prime counting function.
What seems to be the case is that the constraint on $a$ forces all elements in the constrained random model to be strongly correlated. This suggests that perhaps not only the expected value of the constrained random model is $\sim x/\log x$, but that all elements in this model has the same asymptotic mean.