A double series $\frac13 \sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\frac{(i-1)! (j-1)!}{(i+j)!}H_{i+j}$ giving $\zeta(3)$

Here is a symmetric rational double series giving Apery's constant:

$$ \frac13 \sum_{j=1}^{\infty}\sum_{i=1}^{\infty} \displaystyle \frac{(i-1)! (j-1)!}{(i+j)!} H_{i+j} = \zeta(3) $$

where $\displaystyle H_{n}:=\sum_{1}^{n} \frac{1}{k}$, $n=1,2,\cdots,$ are the harmonic numbers.

How would you prove it?

Edit. In 2005, I sent this result to Wolfram MathWorld (see equation 25), I built it as an echo of

$$ \sum_{j=1}^{\infty}\sum_{i=1}^{\infty} \displaystyle \frac{(i-1)! (j-1)!}{(i+j)!} = \zeta(2). $$

See Jack's pretty answer and see my answer below.

Summation by parts seems to work fine. For $j=1$ we have: $$\sum_{i=1}^{+\infty}\frac{H_{i+1}}{i(i+1)}=2,$$ since $\sum_{i=1}^{N}\frac{1}{i(i+1)}=1-\frac{1}{1+N}$ and so: $$\begin{eqnarray*}\sum_{i=1}^{N}\frac{H_{i+1}}{i(i+1)}&=&H_{N+1}\left(1-\frac{1}{N+1}\right)-\sum_{i=1}^{N-1}\left(1-\frac{1}{i+1}\right)\frac{1}{i+2}\\&=&H_2+\sum_{i=1}^{N-1}\frac{1}{(i+1)(i+2)}+O\left(\frac{\log N}{N}\right)\\&=&H_2+\frac{1}{2}+O\left(\frac{\log N}{N}\right),\end{eqnarray*}$$ since: $$\sum_{i=1}^{+\infty}\frac{1}{(i+1)\cdot\ldots\cdot(i+k)}=\frac{1}{(k-1)k!}.$$ For $j=2$ we have: $$\begin{eqnarray*}1!\sum_{i=1}^{N}\frac{H_{i+2}}{i(i+1)(i+2)}&=&H_{N+1}\left(\frac{1}{4}-\frac{1}{2(N+1)(N+2)}\right)-\sum_{i=1}^{N-1}\left(\frac{1}{4}-\frac{1}{2(i+1)(i+2)}\right)\frac{1}{i+3}\\&=&\frac{H_3}{4}+\frac{1}{4\cdot 3!}+O\left(\frac{\log N}{N}\right)\end{eqnarray*}$$ and in the general case we have: $$(j-1)!\sum_{i=1}^{+\infty}\frac{H_{i+j}}{i\cdot\ldots\cdot(i+j)}=(j-1)!\left(\frac{H_{j+1}}{j^2(j-1)!}+\frac{1}{j^2(j+1)!}\right)=\frac{H_{j+1}}{j^2}+\frac{1}{j^3(j+1)},$$ so the original sum equals: $$\sum_{j=1}^{+\infty}\frac{H_{j+1}}{j^2}+\sum_{j=1}^{+\infty}\frac{1}{j^3(j+1)}=\sum_{j=1}^{+\infty}\frac{H_j}{j^2}+\sum_{j=1}^{+\infty}\frac{1}{j^3}=2\zeta(3)+\zeta(3)=3\zeta(3).$$

This is another path.

Theorem. Let $x$ be a real number such that $-1<x<1$ and set $$ Z(x):=\sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\frac{(i-1)! (j-1)!}{(i+j)!} x^{i+j}. $$ Then $$ Z(x)= 2 \:\mathrm{Li_2} (x) - \mathrm{Li_2} \left( x(2-x) \right) \tag1 $$

where $\mathrm{Li_2}$ is the dilogarithm function such that $\mathrm{Li_2}(x)=\displaystyle \sum_{n=1}^{\infty}\displaystyle \frac{x^{n}}{n^2}= - \displaystyle \int_{0}^{x}\,\displaystyle \frac{\log(1-t)}{t}\mathrm{d}t$.

Proof. Observe that $$ {\frac{(i-1)! (j-1)!}{(i+j)!}= \frac{\Gamma(i) \Gamma(j)}{(i+j) \Gamma(i+j)} = \frac{1}{i+j}\,B(i,j)= \int_{0}^{1}\, \frac{t^{i-1}(1-t)^{j-1}}{i+j}\mathrm{d}t.} $$ Let $|x|<1$. We can obtain $$ \begin{align*} \sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\frac{(i-1)! (j-1)!}{(i+j)!} x^{i+j} &= \sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\displaystyle \int_{0}^{1} \frac{t^{i-1}(1-t)^{j-1}}{i+j}dt \:x^{i+j} \\ &= x^2 \int_{0}^{1} \sum_{j=1}^{\infty}\,(x(1-t))^{j-1}\sum_{i=1}^{\infty} \frac{(tx)^{i-1}}{i+j}\:\mathrm{d}t \\ &= x^2 \int_{0}^{1} \sum_{j=1}^{\infty}\,(x(1-t))^{j-1}(tx)^{-1-j}\sum_{i=j+1}^{\infty} \frac{(tx)^{i}}{i}\:\mathrm{d}t\\ &= x^2 \int_{0}^{1} \sum_{j=1}^{\infty}\,(x(1-t))^{j-1}(tx)^{-1-j}\int_{0}^{tx} \frac{u^{j}}{1-u}\:\mathrm{d}u\:\mathrm{d}t \\ &= -x \int_{0}^{1} \frac{\log(1-x + x u)}{1-x u} \:\mathrm{d}u \\ &= -x \int_{1/x - 1}^{1/x} \frac{ \log(1-\frac{x}{2-x }u) + \log(\frac{2-x}{x})}{u} \:\mathrm{d}u\\ &= 2 \: \mathrm{Li_2} (x) - \mathrm{Li_2} \left( x(2-x) \right) \end{align*}$$

Example 1. $$ \sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\frac{(i-1)! (j-1)!}{(i+j)!}= \frac{\pi^2}{6} $$ Put $x=1$ in $(1)$.

Example 2. $$ \sum_{j=1}^{\infty}\sum_{i=1}^{\infty}\frac{(i-1)! (j-1)!}{(i+j)!} \frac{1}{ \varphi^{2(i+j)}}= \frac{\pi^2}{30}-\ln^2 \varphi $$ Put $x=1/\varphi^2$ in $(1)$, with $\varphi=\frac{1+\sqrt{5}}{2}$.

Proposition. $$ \sum_{j=1}^{\infty}\sum_{i=1}^{\infty} \displaystyle \frac{(i-1)! (j-1)!}{(i+j)!} H_{i+j} = 3 \: \zeta(3) \tag2 $$

Proof. From the identity $$ 1+x+\cdots+x^{i+j-1} = \frac{1- x^{i+j}}{1 - x\,}, \qquad i\geq1,\,j\geq1, 0<x<1, $$ we write $$ H_{i+j} = \displaystyle \int_{0}^{1} \displaystyle \frac{1- x^{i+j}}{1 - x\,}\: \mathrm{d}x. $$ Then, using $(1)$, we obtain $$ \begin{align} \sum_{j=1}^{\infty}\sum_{i=1}^{\infty} \frac{(i-1)! (j-1)!}{(i+j)!} \: H_{i+j} & =\int_{0}^{1} \frac{ Z(1) - Z(x)}{1 - x}\: \mathrm{d}x \\ & = \int_{0}^{1} \ln(1-x)\left(Z(1) - 2 \: \mathrm{Li_2} (x) + \mathrm{Li_2} \left( x(2-x) \right) \right)'\: \mathrm{d}x \\& = 2 \int_{0}^{1} \dfrac{\ln^2 (1-x)}{2-x}\: \mathrm{d}x \\ & = 2 \int_{0}^{1} \dfrac{\ln^2 x}{1+x}\: \mathrm{d}x \nonumber \\ & = 2 \int_{0}^{1} \sum_{n=0}^{\infty}(-1)^{n} x^{n} \ln^2 x\: \mathrm{d}x \\ & = 3 \: \zeta(3). \end{align} $$