Specific integral of modified Bessel function of first kind While doing some calculations for a random matrix model I arrived at the following formula
$$ \int_0^x\int_x^∞ e^{−t−s} I_0(2\sqrt{ts}) \mathrm{d}t\mathrm{d}s 
= xe^{−2x} (I_0(2x) + I_1(2x)) $$
where $I_ν$ denotes the modified Bessel function of first kind, i.e.
$$I_ν(x) = \sum_{n=0}^{\infty} \frac{1}{n!Γ(n+ν+1)}\left(\frac{z}{2}\right)^{2n+ν}.$$
I am now looking for some (short) analytic proof for it? I tried a bit around but could not find any.

A few thoughts that may propel you onto a full proof.
Consider taking the following integral over the rectangle $0 \leq u \leq x, 0 \leq t \leq y$
$$I(x, y) = \int_0^x \int_0^ye^{-u-t}I_0(2 \sqrt{ut})\,\mathrm{du}\, \mathrm{dt}$$
I aim to show that
$$I(x, y) = x-\frac{1}{2}e^{-x-y}[(x+y)I_0(\xi)+\xi I_1(\xi)]+(y-x)F(x, y) \tag{1}$$
Where
$$F(x, y) = \frac{1}{2}\kappa\int_\xi^\infty e^{- \sigma t}f(t)\, \mathrm{dt}$$
and
\begin{align}
f(t) &= e^{-t}I_0(t)\\
\xi  &=2\sqrt{xy} \\
\sigma &= \frac{(\sqrt{y}-\sqrt{x})^2}{\xi}\\
\kappa &= \frac{y-x}{\xi}
\end{align}
So, write
\begin{align}
I(x, y)&=x+(y-x)K(x, y)-e^{-x-y}[\frac{1}{2}\xi I_1(\xi)+xI_0(\xi)] \tag{2}\\
K(x, y) &=\int_0^xe^{-(t+y)}I_0(2 \sqrt{ty})\,\mathrm{dt} \tag{3}
\end{align}
This relation comes straight from Lassey On the computation of certain integrals containing the modified Bessel function $I_0(\xi)$.
The second relation here can be verified by applying $\partial^2 / \partial x \partial y$ to the first relation and using the fact that
\begin{align}
\frac{\partial K}{\partial y} &= -e^{-x-y}\sqrt{x/y}I_1(\xi)\\
\frac{\partial}{\partial x}\sqrt{x/y}I_1(\xi)&=I_0(\xi)
\end{align}
Now, proceeding with $(3)$ and replacing the Bessel function with 
$$I_0(2\sqrt{ty}) = \frac{1}{2 \pi i}\int_L e^{t/s + ys}\,\frac{\mathrm{ds}}{s}\tag{4}$$
Where $L$ is the unit circle around $s=0$. Inserting $(4)$ into $(2)$ yields
$$K(x, y) = \frac{e^{-x-y}}{2 \pi i}\int_L \frac{e^{ys}(e^{x/s}-e^x)}{1-s}\, \mathrm{ds}$$
Integrating this along $|s| = \rho = \sqrt{x/y}$ (assuming that $\rho < 1$ then
$$K(x,y)=- \frac{e^{-x-y}}{2 \pi}\int_0^\pi e^{\xi \cos \theta}\frac{2 \rho (\rho-\cos \theta)}{\rho^2 - 2\rho \cos \theta +1}\, \mathrm{d \theta} $$
Writing $2 \rho (\rho - \cos \theta) = \rho^2 - 2 \rho \cos \theta + 1 + \rho^2-1$ gives
$$K(x, y) = F(x, y)-\frac{1}{2}e^{-x-y}I_0(\xi)\tag{5}$$ 
In the above we have integrated $(4)$ along the circle $|s|=\sqrt{t/y}$ to obtain
$$I_0(\xi) = \frac{1}{\pi}\int_0^\pi e^{\xi \cos \theta}\, \mathrm{d \theta}$$
Returning to $(5)$ we define
$$F(x, y) = \frac{(1-\rho^2)e^{-x-y}}{2 \pi} \int_0^\pi \frac{e^{\xi \cos \theta}}{\rho^2-2\rho\cos \theta +1}\, \mathrm{d \theta}$$
Combining $(5)$ and $(2)$ you should get to $(1)$. I hope that this gives some help, please comment and critique as you see fit!
