definite integralssummation
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?
We will look into the integral \begin{align*} I = &\int_{0}^{1}\arcsin^4 x\frac{ \ln x}{\sqrt{1-x^2}}\ \mathrm dx \end{align*} taking the @nospoon's novel approach presented here. Using the MacLaurin series of $\arcsin^4 x$ $$ \arcsin^4 x =\frac 3 2 \sum_{n=1}^\infty \frac{4^{n}H_{n-1}^{(2)}}{n^2{2n \choose n}}x^{2n} $$ and the fact that $$ \small\operatorname{B}(n+\tfrac 1 2,\tfrac 1 2) = \int_0^1 x^{n-1/2}(1-x)^{-1/2}\ \mathrm dx = 2\int_0^{\frac\pi 2} \sin^{2n}\theta\ \mathrm d\theta = \frac{\pi}{4^n}{2n \choose n},\tag{$\small x\mapsto \sin^2\theta$} $$ \begin{align*}\small \psi(n+\tfrac 12 ) -\psi(n+1) =&\small \sum_{k=1}^\infty \frac 1{\scriptsize k+n} - \frac 1{\scriptsize k+n-\tfrac 1 2} \\ =&\small\sum_{k=1}^\infty \left(\frac 1{\scriptsize k} - \frac 1{\scriptsize k-\tfrac 1 2}\right)-\sum_{k=1}^n\frac 1 {\scriptsize k} + \sum_{k=1}^n\frac 1{\scriptsize k-\tfrac 1 2}\\ =&\small-2\ln 2 -H_n +2(H_{2n}-\tfrac 1 2H_n)\\ =&\small 2(H_{2n}-H_n-\ln 2), \end{align*} \begin{align*} \Longrightarrow \ {\int_{ 0}^{1 }x^{2n}\frac{ \ln x}{\sqrt{1-x^2}}\ \mathrm dx} = & \frac 1 4\int_{0 }^{1 }x^{n-1/2} { \ln x \over \sqrt{1-x}}\ \mathrm dx\tag{$\small x^2\mapsto x$}\\ =& \frac 1 4 \left[\frac{\partial }{\partial x}\operatorname{B}(x,y) \right]_{x=n+1/2,y=1/2}\\ =&\frac 1 4\Big[ \operatorname{B}(x,y)\big[\psi(x) -\psi(x+y) \big]\Big]_{x=n+1/2,y=1/2}\\ =& \frac 1 4 \operatorname{B}(n+\tfrac 1 2,\tfrac 1 2)\big[\psi(n+\tfrac 12 ) -\psi(n+1) \big]\\ =& \frac{\pi}2\frac{{2n \choose n}}{4^{n}} \left(H_{2n} - H_n -\ln 2\right), \end{align*} where $\operatorname{B}(x,y)$ and $\psi(x)$ are the Beta and digamma function, respectively, we have \begin{align*} I = &\frac 3 2\sum_{n=1}^\infty \frac{4^{n}H_{n-1}^{(2)}}{n^2{2n \choose n}}\int_{0}^{1}x^{2n}\frac{ \ln x}{\sqrt{1-x^2}}\ \mathrm dx \\ =&\frac {3\pi}4 \sum_{n=1}^\infty \frac{H^{(2)}_{n-1}}{n^2}\left(H_{2n} - H_n -\ln 2\right) \\ =&\frac {3\pi}4\sum_{n=1}^\infty \frac{H^{(2)}_{n-1}H_{2n}}{n^2}-\frac {3\pi}4\underbrace{\sum_{n=1}^\infty \frac{H^{(2)}_{n-1}H_{n}}{n^2}}_{=-2\zeta(5) +2\zeta(2)\zeta(3)}-\frac {3\pi\ln 2}4\underbrace{\sum_{n=1}^\infty \frac{H^{(2)}_{n-1}}{n^2}}_{=\frac{3}4 \zeta(4)}\\ =&\frac{3\pi}{4} \sum_{n=1}^\infty \frac{H^{(2)}_{n}H_{2n}}{n^2} -\frac{3\pi}4\underbrace{\sum_{n=1}^\infty \frac{H_{2n}}{n^4}}_{=\frac{37}{4}\zeta(5)-4\zeta(2)\zeta(3)} +\frac{3\pi}2 \zeta(5) -\frac{\pi^3}4\zeta(3) -\frac{\pi^5\ln 2}{160}\\ =&\boxed{3\pi S -\frac{87\pi}{16} \zeta(5) +\frac{\pi^3}{4}\zeta(3) -\frac{\pi^5\ln 2}{160}} \end{align*} where $S = \sum_{n=1}^\infty \frac{H_{2n}H^{(2)}_{n}}{4n^2}$ is the sum in question, and the known values of several Euler sums $$ \sum_{n=1}^\infty \frac{H^{(2)}_{n-1}H_{n}}{n^2}=-2\zeta(5) +2\zeta(2)\zeta(3),\tag{1} $$ $$\sum_{n=1}^\infty \frac{H^{(2)}_{n}}{n^2}=\frac{7}4 \zeta(4),\tag{2} $$ \begin{align*}\sum_{n=1}^\infty \frac{H_{2n}}{n^4} =& 8\sum_{n=1}^\infty \frac{H_{n}}{n^4}-8\sum_{n=1}^\infty \frac{(-1)^{n-1} H_{n}}{n^4}\\ =&8\big(3\zeta(5)-\zeta(2)\zeta(3)\big)-8\left(\frac{59}{32}\zeta(5)-\frac 1 2\zeta(2)\zeta(3)\right)\\ =&\frac{37}4\zeta(5) - 4\zeta(2)\zeta(3)\tag{3} \end{align*} are used.
Note: $(1)$ is in @nospoon's answer here, $(2)$ can be found here, and for $(3)$ you can see Euler's formula and here.
Evaluation of $I$: By making substitution $x = \sin \theta$ and using the Fourier series of $$ \ln (\sin\theta) = -\ln 2 -\sum_{k=1}^\infty \frac{ \cos(2k \theta)}{k}, $$ we get \begin{align*} I =& \int_{0}^{\frac\pi 2} \theta^4 \ln(\sin\theta)\ \mathrm d\theta\\ =&\int_{0}^{\frac\pi 2} \theta^4\left(-\ln 2 -\sum_{k=1}^\infty \frac{ \cos(2k \theta)}{k}\right)\ \mathrm d\theta\\ =& -\ln 2\int_0^{\frac \pi 2}\theta^4\ \mathrm d\theta-\sum_{k=1}^\infty \frac{1}{k}\underbrace{\int_{0}^{\frac\pi 2}\theta^4 \cos(2k \theta) \ \mathrm d\theta}_{\text{IBP}\times 4}\\ =& -\frac{\pi^5\ln 2}{160}-\sum_{k=1}^\infty \frac{1}{k}\cdot\left(-\frac{\pi^3}{8}\frac{(-1)^{k-1}}{k^2} +\frac{3\pi}{4}\frac{(-1)^{k-1}}{k^4}\right)\\ =&-\frac{\pi^5\ln 2}{160}+\frac{\pi^3}8\underbrace{\sum_{k=1}^\infty \frac{(-1)^{k-1}}{k^3}}_{=\frac 3 4 \zeta(3)} - \frac{3\pi}4\underbrace{\sum_{k=1}^\infty \frac{(-1)^{k-1}}{k^5}}_{=\frac{15}{16}\zeta(5)}\\ =&\boxed{-\frac{\pi^5\ln 2}{160}+\frac{3\pi^3}{32}\zeta(3) -\frac{45\pi}{64}\zeta(5).} \end{align*}
Combining these, we get the equation $$ 3\pi S-\frac{87\pi}{16} \zeta(5) +\frac{\pi^3}{4}\zeta(3) -\frac{\pi^5\ln 2}{160}=-\frac{\pi^5\ln 2}{160} +\frac{3\pi^3}{32}\zeta(3)-\frac{45\pi}{64}\zeta(5), $$hence it follows $$ \boxed{S = \frac{101}{64}\zeta(5) -\frac{5\pi^2}{96}\zeta(3).} $$
Addendum: By considering MacLaurin series of \begin{align*} \ln(1-x)\ln(1+x) =&-\sum_{k=1}^\infty \left(\frac{ H_{2k}}k-\frac{H_k}{k} + \frac1{2k^2}\right)x^{2k} \end{align*} and \begin{align*} \frac{H_k}{k^2} + \frac{H_k^{(2)}}{k} -\frac{\zeta(2)}{k} =& \frac{\partial }{\partial k}\left[-\frac{H_k}{k}\right]\\ =& \int_0^1 x^{k-1}\ln x\ln(1-x)\ \mathrm dx\\ =&4\int_0^1 x^{2k-1}\ln x \ln(1-x^2)\ \mathrm dx \end{align*} we have that \begin{align*} &\int_{0}^{1}\ln(1-x)\ln(1+x) \frac{\ln x\ln(1-x^2)}x \ \mathrm dx \\&=-\sum_{k=1}^\infty \left(\frac{ H_{2k}}k-\frac{H_k}{k} + \frac1{2k^2}\right)\int_{0}^{1}x^{2k-1} \ln x \ln(1-x^2)\ \mathrm dx \\ &=-\frac 1 4\sum_{k=1}^\infty \left(\frac{ H_{2k}}k-\frac{H_k}{k} + \frac1{2k^2}\right)\left(\frac{H_k}{k^2} + \frac{H_k^{(2)}}{k} -\frac{\zeta(2)}{k}\right). \end{align*} The integral can be attacked by considering algebraic identity $$ ab(a+b) = \frac 1 3 (a+b)^3 - \frac {a^3}3 -\frac{b^3}3 $$ with $a=\ln(1-x)$ and $b=\ln(1+x)$, and extant results.
For the sum, after expanding the summand, the only tricky part is $$ \sum_{k=1}^\infty\frac{H_{2k}H_k}{k^3}, $$ which can be found here. Then, the sum $\sum_{k=1}^\infty \frac{H_{2k}H_k^{(2)}}{4k^2}$ can be evaluated by solving the equation obtained.
Perhaps this has something to do with the theory of multiple zeta values. Here is a self-contained solution using basic properties of $\psi_1(n)=\sum_{m\geqslant n}1/m^2$, namely $0<\psi_1(n)-1/n<1/n^2$: $$\small\psi_1(n)-\frac1n=\sum_{m\geqslant n}\left(\frac1{m^2}-\frac1{m(m+1)}\right)=\sum_{m\geqslant n}\frac1{m^2(m+1)}<\sum_{m\geqslant n}\frac1{nm(m+1)}=\frac1{n^2}$$ (which implies $\lim\limits_{n\to\infty}n\psi_1(n)=1$ and the convergence of $S$) and $\sum_{n\geqslant 1}\psi_1(n)/n=2\zeta(3)$, equivalent to the simplest Euler sum a.k.a. $\zeta(2,1)$; an easy proof is as follows: $$\small\Xi=\sum_{n\geqslant 1}\frac1n\sum_{m\geqslant n}\frac1{m^2}=\sum_{m\geqslant 1}\frac1{m^2}\sum_{1\leqslant n\leqslant m}\frac1n=\sum_{m\geqslant 1}\frac1{m^2}\sum_{n\geqslant 1}\left(\frac1n-\frac1{n+m}\right)=\sum_{n,m\geqslant 1}\frac1{nm(n+m)}\\\small=\sum_{n,m\geqslant 1}\left(\frac1n+\frac1m\right)\frac1{(n+m)^2}=2\sum_{n,m\geqslant 1}\frac1{n(n+m)^2}=2\sum_{n\geqslant 1}\frac1n\sum_{m>n}\frac1{m^2}=2\big(\Xi-\zeta(3)\big).$$
We compute (in two ways) the limit $\lim\limits_{N\to\infty}(S_N-H_N)$, where $$\color{gray}{H_N=\sum_{n=1}^N\frac1n,}\qquad S_N=\sum_{n=1}^N n\left(\sum_{m=n}^N\frac1{m^2}\right)^2.$$
First, writing $S_{N-1}=\sum_{n=1}^{N-1}n\big(\psi_1(n)-\psi_1(N)\big)^2$ and expanding the square, we get $$S_{N-1}-H_{N-1}=\sum_{n=1}^{N-1}\left(n\psi_1^2(n)-\frac1n\right)-2\psi_1(N)\sum_{n=1}^{N-1}n\psi_1(n)+\psi_1^2(N)\frac{N(N-1)}{2}.$$
Now $\lim\limits_{n\to\infty}n\psi_1(n)=1$ implies $\lim\limits_{N\to\infty}(1/N)\sum_{n=1}^{N-1}n\psi_1(n)=1$ by Stolz–Cesàro, hence $$\lim_{N\to\infty}(S_N-H_N)=\color{blue}{S-\frac32}.$$
Second, fix $N$ and put $A_n=n(n+1)/2$ and $B_n=\sum_{m=n}^N 1/m^2$ (with $B_{N+1}=0$), then $$S_N=\sum_{n=1}^N(A_n-A_{n-1})B_n^2=\sum_{n=1}^N A_n(B_n^2-B_{n+1}^2)\\\small=\sum_{n=1}^N A_n(B_n-B_{n+1})(B_n+B_{n+1})=\sum_{n=1}^N\frac{n+1}{2n}\left(2B_n-\frac1{n^2}\right)$$ and, since $\sum_{n=1}^N B_n=H_N$, and $B_n\to\psi_1(n)$ monotonically (as $N\to\infty$), we get $$\lim_{N\to\infty}(S_N-H_N)=\sum_{n=1}^\infty\left(\frac{\psi_1(n)}{n}-\frac{n+1}{2n^3}\right)=\color{blue}{\frac32\zeta(3)-\frac12\zeta(2)}.$$
Best Answer
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
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}$$