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!
It is known that
$$
I_\nu (x) \sim \frac{{e^x }}{{\sqrt {2\pi x} }}\left( {1 - \frac{{4\nu ^2 - 1}}{{8x}} + \frac{{(4\nu ^2 - 1)(4\nu ^2 - 9)}}{{128x^2 }} - \cdots } \right)
$$
as $x\to +\infty$ with fixed $\nu$ (cf. http://dlmf.nist.gov/10.40.E1). Thus
$$
I_0 (x) \sim \frac{{e^x }}{{\sqrt {2\pi x} }}\left( {1 + \frac{1}{{8x}} + \frac{9}{{128x^2 }} + \cdots } \right),\;\, I_1 (x) \sim \frac{{e^x }}{{\sqrt {2\pi x} }}\left( {1 - \frac{3}{{8x}} - \frac{{15}}{{128x^2 }} - \cdots } \right),
$$
and therefore,
$$
e^{ - x} (I_0 (x) - I_1 (x)) \sim \frac{1}{{\sqrt {2\pi x} }}\left( {\frac{1}{{2x}} + \frac{3}{{16x^2 }} + \cdots } \right) = \frac{1}{{2\sqrt {2\pi } x^{3/2} }}\left( {1 + \frac{3}{{8x}} + \cdots } \right)
$$
as $x\to +\infty$.
Best Answer
From $$ I_n(x) ~=~\frac{1}{\pi}\int_0^{\pi} \! \mathrm{d}\theta ~\exp\left(x\cos\theta\right)\cos n\theta, \qquad n\in~\mathbb{N}_0,\tag{A} $$ we calculate $$\begin{align} \sqrt{x}\pi e^{-x}I_n(x) &~~=~\sqrt{x} \int_0^{\pi} \! \mathrm{d}\theta ~\exp\left(- x(1-\cos\theta)\right)\cos n\theta \cr &\stackrel{t=\sqrt{x}\theta}{=}~ \int_0^{\pi\sqrt{x}} \! \mathrm{d}t ~\exp\left(- x(1-\cos\frac{t}{\sqrt{x}})\right)\cos \frac{nt}{\sqrt{x}} \cr &~~=~ \int_0^{\infty} \! \mathrm{d}t ~\exp\left(- \frac{t^2}{2} + \frac{t^4}{24 x} + O(x^{-2})\right) \left(1- \frac{(nt)^2}{2x} + O(x^{-2})\right) \cr &~~=~ \int_0^{\infty} \! \mathrm{d}t ~\exp\left(- \frac{t^2}{2}\right) \left(1- \frac{(nt)^2}{2x}+ \frac{t^4}{24 x} + O(x^{-2})\right) \cr &\stackrel{u=t^2/2}{=}~ \int_0^{\infty} \! \frac{\mathrm{d}u}{\sqrt{2u}} ~\exp\left(- u\right) \left(1- \frac{n^2u}{x}+ \frac{u^2}{6x} + O(x^{-2})\right) \cr &~~=~\frac{1}{\sqrt{2}}\left( \Gamma(\frac{1}{2}) - \frac{n^2}{x}\Gamma(\frac{3}{2}) +\frac{1}{6x}\Gamma(\frac{5}{2}) + O(x^{-2})\right)\cr &~~=~\sqrt{\frac{\pi}{2}}\left( 1+ \frac{1-4n^2}{8x} + O(x^{-2})\right), \end{align}\tag{B}$$ which agrees with OP's sought-for formulas (1).