I am interested in the function ${f_n}\left( m \right)$$f_n(m)$ which can be defined by the Dirichlet generating function $$\zeta \left( s \right)\zeta \left( {s - 1} \right)...\zeta \left( {s - n + 1} \right) = \sum\limits_{m = 1}^\infty {\frac{{{f_n}\left( m \right)}}{{{m^s}}}} $$$$\zeta(s)\zeta(s - 1) \cdots\zeta(s - n + 1) = \sum\limits_{m = 1}^\infty \frac{f_n(m)}{m^s} $$
This is called the number-of-sublattices function. It is straightforward to get the standard estimate of $\sum\limits_{m \leqslant x}^{} {{f_n}\left( m \right)}$$\sum\limits_{m \leqslant x} f_n(m)$:
$$\sum\limits_{m \leqslant x}^{} {{f_n}\left( m \right)} = \frac{{{x^n}}}{n}\prod\limits_{\ell = 2}^{n - 1} {\zeta \left( \ell \right)} + O\left( {{x^{n - 1}}\log x} \right)\quad \quad \quad \left( 1 \right)$$$$\sum\limits_{m \leqslant x} f_n(m) = \frac{x^n}{n} \prod\limits_{\ell = 2}^{n - 1} \zeta(\ell) + O\left( x^{n - 1}\log x \right) \tag 1$$
I am looking for an improvement in the error term. The sharpest known estimate for $\sum\limits_{m \leqslant x}^{} {{\sigma _1}\left( m \right)}$$\sum\limits_{m \leqslant x} \sigma_1(m) $ is known due to Walfisz (Chapter III in Weylsche Exponentialsummen in der neueren Zahlentheorie):
$$\sum\limits_{m \leqslant x}^{} {{\sigma _1}\left( m \right)} = \frac{{{\pi ^2}}}{{12}}{x^2} + O\left( {x{{\log }^{2/3}}x} \right)\quad \quad \quad \left( 2 \right)$$$$\sum\limits_{m \leqslant x} \sigma_1 (m) = \frac{\pi^2}{12} x^2 + O(x\log^{2/3}x) \tag 2$$
Since $\zeta \left( s \right)\zeta \left( {s - 1} \right) = \sum\limits_{m = 1}^\infty {\frac{{{\sigma _1}\left( m \right)}}{{{m^s}}}}$$\zeta(s)\zeta(s - 1) = \sum\limits_{m = 1}^\infty \frac{\sigma_1(m)}{m^s}$, I have been wondering if it is possible to get a sharper bound for the error term of (1)$(1)$ using (2).$(2).$
So far, I have tried estimating ${f_3}\left( m \right) = \sum\limits_{d|m}^{} {{\sigma _1}\left( d \right){{\left( {\tfrac{m}{d}} \right)}^2}}$$f_3(m) = \sum\limits_{d\mid m} \sigma_1(d)\left( \tfrac{m}{d} \right)^2 $ but without success. I have tried summation interchange identities as well as the Dirichlet hyperbola method. I also tried to do some work involving the Fourier series representation of the fractional part, again unsuccessfully. Does anybody know of any theorems or techniques which could be use (2)$(2)$ to obtain a better error term in (1)$(1)$?