$\sum_{i=1}^x\sum_{j=1}^xf(i\cdot j)$ Double Summing a (Not Completely) Multiplicative Function Let $f(n)$ be a multiplicative function that is not completely multiplicative, i.e $f(m)\cdot f(n)= f(m\cdot n)$ only if $gcd(m,n)=1$. Let $S(x)$ be the double sum over $f$, that is:
$$S(x)=\sum_{i=1}^x\sum_{j=1}^xf(i\cdot j)$$
It is not difficult to see that if $f(n)$ were completely multiplicative, then $S(x)$ could be simplified: 
$$S(x)=\sum_{i=1}^x\sum_{j=1}^xf(i\cdot j)= \sum_{i=1}^xf(i)\sum_{j=1}^xf(j)= \biggl(\sum_{k=1}^xf(k)\biggr)^2$$
But since $f(n)$ is not completely multiplicative, this simplification is not completely true, and it fails in every combination where $gcd(i,j)\neq1$. Hence, $S(x)$ can be written this way provided we add some additional error term, let's call it $E$:
$$S(x)=\sum_{i=1}^x\sum_{j=1}^xf(i\cdot j)= \biggl(\sum_{k=1}^xf(k)\biggr)^2+E$$
$E$ is either negative or positive, I'm not sure. Obviously, $E$ is comprised of all the small errors generated by the initial sum term, when $gcd(i,j)\neq1$. I am mainly interested in the cases where $f(n)$ takes the form of:


*

*Euler totient function: $$S_{\varphi}(x)=\sum_{i=1}^x\sum_{j=1}^x\varphi(i\cdot j)$$

*Sum of divisors function: $$S_{\sigma_1}(x)=\sum_{i=1}^x\sum_{j=1}^x\sigma_1(i\cdot j)$$

*Moebius function: $$S_{\mu}(x)=\sum_{i=1}^x\sum_{j=1}^x\mu(i\cdot j)$$
My question is, what is this error term $E$ exactly? how can I calculate it? How can I properly sum all those small errors to get a correct evaluation of $S(x)$? For clarification, I am concerned with evaluating $S(x)$, but I think I must evaluate $E$ first in order to do it. I am taking this approach because I can compute $\biggl(\sum_{k=1}^xf(k)\biggr)^2$ very efficiently, and so, finding the error term $E$ will solve my question.
 A: One underappreciated but useful fact about multiplicative functions is the following: if $f(n)$ is multiplicative, and $k$ is any positive integer such that $f(k)\ne0$, then the function $g(n) = f(nk)/f(k)$ is a multiplicative function of $n$. (You'll get the proof correct on the first try.) In particular, we can write
\begin{align*}
\sum_{m\le x} \sum_{n\le x} f(mn) &= \sum_{m\le x} f(m) \sum_{n\le x} \frac{f(mn)}{f(m)}
\end{align*}
and use whatever techniques we want for sums of multiplicative functions on the inner sum. Using the Euler totient function as an example (quickly sketching the computation here):
\begin{align*}
\sum_{m\le x} \sum_{n\le x} \phi(mn) &= \sum_{m\le x} \phi(m) \sum_{n\le x} \frac{\phi(mn)}{\phi(m)} \\
&= \sum_{m\le x} \phi(m) \sum_{n\le x} \sum_{\substack{d\mid n \\ (d,m)=1}} \mu(d) \frac nd \\
&= \sum_{m\le x} \phi(m) \sum_{\substack{d\le x \\ (d,m)=1}} \frac{\mu(d)}d \sum_{\substack{n\le x \\ d\mid n}} n \\
&= \sum_{m\le x} \phi(m) \sum_{\substack{d\le x \\ (d,m)=1}} \mu(d) \sum_{m\le x/d} m \\
&\sim \sum_{m\le x} \phi(m) \sum_{\substack{d\le x \\ (d,m)=1}} \mu(d) \frac12 \bigg( \frac xd \bigg)^2 \\
&\sim \frac{x^2}2 \sum_{m\le x} \phi(m) \sum_{\substack{d\in\Bbb N \\ (d,m)=1}} \frac{\mu(d)}{d^2} \\
&= \frac{x^2}2 \frac1{\zeta(2)} \sum_{m\le x} \phi(m) \prod_{p\mid m} \bigg( 1-\frac1{p^2} \bigg)^{-1} \\
&= \frac{3x^2}{\pi^2} \sum_{m\le x} m \prod_{p\mid m} \frac p{p+1} \\
&\sim \frac{3x^2}{\pi^2} \frac{x^2}2 \prod_{p} \bigg( 1-\frac1p \bigg) \bigg( 1+ \frac p{p+1} \frac1p+ \frac p{p+1} \frac1{p^2} + \cdots \bigg) \\
&\sim \frac{3x^4}{2\pi^2} \prod_{p} \bigg( 1-\frac1{p(p+1)} \bigg) \approx 0.107062 x^4
\end{align*}
which is a good fit with experimental data.
