Asymptotic expansion of multiple integral at $\infty$ I am attempting to figure out the first term of the asymptotic expansion at $t \to \infty$ of the following triple integral:
$$
I(t) = \int_0^1 \int_0^\pi \int_0^\pi e^{[r\cdot(\cos(x)-\cos(y))+2 \cos(y)]\cdot t} \sin^2(x) \sin^2(y) \mathrm d x \mathrm dy \mathrm dr
$$
Intuitively, I am expecting to find something like $I(t)\sim \gamma e^{2t} t^\alpha$ with some constants $\alpha, \gamma$, as the exponential is most certainly going to dominate for such values (for instance when $y=0, r=0$)
As for solving it, I am considering using Laplace's method - or at least a variant of it. Now, it's a bit puzzling because the inward term in the exponential function : $f(r,x,y) = r\cdot(\cos(x)-\cos(y))+2 \cos(y)$ has some issues: its maximum value 2 is reached for all values for which we have $y=0$ and either $x=0$ or $r=0$ - so that means potentially two segments on which the exponential is maximal (while the link provides inputs for a single point). Secondly, these points are on the border of the integration surface.
Would you know another variant of this method to solve this properly?
Thank you in advance!
 A: As remarked by @user2673062 result in the comments, the triple integral simplifies. By using the integral representation for the modified Bessel function
$$I_{\nu}\left(z\right)=\frac{(\frac{1}{2}z)^{\nu}}{\pi^{\frac{1}{2}}\Gamma\left(\nu+\frac{1}{2}\right)}\int_{0}^{\pi}e^{\pm z\cos\theta}(\sin\theta)^{2\nu}\mathrm{d}\theta$$
for $\nu>-1/2$, we can perform the integrations over $x$ and $y$ with $\nu=1$ to obtain:
\begin{align} I(t) &= \int_0^1 \int_0^\pi \int_0^\pi e^{[r\cdot(\cos(x)-\cos(y))+2 \cos(y)]\cdot t} \sin^2(x) \sin^2(y) \mathrm d x \mathrm dy \mathrm dr\\
&=\frac{\pi^2}{t^2}\int_0^1\frac{I_1(tr)}r\frac{I_1\left(t(2-r)\right)}{2-r}\,dr
\end{align}
By changing $r\to 2-r$, this integral can be written as
\begin{align}I(t)&= \frac{\pi^2}{2t^2}\int_0^2\frac{I_1(tr)}r\frac{I_1\left(t(2-r)\right)}{2-r}\,dr\\
&=\frac{\pi^2}{2t^2}\Phi(2)
\end{align}
where
$$\Phi(\theta)=\int_0^\theta\frac{I_1(tr)}r\frac{I_1\left(t(\theta-r)\right)}{\theta-r}\,dr$$
can be considered as a convolution integral. Its Laplace transform is then
\begin{align}\mathcal L[\Phi] &=\left\lbrace \mathcal L\left[\frac{I_1(tx)}x\right]\right\rbrace^2\\
&=\frac1{t^2}(u-\sqrt{u^2-t^2})^2
\end{align}
By inversion, we find
$$ \Phi(\theta)=\frac2\theta I_0(\theta t)-\frac4{\theta^2 t}I_1(\theta t)$$
These integral transforms were found using a CAS but could certainly be found from classical tables of Laplace transforms (Ederlyi TI, for example).
Finally,
$$ I(t)=\frac{\pi^2}{2t^2}\left[I_0(2 t)-\frac1{ t}I_1(2t)\right]$$
which seems to be numerically correct.
Using the asymptotic expansions for the Bessel functions, we find
$$ I(t)\sim \frac{\pi^{3/2}}{4}\frac{e^{2t}}{t^{5/2}}$$ which is not very far from @Raffaele estimation.
A: Using the representation of the integral in terms of modified Bessel functions we have:
$$
I(t)=\frac{\pi^2}{t^2}\int_0^1dr \frac{I_1(tr)I_1(t(2-r))}{r(2-r)}
$$
Observation: $t(2-r)\gg1$ for $r\in(0,1)$ so we won't do much harm if we just replace $$
I_1(t(2-r))\approx\frac{e^{2t-r}}{\sqrt{2t\pi(2-r)}}$$
(I am ignoring $O$'s for the moment but all statements can be made precise with a bit of extra work),
We get
$$
I(t)=\underbrace{\frac{\pi^{3/2} e^{2t}}{\sqrt{2}t^{5/2}}}_{C(t)}\underbrace{\int_0^1dr \frac{I_1(tr)e^{-tr}}{r(2-r)^{3/2}}}_{\mathcal{J(t)}}
$$
We now split this integral at $1/t \ll\delta \ll 1$:
$$
\mathcal{J(t)}=\mathcal{J_1(t)}+\mathcal{J_2(t)}
$$
We get (Lemma 1 plus large argument asymptotics):
$$
\mathcal{J_1(t)}\approx \frac{1}{2^{3/2}}\int_0^{\delta} I_1(tr)e^{-tr}/r=\frac{1}{2^{3/2}}((1+O(\sqrt{\delta t}^{-1}))
$$
for $\mathcal{J_2(t)}$ we can use again large argument asymptotics for $I_1(rt)$:
$$\mathcal{J_2(t)}\approx\frac{1}{\sqrt{2\pi t}}\int_{\delta}^1\frac{dr}{r^{3/2}(2-r)^{3/2}}=\frac{\eta}{\sqrt{2\pi t}}-O(\delta^{5/2})
$$
where $\eta$ is a constant independet of $t$.
Finally since $\mathcal{J_2(t)}\ll\mathcal{J_1(t)}$

$$
I(t)\approx C(t)\mathcal{J_1(t)} \approx \frac{C(t)}{2^{3/2}}\\
=\frac{e^{2t}\pi^{3/2}}{4t^{5/2}}
$$


Lemma 1:
$$ \int \frac{dz}ze^{-z}I_1(z)=Const-e^{-z}(I_1(z)+I_0(z))
$$
Proof:
Differentiate w.r.t. $z$ and use $I_0'(z)=I_1(z), \,\, 2 I_1'(z)=I_0(z)+I_2(z)$ as well as $I_0(z)-I_2(z)=I_1(z)/(2z)$
