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)$$?

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 $$$$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$$$$ The orthogonality relation for the Laguerre polynomials reads $$$$\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}$$$$ It can be differentiated with respect to $$x$$ to obtain $$$$\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}$$$$ From the differentiation relation $$$$\frac{d}{dx}L_n^{(x-1)}(t)=\sum_{k=0}^{n-1}\frac{L_k^{(x-1)}(t)}{n-k}$$$$ 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$$: $$$$\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}$$$$ which is the proposed expression.