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

1. Euler totient function: $$S_{\varphi}(x)=\sum_{i=1}^x\sum_{j=1}^x\varphi(i\cdot j)$$
2. Sum of divisors function: $$S_{\sigma_1}(x)=\sum_{i=1}^x\sum_{j=1}^x\sigma_1(i\cdot j)$$
3. 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.

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.
• Thank you very much for your answer. Different forms $f$ can take involve additional approaches for simplifying the sum, as can be shown from your $\phi$ example. Following on your excellent explanation, I am trying to extend this to the sum of divisors function, $\sigma(m\cdot n)$. Summing the divisors of different $n\leq x$ can be simplified greatly by running the sum over the divisors rather than the numbers themselves, with each divisor repeating $\lfloor \frac xd \rfloor$ times. This gives $\sum_{d=1}^x d\cdot \lfloor \frac xd \rfloor$. Can you show how this fits into your observation? – MC From Scratch May 31 '19 at 10:01
• The reason I am interested in this, is that $\sum_{d=1}^x d\cdot \lfloor \frac xd \rfloor$ has an $O(n^{1/2+\epsilon})$ computation, using Dirichlet convolution. – MC From Scratch May 31 '19 at 10:02
• The key starting step (very often, in my mind) is to find a way to write the multiplicative function $f(n)$ as $\sum_{d\mid n} g(d)$, then switch the order of summation to get $\sum_{n\le x} f(n) = \sum_{d\le x} g(d) \lfloor x/d\rfloor$. When $f=\sigma$ it's very easy, since $\sigma(n) = \sum_{d\mid n} g(d)$ with $g(d)=d$; that's how the formula from your comment is proved. When $f(n) = \sigma(mn)/\sigma(m)$, the corresponding formula is $f(n) = \sum_{d\mid n} g(d)$ with $g(d) = d\prod_{p^k\|m,\,p\mid d} (1-(p^k-1)/(p^{k+1}-1))$ (found by looking at prime power inputs). – Greg Martin May 31 '19 at 17:39
• Also, the first line of the sum gives a different result from the 2nd line and forth. I tried both paper and computer calculations (which were the same, but differed between the first and second lines of the sum). For example, when $f=\sigma$, $S(5)$ results in $401$ on the first line, but $210$ on the second/third/fourth lines. Might there be a mistake? – MC From Scratch May 31 '19 at 18:26