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.