I am interested in the function $f_n(m)$ which can be defined by the Dirichlet generating function $$\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(m)$:

$$\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(m) $ is known due to Walfisz (Chapter III in Weylsche Exponentialsummen in der neueren Zahlentheorie):

$$\sum\limits_{m \leqslant x} \sigma_1 (m) = \frac{\pi^2}{12} x^2 + O(x\log^{2/3}x) \tag 2$$

Since $\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)$ using $(2).$

So far, I have tried estimating $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)$ to obtain a better error term in $(1)$?