Reading about the abc conjecture and the interesting definition and properties of radicals of an integer, I wondered about the properties and bounds of the sum $$S_{\text{rad}}(k)=\sum_{n=1}^k \text{rad}(n)$$
Note that, based on the property that $\text{rad}(n)=n$ if $n$ is square free, and $\text{rad}(n)<n$ otherwise, we have that $$\sum_{n_{\text{square free}}\leq k} n <S_{\text{rad}}(k)<\sum_{n=1}^k n$$
We have that $\sum_{n=1}^k n = \frac{k^2+k}{2}$, and it is not difficult to show that $\sum_{n_{\text{square free}}\leq k} n \sim \frac{6}{\pi^2}\left(\frac{k^2+k}{2}\right)$; therefore, a first tentative bound could be $$\frac{6}{\pi^2}\left(\frac{k^2+k}{2}\right)<S_{\text{rad}}(k)<\left(\frac{k^2+k}{2}\right)$$
Now it comes the interesting part. I wondered if $S_{\text{rad}}$ could converge to some constant $\frac{6}{\pi^2}<C<1$, and that is indeed what numerical results seem to yield. Concretely, $C\approx 0.705...$ I have found that those are the first digits of $$\sum_{n=1}^\infty \frac{1}{p_n\#} = \frac{1}{2}+\frac{1}{2\times3}+\frac{1}{2\times3\times5}+\dots$$
where $p_n\#$ is the nth Primorial.
Here I have two graphs (forgive me the spanish legends) showing (i) the comparison between $S_{\text{rad}}(k)$ and $C \cdot \left(\frac{k^2+k}{2}\right)$ for every $k\leq 10^4$, and (ii) the numerical differences between them (that can not be appreciated in the first graph due to scale). It seems that the differences oscilate around $0$ (showing many values of $k$ for which $S_{\text{rad}}(k)=C \cdot \left(\frac{k^2+k}{2}\right)$).
And here comes my
Question. How can be proved / disproved that $C=\sum_{n=1}^\infty \frac{1}{p_n\#}$?


It turns out that the correct constant is $$ C = \prod_p \biggl( 1 - \frac1{p(p+1)} \biggr) \approx 0.70444, $$ which is already smaller than $\frac12+\frac16+\frac1{30}+\frac1{210}$. Here's a justification (with many steps missing).
When summing a multiplicative function, the main term is very often related to the residue of the rightmost pole of the corresponding Dirichlet series. For example, Perron's formula tells us that $$ \sum_{n\le x} n = \int_{c-i\infty}^{c+i\infty} \biggl( \sum_{n=1}^\infty n\cdot n^{-s} \biggr) \frac{x^s}s \, ds = \int_{c-i\infty}^{c+i\infty} \zeta(s-1) \frac{x^s}s \, ds, $$ where $c$ is any constant larger than $2$. The integrand has its rightmost pole at $s=2$, where its residue is $\frac{x^2}2$—this is the asymptotic main term of the left-hand side.
Similarly, $$ \sum_{\substack{n\le x \\ n\text{ squarefree}}} n = \int_{c-i\infty}^{c+i\infty} \biggl( \sum_{\substack{n\ge1 \\ n\text{ squarefree}}} n\cdot n^{-s} \biggr) \frac{x^s}s \, ds = \int_{c-i\infty}^{c+i\infty} \frac{\zeta(s-1)}{\zeta(2s-2)} \frac{x^s}s \, ds, $$ whose residue at $s=2$ is $\frac1{\zeta(2)} \frac{x^2}2 = \frac6{\pi^2}\frac{x^2}2$, again giving the main term.
So we can look at $$ \sum_{n\le x} \mathop{\rm rad}(n) = \int_{c-i\infty}^{c+i\infty} \biggl( \sum_{n=1}^\infty \mathop{\rm rad}(n) n^{-s} \biggr) \frac{x^s}s \, ds = \int_{c-i\infty}^{c+i\infty} \prod_p \biggl( 1 + \sum_{k=1}^\infty \frac p{p^{ks}} \biggr) \frac{x^s}s \, ds, $$ where we've converted the sum into its Euler product using multiplicativity. Some calculation reveals that $$ 1 + \sum_{k=1}^\infty \frac p{p^{ks}} = \frac{1-p^{-s}-p^{2-2 s}-p^{1-2 s}}{(1-p^{-(s-1)})(1-p^{-s})}, $$ and so $$ \prod_p \biggl( 1 + \sum_{k=1}^\infty \frac p{p^{ks}} \biggr) = \zeta(s-1) \prod_p \frac{1-p^{-s}-p^{2-2 s}+p^{1-2 s}}{1-p^{-s}} . $$ The residue of this function times $\frac{x^s}s$ at $s=2$ equals $$ \frac{x^2}2 \prod_p \frac{1-p^{-2}-p^{-2}+p^{-3}}{1-p^{-2}} = \frac{x^2}2 \prod_p \biggl( 1-\frac1{p(p+1)} \biggr). $$
(There is probably a way to get this main term without using complex analysis, and I suspect that it has already been done in the math literature.)