Proving interesting identities for an integral involving product of confluent hypergeometric functions. While working on a problem I came across the following interesting result.
Let:
$$
H_{nm}(x)=\int_0^\infty t^{x-1}e^{-t}\log t\;F(-n;x;t)F(-m;x;t)\;dt,
$$
where $n,m$ are non-negative integer numbers, $x$ is positive real number and
$$F(a;b;t)=\sum_{k\ge0}\frac{a^{\overline k}}{b^{\overline k}}\frac{t^k}{k!}$$
is the confluent hypergeometric function.
Since the integral is symmetric with respect to permutation of $n$ and $m$ in what follows $n\ge m$ is assumed.
By numerical evidence the integral evaluates to the following simple values:
$$
H_{nm}(x)=\frac{n!\;\Gamma^2(x)}{\Gamma(x+n)}\times
\begin{cases}
\displaystyle\psi(x+n),& n=m; \\
\displaystyle\frac1{m-n},&n\ne m.
\end{cases}\tag1
$$
where $\Gamma(x)$ and $\psi(x)$ are the gamma and digamma functions, respectively.
Is there a simple way to prove the relations $(1)$? 
 A: The confluent hypergeometric functions are related to the generalized Laguerre polynomials:
\begin{align}
F(-n;x;t)&=\frac{\Gamma(n+1)\Gamma(x)}{\Gamma(x+n)}L_n^{(x-1)}(t)
\end{align}
so
\begin{equation}
H_{n,m}(x)=\frac{n!m!\Gamma^2(x)}{\Gamma(x+n)\Gamma(x+m)}\int_0^\infty t^{x-1}e^{-t}\ln t L_n^{(x-1)}(t)L_m^{(x-1)}(t)\,dt
\end{equation} 
The orthogonality relation for the Laguerre polynomials reads
\begin{equation}
\int_0^\infty t^{x-1}e^{-t}L_n^{(x-1)}(t)L_m^{(x-1)}(t)\,dt=\frac{\Gamma(n+x)}{n!}\delta_{n,m}
\end{equation} 
It can be differentiated with respect to $x$ to obtain
\begin{equation}
\int_0^\infty t^{x-1}e^{-t}\ln t L_n^{(x-1)}(t)L_m^{(x-1)}(t)\,dt+\int_0^\infty t^{x-1}e^{-t}\frac{d}{dx}\left[L_n^{(x-1)}(t)L_m^{(x-1)}(t)\right]\,dt=\frac{\Psi(n+x)\Gamma(n+x)}{n!}\delta_{n,m}
\end{equation}
From the differentiation relation
\begin{equation}
\frac{d}{dx}L_n^{(x-1)}(t)=\sum_{k=0}^{n-1}\frac{L_k^{(x-1)}(t)}{n-k}
\end{equation} 
and recognizing the definition of $H_{n,m}(x)$, we have thus
\begin{align}
\frac{\Gamma(n+x)\Gamma(m+x)}{n!m!\Gamma^2(x)}H_{n,m}(x)&+\sum_{k=0}^{n-1}\frac{1}{n-k}\int_0^\infty t^{x-1}e^{-t}L_k^{(x-1)}(t)L_m^{(x-1)}(t)\,dt\\
&+\sum_{k=0}^{m-1}\frac{1}{m-k}\int_0^\infty t^{x-1}e^{-t}L_k^{(x-1)}(t)L_n^{(x-1)}(t)\,dt\\
&=\frac{\Psi(n+x)\Gamma(n+x)}{n!}\delta_{n,m}
\end{align} 
using the orthogonality relation and supposing that $n> m$, only one term in the sums survives, while there is no if $n=m$:
\begin{equation}
\frac{\Gamma(n+x)\Gamma(m+x)}{n!m!\Gamma^2(x)}H_{n,m}+\frac{1}{n-m}\frac{\Gamma(m+x)}{m!} \left( 1-\delta_{n,m} \right)=\frac{\Psi(n+x)\Gamma(n+x)}{n!}\delta_{n,m}
\end{equation} which is the proposed expression.
