Average decrease of a random number in "subtracting every prime by one in factorization"

Here’s a somewhat more pedestrian approach:

For a given prime $p$, a proportion $\frac{p-1}{p^{k+1}}$ of integers have exactly $k$ factors of $p$, and for these integers $a(n)$ is reduced by a factor $\left(\frac{p-1}p\right)^k$. Thus the average reduction is

\begin{eqnarray} \sum_{k=0}^\infty\frac{p-1}{p^{k+1}}\left(\frac{p-1}p\right)^k &=& \frac{p-1}p\sum_{k=0}^\infty\left(\frac{p-1}{p^2}\right)^k \\ &=& \frac{p-1}p\frac1{1-\frac{p-1}{p^2}} \\ &=& \frac{p^2-p}{p^2-p+1} \\ &=& \left(1+\frac1{p(p-1)}\right)^{-1}\;. \end{eqnarray}

Since $a(n)$ is completely multiplicative, the reductions due to different primes are multiplied, so the overall reduction factor, and thus the desired limit, is, as derived in Daniel Fischer’s answer, the reciprocal of Landau’s totient constant

$$ \prod_p\left(1+\frac1{p(p-1)}\right)=\frac{315}{2\pi^4}\zeta(3)\approx1.943596\;. $$

In "But the Dirichlet g.f. can also be written as" you have a spurious $1$, $$\sum_{n = 1}^{\infty} \frac{a(n)}{n^s} = \sum_{n = 1}^{\infty} A(n)\biggl(\frac{1}{n^s} - \frac{1}{(n+1)^s}\biggr) = s \int_1^{\infty} \frac{A(x)}{x^{s+1}}\,dx\,.$$ That doesn't affect the asymptotics of course.

If we look at the Euler product, we see that $$F(s) = \prod_p \frac{1}{1 - \frac{p-1}{p^s}}$$ is very close to $\zeta(s-1)$, hence it is promising to write $F(s) = \zeta(s-1)\cdot H(s)$, where $$H(s) = \prod_p \frac{1 - \frac{p}{p^s}}{1 - \frac{p-1}{p^s}} = \prod_p\frac{1 - \frac{p-1}{p^s} - \frac{1}{p^s}}{1 - \frac{p-1}{p^s}} = \prod_p\biggl(1 - \frac{1}{p^s - p + 1}\biggr)\,.$$ This product converges absolutely for $\operatorname{Re} s > 1$ and hence the principal part of $F(s)$ at the pole $s = 2$ is $$\frac{H(2)}{s-2}\,.$$ Now $A(x) = \frac{C}{2}x^2 + O(x^{2 - \varepsilon})$ implies that the principal part of $F(s)$ at $s = 2$ is the principal part of $$\frac{Cs}{2}\int_1^{\infty} \frac{dx}{x^{s-1}} = \frac{Cs}{2(s-2)} = \frac{C}{s-2} + \frac{C}{2}$$ at $s = 2$, i.e. $C = H(2)$.

We can evaluate $H(2)$ in closed form: \begin{align} H(2) &= \prod_p \biggl(1 - \frac{1}{p^2-p+1}\biggr) \\ &= \prod_p \frac{p(p-1)}{p^2-p+1} \\ &= \prod_p \frac{p(p-1)(p+1)}{p^3+1} \\ &= \prod_p \frac{p(p^2-1)(p^3-1)}{p^6-1} \\ &= \prod_p \frac{(1 - p^{-2})(1 - p^{-3})}{1 - p^{-6}} \\ &= \frac{\zeta(6)}{\zeta(2)\zeta(3)} \\ &= \frac{2\pi^4}{315\zeta(3)} \\ &\approx 0.5145\,. \end{align} We have not proved that indeed $A(x) = \frac{C}{2}x^2 + O(x^{2-\varepsilon})$, hence we have not yet proved that $$\lim_{n \to \infty} \frac{A(n)}{n^2/2} = H(2)\,,$$ and that thus the average order of $a(n)$ is $H(2)\cdot n$. This however follows either by Tauberian theorems ($A$ is monotonic, so a variant of the Wiener–Ikehara theorem is applicable) or by estimates for $\zeta(s-1)$ near the line $\operatorname{Re} s = 2$ and an argument similar to Landau's proof of the prime number theorem. Since $F(s)$ is, except for the pole at $s = 2$, holomorphic in the half-plane $\operatorname{Re} s > 1$ and decent estimates for $\lvert\zeta(s)\rvert$ in the critical strip are known, we can in this way establish $A(x) = \frac{1}{2}H(2)x^2 + O(x^{2-\varepsilon})$ for some $\varepsilon > 0$. Getting the remainder term to $O(x^{1+\varepsilon})$ or possibly even $O(x\log x)$ requires more careful estimates.