An identity involving $[\sigma(n)]^2$ 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?
 A: Here's a direct-ish proof. at least, without using that the various functions are multiplicative.
First, let's prove a pair of counting lemmas.
Lemma 1: If $m\mid n^2$, then the number of ways to write $m=d_1d_2$ with $d_1\mid n$ and $d_2\mid n$ is equal to the number of divisors of $\gcd\left(m,\frac{n^2}{m}\right)$.
Proof: We need both that $d_1\mid m$ and $d_1\mid n$. So we need $d_1\mid \gcd(m,n)$. We also need $d_2=\frac{m}{d_1}\mid n$, so $\frac{m}{d_1}\mid\gcd(m,n)$ or $\frac{m}{\gcd(m,n)}\mid d_1 $. Therefore, we can choose any $d_1$ with:
$$\frac{m}{\gcd(m,n)}\mid d_1\mid \gcd(m,n)$$
If we let $d=\frac{d_1\gcd(m,n)}{m}$, we can pick any $d$ with $$d\mid\frac{\gcd(n,m)^2}{m}=\gcd\left(\frac{n^2}{m},m\right).$$
Lemma 2: If $m\mid n^2$, then the number of ways to write $m=\frac{n}{d_1}d_2$ with $d_1\mid n$ and $d_2\mid d_1^2$ is also equal to the number of divisors of $\gcd\left(\frac{n^2}{m},m\right)$.
Proof: Note that $\gcd\left(\frac{n^2}{m},m\right)\mid n$ (this follows from $\gcd(m^2,n^2)=\gcd(m,n)^2\mid mn$.)
So, if $d\mid \gcd\left(\frac{n^2}{m},m\right)$ then $d\mid n$ so let $d_1=\frac nd$ and let $d_2=\frac{md_1}{n}=\frac{m}{d}$, which is an integer. 
Now, $\frac{d_1^2}{d_2} = \frac{n^2}{d^2}\frac{d}{m} = \frac{n^2}{md}$. But $d$ is a divisor of $\frac{n^2}m$ by definition, so this is an integer, and $d_2\mid d_1^2$.
On the other hand, if we have a pair $d_1,d_2$ with $\frac{n}{d_1}d_2=m$ and $d_1\mid n$ and $d_2\mid d_1^2$, then let $d=\frac{n}{d_1}$. We clearly have $d\mid n$.
Writing $\frac{n^2}{m}=\frac{n^2}{dd_2} = d\frac{d_1^2}{d_2}$ which is divisible by $d$.
Theorem: Let $\sigma_k(n)=\sum_{d\mid n} d^k$. Then:
$$\left[\sigma_k(n)\right]^2 = \sum_{d\mid n} \left(\frac{n}{d}\right)^k\sigma_k(d^2)$$
Proof: Both sides can be written as a sum of $k$th powers of divisors $m$ of $n^2$. The number of times $m^k$ occurs on the left side is counted in Lemma 1. The number of times $m$ occurs in the sum on the right side is counted by Lemma 2. These two values are the same.

It still feels almost accidental.
A: I don't know if this is what you're looking for but I believe there is a relatively simple way of seeing that 
$$(1 + p + ... + p^k)^2 = \displaystyle\sum_{m=0}^k p^{k-m}(1 + p + ... + p^{2m})$$
Namely, let us consider the coefficient of $p^j$ on both sides of the above equation. 
On the left hand side, this is clearly equal to the number of ways of writing $p^j = p^rp^s = p^{r+s}$ with $0 \leq r,s \leq k$, so it's the number of ways of writing $j = r+s$ with $0 \leq r,s \leq k$. Now if $j \leq k$, then this is possible for $r = 0,...,j$ and if $j \geq k$, it is possible for $r = j-k,...,k$.
On the right hand side, this coefficient is the number of $m$ such that the above sum contains $p^j$, which is clearly the number of $m$ such that $k-m \leq j \leq k+m$, with $0 \leq m \leq k$. Again, if $j \leq k$, then this is true
for $m = k-j,...,k$ and if $j \geq k$, it is true for $m = j-k,...,k$.
Thus the coefficient of $p^j$ is the same on both sides so they must be equal.
A: Another way to show the equality below that @jwsiegel proved:

$$(1 + p + \ldots + p^k)^2 = \sum_{m=0}^k p^{k-m}(1 + p + \ldots + p^{2m})$$

Group the terms in two different ways. (I think this is a nice direct visual proof.)
$$
\begin{array}{r}
\big(1 &+p &+p^2 &+\ldots &+p^k\big)\\
&+\big(p &+p^2 &+\ldots &+p^k &+p^{k+1}\big)\\
& &+\big(p^2 &+\ldots &+p^k &+p^{k+1} &+p^{k+2}\big)\\
& & & &\cdots\\
& & & &+\big(p^k &+p^{k+1} &+p^{k+2} &+\ldots &+p^{2k-1} &+p^{2k}\big)\\
\hline
& & & &=\\
& & & &p^k\\
& & &p^{k-1} &+p^k &+p^{k+1}\\
& & & &\cdots\\
& &p^2 &+\ldots+p^{k-1} &+p^k &+p^{k+1} &+p^{k+2} &+\ldots+p^{2k-2}\\
&p &+p^2 &+\ldots+p^{k-1} &+p^k &+p^{k+1} &+p^{k+2} &+\ldots+p^{2k-2} &+p^{2k-1}\\
1 &+p &+p^2 &+\ldots+p^{k-1} &+p^k &+p^{k+1} &+p^{k+2} &+\ldots+p^{2k-2} &+p^{2k-1} &+p^{2k}\\
\ \\
\end{array}
$$
