For a positive integer $n$, let $\sigma(n)$ denote the sum of the divisors of $n$. For example, $\sigma(1)=1$, $\sigma(2)=3$, $\sigma(4)=7$, etc.
I would like to prove the following identity:
For every positive integer $n$, $$ [\sigma(n)]^2=\sum_{d\mid n} \frac{n}{d} \sigma(d^2)$$ where the sum ranges over all divisors $d$ of $n$.
Here is one attempt: We know $\sigma(n)$ is a multiplicative function, that is, $\sigma(ab)=\sigma(a)\sigma(b)$ for $\gcd(a, b)=1$. And I think one can show the right hand side is also multiplicative, so it suffices to check the equality when $n=p^k$ is a prime power. This last task boils down to showing that: $$ (1+p+\cdots + p^k)^2 = \sum_{0\leq m\leq k} p^{k-m} (1+p+\cdots + p^{2m}) $$ Using the geometric series, the left hand side is $\left(\dfrac{p^{k+1}-1}{p-1}\right)^2$ and similarly the right hand side is $$ \sum_{0\leq m\leq k} p^{k-m} \frac{p^{2m+1}-1}{p-1} $$ So we just need to show: $$ (p^{k+1}-1)^2 = (p-1)\sum_{0\leq m\leq k} p^{k-m}(p^{2m+1}-1) = (p-1) p^{k} \sum_{0\leq m\leq k} (p^{m+1}-p^{-m}) $$ and indeed, we can compute that: $$ (p-1)p^{k} \sum_{0\leq m\leq k} (p^{m+1}-p^{-m}) = (p-1)p^{k}\left(\frac{p^{k+2}-p}{p-1} + \frac{p^{-k+1}-p}{p-1}\right) = p^{2k+2} -2 p^{k+1} + 1 = (p^{k+1}-1)^2 $$ as desired. Is this proof correct? I am not really satisfied by it, because I don't seem to get any insight out of it. Does anyone have a more conceptual (perhaps combinatorial) proof?