A similar problem to this problem (ccorn has given a nice answer to it). Prove that $$\sum_{\substack{0<k<3^n\\3\nmid k}}\sigma{(3^n-k)}\sigma{(k)}=6\cdot27^{n-1},$$ where $\sigma(N)$ is defined as $$\sigma{(N)}=\sum_{d\mid N}d.$$
I find this:(From Elementary Evaluation of Certain Convolution Sums Involving Divisor Functions by Huard et al. On page 20)
so we are done.
For $p$ prime and $n$ a nonnegative integer, define
$$f_p(n) = \sum_{\substack{k=1,\ldots,p^n-1\\p\not\mid k}} \sigma(k)\,\sigma(p^n - k)$$
Proposition:
$$f_p(n)=\begin{cases} 0 & \text{for }n=0 \\[1ex] \frac{1}{12}\,\frac{p^2 - 1}{p^2 + 1} \left(5 (p - 1)\,p^{3n-1} - (p - 2)(p - 3)(-p)^{n-1}\right) & \text{for }n \geq 1 \end{cases}$$ This confirms the previous result for $p=2$ and the claim for $p=3$ and generalizes the formula to arbitrary primes $p$. I actually had hoped for a formula that is as simple for $p>3$ as it is for $p=2$ and $p=3$, but as you can see, this is not the case.
To establish that result, I will use generating functions, along with a formula supplied by Ethan, which I will prove first.
As before, we will use the Ramanujan functions $P(q)$ and $Q(q)$ and the corresponding normalized Eisenstein series $\mathrm{E}_2(\tau)$ and $\mathrm{E}_4(\tau)$ via $q = \mathrm{e}^{2\pi\mathrm{i}\tau}$ where $\tau$ ranges over the complex upper half plane $\mathbb{H}$:
$$\begin{aligned} \mathrm{E}_2(\tau) &= P(q) = 1 - 24\sum_{n=1}^\infty \sigma(n)\,q^n \\ \mathrm{E}_4(\tau) &= Q(q) = 1 + 240\sum_{n=1}^\infty \sigma_3(n)\,q^n \end{aligned}$$
When looking at the formula supplied by Ethan,
$$\sum_{k=1,\ldots,n-1} \sigma(k)\,\sigma(n - k) = \frac{1}{12}\left(5\sigma_3(n) - (6n - 1)\,\sigma(n)\right)$$ it is helpful to think about generating functions. You will recognize the $\sigma_3$ as coming from $Q$ and the $\sigma$ as coming from $P$. Finally, the $n$ in the factor $(6n-1)$ comes from... differentiation. Yes.
In terms of generating functions, the formula is actually quite famous:
$$P^2 = Q + 12q\frac{\mathrm{d}P}{\mathrm{d}q}$$ It is part of a triplet of equations that relate Ramanujan's functions $P$, $Q$, $R$ with their derivatives. But how can we prove it, elegantly?
Key facts are again the modular symmetries of the Eisenstein series. Let $\tau' = \frac{-1}{\tau}$, then
$$\begin{aligned} \mathrm{E}_2(\tau+1) &= \mathrm{E}_2(\tau) & \mathrm{E}_2(\tau') &= \frac{6\tau}{\pi\,\mathrm{i}} + \tau^2\,\mathrm{E}_2(\tau) \\ \mathrm{E}_4(\tau+1) &= \mathrm{E}_4(\tau) & \mathrm{E}_4(\tau') &= \tau^4\,\mathrm{E}_4(\tau) \end{aligned}$$ These symmetries are demonstrated in introductory textbooks such as
Differentiating the formula for $\mathrm{E}_2(\tau')$ with respect to $\tau$, you will find
$$\mathrm{E}_2'(\tau') = \frac{6\tau^2}{\pi\,\mathrm{i}} + 2\tau^3\,\mathrm{E}_2(\tau) + \tau^4\,\mathrm{E}_2'(\tau)$$ Now let
$$L = \mathrm{E}_2^2 - \frac{6}{\pi\mathrm{i}}\,\mathrm{E}_2'$$ and examine what happens to $L$ under the transformations $\tau\to\tau+1$ and $\tau\to\tau'$. You will find:
$$\begin{aligned} L(\tau+1) &= L(\tau) \\ L(\tau') &= \left(\frac{6\tau}{\pi\mathrm{i}} + \tau^2\mathrm{E}_2(\tau)\right)^2 - \frac{6}{\pi\mathrm{i}} \left(\frac{6\tau^2}{\pi\,\mathrm{i}} + 2\tau^3\,\mathrm{E}_2(\tau) + \tau^4\,\mathrm{E}_2'(\tau)\right) \\ &= \tau^4\left(\mathrm{E}_2^2(\tau) - \frac{6}{\pi\mathrm{i}}\,\mathrm{E}_2'(\tau)\right) \\ &= \tau^4\,L(\tau) \end{aligned}$$ That is, $L$ has the same modular symmetries as $\mathrm{E}_4$. Furthermore, from its definition in terms of $\mathrm{E}_2$, we see that $L$ is holomorphic over $\mathbb{H}$ and has a Maclaurin series in $q$. Thus, $L$ is an entire modular form of weight $4$.
As Apostol shows in his modular functions text (chapter 6, section 5), there exists, up to a constant factor, just one entire modular form of weight $4$, and that is $\mathrm{E}_4$. In fact, by comparing limits for $\Im\tau\to\infty$ (thus forcing $q\to0$), we find that the factor is $1$, so $L = \mathrm{E}_4$.
There it is again: Finding modular symmetries with low weights can suffice to identify functions. You must make sure however that the function you examine has no poles, not even for $\Im\tau\to\infty$ where $q\to0$. This necessitates the above checks about being holomorphic and having a Maclaurin series. However, all this turns out to be amazingly easy.
Knowing that $L=\mathrm{E}_4$, we can now express the derivative of $\mathrm{E}_2$ in terms of $\mathrm{E}_2$ and $\mathrm{E}_4$. Thus we find
$$\mathrm{E}_2^2 = \mathrm{E}_4 + \frac{6}{\pi\mathrm{i}}\,\mathrm{E}_2'$$ Note that $q\frac{\mathrm{d}}{\mathrm{d}q} = \frac{1}{2\pi\mathrm{i}}\frac{\mathrm{d}}{\mathrm{d}\tau}$. This allows us to translate $\frac{6}{\pi\mathrm{i}}\mathrm{E}_2'$ to $12q\frac{\mathrm{d}P}{\mathrm{d}q}$. Altogether,
$$P^2 = Q + 12q\frac{\mathrm{d}P}{\mathrm{d}q}$$ which is what I wanted to prove.
We will now express the above result using $q$-series. For positive integer $m$, define
$$s(m) = \sum_{k=1,\ldots,m-1} \sigma(k)\,\sigma(m - k)$$ Its generating $q$-series is
$$\begin{aligned} \sum_{m=1}^\infty s(m)\,q^m &= \left(\frac{1-P(q)}{24}\right)^2 \\ &= \frac{1 - 2P(q) + P^2(q)}{24^2} = \frac{1 - 2P(q) + Q(q) + 12q\frac{\mathrm{d}}{\mathrm{d}q}P(q)}{24^2} \\ &= \frac{1}{12} \sum_{m=1}^\infty \left((1 - 6m)\,\sigma(m) + 5\sigma_3(m)\right) q^m \end{aligned}$$ Comparing coefficients yields the formula supplied by Ethan.
The modular stuff ends here. Below, I will merely use generating functions such as $\frac{1}{1-X} = 1 + X + X^2 + \cdots$ which should be familiar.
Now we focus on those $m$ that we are interested in. For prime $p$ and nonnegative integer $n$, let
$$\begin{aligned} h_p(n) = s(p^n) &= \frac {5 (p^{3(n+1)} - 1) - (p^2+p+1)(6 p^n - 1)(p^{n+1} - 1)} {12 (p^3 - 1)} \\ &= \frac {(5p^3)\,p^{3n} - \left(6 p\,(p^2+p+1)\right) p^{2n} + \left((p+6)(p^2+p+1)\right) p^n - (p^2+p+6)} {12 (p^3 - 1)} \end{aligned}$$ (To ease things, separate factors $p^{3n}, p^{2n}, p^n$ from factors with constant exponents.)
It is time to note a relationship of $h_p$ with the $f_p$ we are looking for:
$$\begin{aligned} h_p(n) = s(p^n) &= \sum_{k=1,\ldots,p^n-1} \sigma(k)\,\sigma(p^n - k) \\ &= \sum_{j=0}^{n-1} \sum_{\substack{ k=1,\ldots,p^{n-j}-1\\p\not\mid k}} \sigma(k\,p^j)\,\sigma(p^n - k\,p^j) \\ &= \sum_{j=0}^{n-1} \sigma(p^j)^2 \sum_{\substack{ k=1,\ldots,p^{n-j}-1\\p\not\mid k}} \sigma(k)\,\sigma(p^{n-j} - k) \\ &= \sum_{j=0}^{n-1} \sigma(p^j)^2\,f_p(n-j) = \sum_{j=0}^{n} \sigma(p^j)^2\, \underbrace{f_p(n-j)}_{0\text{ for }j=n} \end{aligned}$$ In pulling out $\sigma(p^j)^2$, we use the fact that $\sigma$ is multiplicative for coprime factors of its arguments.
Consequently, when using generating functions
$$\begin{aligned} F_p(X) &= \sum_{n=0}^\infty f_p(n)\,X^n \\ G_p(X) &= \sum_{n=0}^\infty \sigma(p^n)^2\,X^n \\ H_p(X) &= \sum_{n=0}^\infty h_p(n)\,X^n \end{aligned}$$ we find
$$\begin{aligned} H_p(X) &= G_p(X)\,F_p(X) \end{aligned}$$ Henceforth, expressions with $X$ shall implicitly be understood as the corresponding formal power series. We will not evaluate the series at any particular $X$, so we do not need to address issues of domain and convergence.
The generating function $H_p(X)$ is
$$\begin{aligned} H_p(X) &= \sum_{n=0}^\infty h_p(n)\,X^n \\ &= \frac{1}{12 (p^3 - 1)}\left( \frac{5p^3}{1 - p^3 X} - \frac{6\,p(p^2+p+1)}{1 - p^2 X} + \frac{(p+6)(p^2+p+1)}{1 - p X} - \frac{(p^2+p+6)}{1 - X}\right) \\ &= \frac{(p^2 - 1)}{12}\,\frac{(5 p - 6) X + p^3 X^2} {(1 - p^3 X)(1 - p^2 X)(1 - p X)(1 - X)} \end{aligned}$$
The generating function $G_p(X)$ is
$$\begin{aligned} G_p(X) &= \sum_{n=0}^\infty \sigma(p^n)^2\,X^n \\ &= \frac{1}{(p - 1)^2}\sum_{n=0}^\infty \left(p^{n+1} - 1\right)^2\,X^n \\ &= \frac{1}{(p - 1)^2}\sum_{n=0}^\infty \left(p^2\,p^{2n} -2p\,p^n + 1\right)\,X^n \\ &= \frac{1}{(p - 1)^2} \left(\frac{p^2}{1 - p^2 X} - \frac{2p}{1 - p X} + \frac{1}{1 - X}\right) \\ &= \frac{1 + p X}{(1 - p^2 X)(1 - p X)(1 - X)} \end{aligned}$$
Therefore
$$\begin{aligned} F_p(X) &= \frac{H_p(X)}{G_p(X)} = \frac{p^2 - 1}{12} \,\frac{(5 p - 6) X + p^3 X^2}{(1 - p^3 X)(1 + p X)} \\ &= \frac{1}{12}\,\frac{p^2 - 1}{p^2 + 1} \left(\frac{p^2}{1 - p^3 X} + \frac{1}{1 + p X}\right) \left((5 p - 6) X + p^3 X^2\right) \\ &= \frac{1}{12}\,\frac{p^2 - 1}{p^2 + 1} \Biggl((p^2 + 1)(5 p - 6) X + \\ &\qquad\qquad\qquad \sum_{n=2}^\infty \left((5 p - 6) \left(p^{3(n-1)+2} + (-p)^{n-1}\right) + p^3 \left(p^{3(n-2)+2} + (-p)^{n-2}\right)\right) X^n\Biggr) \\ &= \frac{1}{12}\,\frac{p^2 - 1}{p^2 + 1} \left((p^2 + 1)(5 p - 6) X + \sum_{n=2}^\infty \underbrace{\left(5 (p - 1)\,p^{3n-1} - (p - 2)(p - 3)(-p)^{n-1}\right)}_{=(p^2+1)(5p-6)\text{ for }n=1} X^n\right) \\ &= \frac{1}{12}\,\frac{p^2 - 1}{p^2 + 1} \sum_{n=1}^\infty \left(5 (p - 1)\,p^{3n-1} - (p - 2)(p - 3) (-p)^{n-1}\right) X^n \end{aligned}$$ This proves the proposition. Let me point out that the modular part was almost effortless in comparison.
