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$?

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\} $$