What is the density of squarefree numbers in $p+n$ for prime $p$ and natural $n$? Let $p_n$ be the $n$-th prime and $d$ be any positive integer. I was investigating for which $d$ is the density of square-free numbers in the sequence $p_n + d$ the highest or the lowest. Based on numerical data, I observed that $p_n+1$ had the lowest density of about $0.374$, while the primorials $d = 2, 6, 30, 210, 2310, \ldots$ had higher densities than any smaller value of $d$. For the primorial $6469693230$, which is the product of the first $10$ primes, the density of square-free numbers in the sequence $p_n + 6469693230$ is about $0.993$.
Question: Let $q_k$ be the product of the first $k$ primes. Is it true that the density of square-free numbers in the sequence $p_n + q_k$ strictly increases with $k$ and is higher than that of any other sequence $p_n + d$ where $0 \le d < q_k$?
 A: Let $S(x,N)$ denote the number of primes $p\le x$ such that $p+N$ is squarefree, then
\begin{aligned}
S(x,N)
&=\sum_{p\le x}\sum_{d^2|(p+N)}\mu(d)=\sum_{d\le\sqrt{x+N}}\mu(d)\sum_{\substack{p\le x\\p\equiv-N(d^2)}}1 \\
&=\sum_{d\le\sqrt{x+N}}\mu(d)\pi(x;d^2,-N)
\end{aligned}
Now, pick $1\le Q\le\sqrt{x+N}$, so that
$$
S(x,N)=\sum_{d\le Q}+\sum_{Q<d\le\sqrt{x+N}}
$$
For the second sum, we can apply the rough upper bound $\pi(x;q,a)\ll x/q$ to obtain:
$$
\sum_{Q<d\le\sqrt{x+N}}\ll\sum_{d>Q}{x\over d^2}\ll\frac xQ
$$
For the part $q\le Q$, we first need to isolate the case where $(d,N)>1$:
$$
\sum_{d\le Q}=\sum_{\substack{d\le Q\\(d,N)=1}}\mu(d)\pi(x;d^2,-N)+\mathcal O(Q)
$$
In this situation, we invoke Siegel-Walfisz theorem, which states that for all $(a,q)=1$ and all $H>0$ there is
$$
\pi(x;q,a)={\operatorname{li}x\over\varphi(q)}+\mathcal O\left(x\over\log^Hx\right)
$$
where the O constant only depends on $A$ but neither $a$ nor $q$. Thus, we have
$$
\sum_{d\le Q}=\operatorname{li}x\color{blue}{\sum_{\substack{d\le Q\\(d,N)=1}}{\mu(d)\over\varphi(d^2)}}+\mathcal O(Qx\log^{-H}x)
$$
The remaining task for us is to handle the blue series. First, using the fact that $\sum_{q\le z}q/\varphi(q)\ll z$, we have
$$
\sum_{\substack{d>Q\\(d,N)=1}}{\mu(d)\over\varphi(d^2)}\ll\frac1Q
$$
and by the properties of Dirichlet series, we get
$$
\sum_{\substack{d\ge1\\(d,N)=1}}{\mu(d)\over\varphi(d^2)}=\prod_{p\nmid N}\left\{1-{1\over p(p-1)}\right\}
$$
Synthesizing what we have, we get
$$
S(x,N)=\operatorname{li}x\prod_{p\nmid N}\left\{1-{1\over p(p-1)}\right\}+\mathcal O\left({x\log x\over Q}+{Qx\over\log^Hx}\right)
$$
Finally, let $Q=(\log x)^{A+1}$ and $H=2A+1$, the formula above gets simplified into
$$
S(x,N)=\operatorname{li}x\prod_{p\nmid N}\left\{1-{1\over p(p-1)}\right\}+\mathcal O\left(x\over\log^Ax\right)
$$
This asymptotic formula suggests OP's observation is correct. The density formula is given by
$$
\rho_N=\lim_{x\to+\infty}{S(x,N)\over\pi(x-N)}=\prod_{p\nmid N}\left\{1-{1\over p(p-1)}\right\}
$$
