Evaluate the sum $\sum_{m,n\geq 1}\frac{1}{m^2n+n^2m+kmn} $ Motivated by the double summation posted here and recently here too. I came up with the general   closed form for one parameter $k>0$ where $k$ being a positive integer.

$$\sum_{m=1}^{\infty}\sum_{n=1}^{\infty} \frac{1}{m^2n+n^2m+kmn}=\frac{H^2_{k}-\psi_1(k+1)}{k}+\frac{\zeta(2)}{k}$$ where $H_n$ is nth Harmonic number and $\psi_1(x)$ is trigamma  function.

How do I prove  the above equality is true?

Note that $$\frac{1}{m^2n+n^2m+kmn}=\frac{1}{mn}\int_0^1 x^{m+n+k-1}dx $$ and thus $$\sum_{m,n\geq 1}\frac{1}{mn}\int_0^1 x^{k-1}\cdot x^{m+n} dx=\int_0^1 x^{k-1} \sum_{m,n\geq 1}\frac{x^{m+n}}{mn}dx=\int_0^1x^{k-1} \ln^2(1-x)dx \cdots(1) $$ where we  exploit  the series of  $\ln(1-x)$.
Now by  definition of beta function  we have $$B(k,y)=\int_0^1 x^{k-1} (1-x)^{y-1}dx\cdots(2)$$  having second derivatives  of  $(2)$ at $y=1$ we obtained the right hand expression of $(1)$.  That $$\lim_{y\to 1^{+}}B(k,y)=\int_0^1 x^{k-1} \ln^2(1-x)dx$$ We evaluate  left hand side of $(2)$ $$\lim_{y\to 1^+}\frac{\partial^2}{\partial y^2} B(k,y)= \lim_{y\to 1^+}\frac{\partial }{\partial y}B(k,y)\left(\psi_0(y)-\psi_0(k+y)\right)$$ on further differentiation  and setting $y=1$ we have $$ B(k,1)\left((-\psi_0(k+1)-\gamma)^2-\psi_1(k+1)+\frac{\pi^2}{6}\right)=\frac{H_k^2-\psi_1(k+1)}{k}+\frac{\zeta(2)}{k}$$ now  we have desired closed form as

$$\frac{H_{k+1}^2-\psi_1(k+2)}{k+1}+\frac{\zeta(2)}{k+1}=\frac{H_k^2+H_{k}^{(2)}}{k}$$


Alternatively
We show  that $$ \int_0^1 x^{k-1} \ln^2(1-x)dx =\frac{{H_k^2}+H_k^{(2)}}{k}$$
Using the integration trick  used in the book, (Almost) Impossible integrals, sums and series by Sir Cornel loan Vălean, page no 59-60.
Consider the integral  of the form  for all $n\geq  1$.$$I(n) = \displaystyle \int_0^1 x^{n-1} \ln(1-x) dx =\frac{1}{n}\displaystyle \int_0^1 \frac{d}{dx}(x^{n}-1)\ln(1-x) dx$$  and by integration by parts we yield
$$I(n)=\frac{1}{n}\left[\underbrace{(x^{n}-1)\ln(1-x)}_{0}\right]_0^1-\frac{1}{n}\int_0^1\frac{1-x^{n}}{1-x}dx $$ $$=-\frac{1}{n}\int_0^1\sum_{j=1}^{n} x^{j-1}
=-\frac{1}{n}\sum_{j=1}^{n}\int_0^1 x^{j-1}dx=-\frac{1}{n} \sum_{j=1}^n\frac{1}{j} =-\frac{H_n}{n}\cdots  (3)$$
Further, consider the  integral $ I(k)=\displaystyle \int_0^1 x^{k}\ln^2(1-x) dx$ for $ k\geq 0 $ which further can be written as $\displaystyle \frac{1}{k+1}\int_0^1 \frac{d}{dx}(x^{k+1}-1)\ln^2(1-x) dx$ and hence on Integration by parts we see that $$I(k)=\frac{1}{k+1}\left[\underbrace{(x^{k+1}-1)\ln^2(1-x)}_{0}\right]_0^1-\frac{2}{k+1}\int_0^1\frac{1-x^{k+1}}{1-x}\ln(1-x)dx $$ $$=-\frac{2}{k+1}\int_0^1\sum_{n=1}^{k+1} x^{n-1}\ln(1-x)=-\frac{2}{k+1}\sum_{n=1}^{k+1}\int_0^1 x^{n-1}\ln(1-x)dx$$ Plugging the result from $(3)$ to last integral we have $$\frac{2}{k+1}\sum_{n=1}^{k+1}\frac{H_n}{n}=\frac{2}{k+1}\left(\frac{H_{k+1}^2+H_{k+1}^{(2)}}{2}\right)=\frac{H_{k+1}^2 +H_{k+1}^{(2)}}{k+1}$$ Further note that nth partial sum of $\displaystyle H_{k+1}^{(2)} = \zeta(2) -\psi^1(k+2)$  giving us required result $$\frac{H_{k+1}^2-\psi^1(k+2)}{k+1}+\frac{\pi^2}{6(k+1)} $$ to get the desired  integral we let $k\to k-1$ which finally result us $$\frac{H_{k}^2+H_{k}^{(2)}}{k}=\frac{H_{k}^2-\psi_{1}(k+1)}{k}+\frac{\zeta(2)}{k}$$
