The number of representations of an integer as the inner product of integral lattice points I was looking through some old notes of mine and I came across a couple lemmas/identities I wrote down in regards to a question I asked about four years ago. In particular I wrote that for an arbitrary fixed integer $k>1$ we have the following asymptotic expansion:
$$\Psi_k(n)=\left|\{(\mathbf{u},\mathbf{v})\in \mathbb{N}^k\times \mathbb{N}^k:\mathbf{u}\cdot\mathbf{v}=n\}\right|\sim \frac{\sigma_{k-1}(n)\text{log}(n)^{k}}{\zeta(k)(k-1)!}$$
With $\sigma_{k-1}(n)=\sum_{d\mid n}d^{k-1}$ and $\zeta(k)=\sum_{n=1}^{\infty}\frac{1}{n^k}$ the Riemann zeta function. Now alternatively if one expands the inner product of $\mathbf{u}\cdot\mathbf{v}$ we see $\Psi_k$ can be expressed in a number of other ways:
$$\Psi_k(n)=\sum_{\substack{\mathbf{u}\cdot\mathbf{v}=n\\(\mathbf{u},\mathbf{v})\in \mathbb{N}^k\times \mathbb{N}^k}}1=\left|\{(u_1,v_1,\ldots u_k,v_k)\in \mathbb{N}^{2k}:n=u_1v_1+\cdots +u_kv_k\}\right|\\=\sum_{\substack{m_1+m_2+m_3+\cdots +m_k=n\\(m_1,m_2,m_3,\ldots ,m_k)\in \mathbb{N}^k}}d(m_1)d(m_2)d(m_3)\cdots d(m_k)$$
Where $d(m)=\sum_{d\mid m}1=\sigma_{0}(m)$ counts the divisors of any natural number $m$. Which unearths a number of $q$-series esque representation for the ordinary generating function of $\Psi_k$, for example:
$$\sum_{n=1}^{\infty}\Psi_k(n)q^n=\left(\sum_{n=1}^{\infty}q^{n^2}\frac{1+q^n}{1-q^n}\right)^k=\frac{\text{log}(1-q)^k}{\text{log}(q)^k}\sum_{j=0}^k\binom{k}{j}\left(\frac{\psi_{q}(1)}{\text{log}(1-q)}\right)^j$$
Where $\psi_{q}(z)=\frac{1}{\Gamma_{q}(z)}\frac{d}{dz}\Gamma_{q}(z)$ is the $q$-analog of the digamma function defined analogously in terms of the $q$-gamma function, expressible as $\Gamma_{q}(z)=(1-q)^{1-z}\prod_{n=0}^{\infty}\frac{1-q^{n+1}}{1-q^{n+z}}$ for $|q|<1$.

Now working with just some of these alternate representations, as well as fiddling with the order of the summands involved I was able to prove by induction on the integer $k>1$ that we have both:
$$\sum_{n\leq N}\Psi_k(n)=\frac{N^k\text{log}(N)^{k}}{k!}+\mathcal{O}(N^k\text{log}(N)^{k-1})$$
$$\sum_{n\leq N}\frac{\sigma_{k-1}(n)\text{log}(n)^{k}}{\zeta(k)(k-1)!}=\frac{N^k\text{log}(N)^{k}}{k!}+\mathcal{O}(N^k\text{log}(N)^{k-1})\\$$
So using the same heuristics as before it seems reasonable that $\Psi_k(n)\sim  \frac{\sigma_{k-1}(n)\text{log}(n)^{k}}{\zeta(k)(k-1)!}$ which at least for the case at $k=2$ would coincide with the answer to my previous question. However I'm unable to find a concrete proof of this result and would therefore appreciate any help in the matter.

 A: The asymptotic formula is true for even dimensions $k\geq 2$. We can prove this by induction on $k$, inspired by Rodrigo's observation on Eisenstein series. 
The case $k=2$ is classical and addressed in the OP's previous post that he linked. Now it suffices to show that if the formula is true for even dimensions $k,\ell\geq 2$, then it is also true for dimension $k+\ell$. Moreover, we can assume in this implication that either $k,\ell\geq 4$ or $k+\ell\in\{4,6\}$. Let us start from the obvious identity
$$ \Psi_{k+\ell}(n)=\sum_{r+s=n}\Psi_k(r)\Psi_\ell(s).$$
Using the fact that the logarithm is a slowly changing function, we are left with proving that
$$ \frac{\sigma_{k+\ell-1}(n)}{\zeta(k+\ell)\Gamma(k+\ell)}\sim\sum_{r+s=n}\frac{\sigma_{k-1}(r)}{\zeta(k)\Gamma(k)}\cdot\frac{\sigma_{\ell-1}(r)}{\zeta(\ell)\Gamma(\ell)}.$$
If $k,\ell\geq 4$, then the left hand side is $(2\pi i)^{-k-\ell}$ times the $n$-th Fourier coefficient of the standard holomorphic Eisenstein series $E_{k+\ell}$ of weight $k+\ell$ and full level, while the right hand side is $(2\pi i)^{-k-\ell}$ times the $n$-th Fourier coefficient of $E_kE_\ell$. As $E_{k+\ell}-E_kE_\ell$ is a cusp form, the result follows from standard (or even weaker) upper bounds for the Fourier coefficients of cusp forms. 
If $k+\ell\in\{4,6\}$, then the result follows from the explicit identities
$$5\sigma_3(n)+(1-6n)\sigma_1(n)=12\sum_{r+s=n}\sigma_1(r)\sigma_1(s),$$
$$21\sigma_5(n)+(10-30n)\sigma_3(n)-\sigma_1(n)=240\sum_{r+s=n}\sigma_1(r)\sigma_3(s).$$
These identities can also be proved with the help of Eisenstein series, or by elementary means, see e.g. (3.10) and (3.12) in the paper of Huard et al.
The proof is complete.
