Different ways to evaluate $\int_0^\infty xI_0(2x)e^{-x^2}\,dx$ Evaluate $$\int_0^\infty xI_0(2x)e^{-x^2}\,dx$$ where $$I_0(x) = \frac 1\pi \int_0^\pi e^{x\cos\theta}\,d\theta$$ is a Bessel Function. 
Source: Berkeley Math Tournament 
This question was on a math contest for high school students, so I am looking for other methods that preferably do not involve higher mathematics than Calc II. However, I am also interested in other ways to solve this problem that goes beyond the normal calculus curriculum. My solution is posted below as an answer.
 A: A simpler way is to exploit Fubini's theorem and polar coordinates:
$$ \int_{0}^{+\infty}x\,I_0(2x) e^{-x^2}\,dx = \frac{1}{\pi}\int_{0}^{\pi}\int_{0}^{+\infty}\rho e^{-2\rho\cos\theta} e^{-\rho^2}\,d\rho\,d\theta $$
where the last integral equals:
$$ \frac{1}{\pi}\int_{0}^{+\infty}\int_{-\infty}^{+\infty} e^{-2x} e^{-x^2-y^2}\,dx\,dy \stackrel{x\mapsto x-1}{=}\frac{e}{\pi}\int_{0}^{+\infty}\int_{-\infty}^{+\infty}e^{-x^2-y^2}\,dx\,dy$$
that simplifies to:
$$ \frac{e}{2\pi}\iint_{\mathbb{R}^2} e^{-x^2-y^2}\,dx\,dy = \color{red}{\frac{e}{2}}.$$

This also gives the straightforward generalization
$$ \int_{0}^{+\infty} x\,I_0(\rho x)e^{-x^2}\,dx = \color{red}{\frac{1}{2}\,e^{\rho^2/4}}.$$
A: Note that $$I=\int_{0}^{\infty}xI_{0}\left(2x\right)e^{-x^{2}}dx\stackrel{e^{-x^{2}}=u}{=}\frac{1}{2}\int_{0}^{1}I_{0}\left(2\sqrt{-\log\left(u\right)}\right)du$$ $$=\frac{1}{2}\sum_{k\geq0}\frac{\left(-1\right)^{k}}{\left(k!\right)^{2}}\int_{0}^{1}\log^{k}\left(u\right)du=\frac{1}{2}\sum_{k\geq0}\frac{1}{k!}=\color{red}{\frac{e}{2}}.$$
A: $$I = \int_0^\infty xI_0(2x)e^{-x^2}\,dx$$
$$I_0(x)=\frac 1\pi\int_0^\pi e^{x\cos\theta}d\theta \stackrel{\theta\rightarrow\theta/2}{=}\frac 1{2\pi}\int_0^{2\pi} e^{x\cos(\theta/2)}d\theta \stackrel{z=\exp(i\theta/2)}{=}\frac 1{2\pi i}\int_{|z|=1} \frac{xe^{[z+1/z]/2}}{z}dz$$
To evaluate this, we proceed by expanding the function as a series and finding the residue, using the residue theorem. We have that $$I_0(x)=\text{Res}(\frac{{e^{x[z+1/z]/2}}}z,z=0) = \text{constant term of } 1+x\frac{z+z^{-1}}2+x^2\frac{(z+z^{-1})^2}{2^{2} 2!}+x^3\frac{(z+z^{-1})^3}{2^3 3!} +\cdots$$
Noting that we only get a constant term contribution from the even terms, and the constant of $(z+z^{-1})^{2n}$ is $\binom{2n}n,$ we have that $$I_0(x)=\text{Res}(\frac{{e^{z+1/z}}}z,z=0) = \sum_{n\geq0} \frac{x^{2n}}{2^{2n}(2n)!}\binom{2n}{n}=\sum_{n\geq0} \frac{x^{2n}}{2^{2n}(n!)^2}$$
Then $$I_0(2x)=\sum_{n\geq0} \frac{x^{2n}}{(n!)^2}$$
So $$I = \int_0^\infty \sum_{n\geq0} \frac{x^{2n}}{(n!)^2} xe^{-x^2}\,dx\stackrel{\text{terms are all positive}}= \sum_{n\geq0}\frac{1}{(n!)^2} \int_0^\infty{x^{2n+1}}e^{-x^2}\,dx\stackrel{u=x^2}=\frac12 \sum_{n\geq0}\frac{1}{(n!)^2} \int_0^\infty{x^{n}}e^{-u}\,du = \frac 12\sum_{n\geq0}\frac{1}{(n!)^2} \Gamma(n+1) = \frac 12\sum_{n\geq0}\frac{1}{(n!)} =\boxed{\frac e2}$$
A: 
I thought it might be instructive to present a way forward the relies on only the series expansion of the exponential function, evaluating two integrals using reduction formulae, and straightforward arithmetic.  To that end, we proceed.

Using the Taylor series for $e^t=\sum_{n=0}^\infty \frac{t^n}{n!}$, with $t=2x\cos(\theta)$, we can write $I_0(2x)$ as
$$\begin{align}
I_0(2x)&=\frac1\pi \int_0^\pi e^{2x\cos(\theta)}\,d\theta\\\\
&=\frac1\pi \sum_{n=0}^\infty \frac{(2x)^n}{n!}\int_0^\pi \cos^n(\theta)\,d\theta\tag 1
\end{align}$$
Next, using the reduction formula for $\int_0^\pi \cos^n(\theta)\,d\theta=\frac{n-1}{n}\int_0^\pi \cos^{n-2}(\theta)\,d\theta$ reveals
$$\int_0^\pi \cos^n(\theta)\,d\theta=\begin{cases}\pi\frac{n!}{(n!!)^2}&,n\,\text{even}\\\\0&,n\,\text{odd}\tag2\end{cases}$$
Using $(1)$ and $(2)$, we find that 
$$\begin{align}
\int_0^\infty xe^{-x^2}I_0(2x)\,dx&= \frac1\pi\sum_{n=0}^\infty \underbrace{\frac{4^n}{(2n)!}\left(\pi\,\frac{(2n)!}{((2n)!!)^2}\right)}_{=\frac{\pi}{(n!)^2}}\,\,\underbrace{\int_0^\infty x^{2n+1}e^{-x^2}\,dx}_{=\frac12 n!}\\\\
&=\frac12\sum_{n=0}^\infty \frac{1}{n!}\\\\
&=\frac{e}{2}
\end{align}$$
where we used the reduction formula $\int_0^\infty x^{2n+1}e^{-x^2}\,dx=n\int_0^\infty x^{2n-1}e^{-x^2}\,dx$, along with the elementary integral $\int_0^\infty xe^{-x^2}\,dx=\frac12$, to establish $\int_0^\infty x^{2n+1}e^{-x^2}\,dx=\frac12 n!$. 

Tools Used:  The Taylor Series for $e^x$, the reduction formula for $\int_0^\pi \cos^n(x)\,dx$, and the reduction formula for $\int_0^\infty x^{2n+1}e^{-x^2}\,dx$

