Expansion of double integral I've come across this deceptively difficult integral:
$$I(t)=\frac{e^{4t}}{t^2}\int_1^\infty dx\int_1^\infty dy\frac{e^{-2t(x+y)}}{(4xy-1)^2},$$
and I want to study its behavior around $t\rightarrow0^+$. I found a simple series expression in terms of exponential integrals, but it doesn't seem to help much:
$$I(t)=\frac{e^{4t}}{t^2}\sum_{n=0}^\infty\frac{n+1}{4^{n+2}}E_{n+2}^2(2t).$$
What are the typical tricks and techniques to handle this sort of problem? I expect that there will be terms like $1/t^2$ and $1/t$ with some factors of $\log(2t)$ and constant bits, but precisely extracting them has proven to be a challenge. Thanks in advance!
 A: The following is not rigorous, and maybe someone can improve it.  The approximation will come from an alternative representation,
$$H(a,t):=\int_1^\infty dx \, \int_1^\infty dy \, \frac{ \exp{\big(-2t(x+y)\big)} }{ (a+4xy)^2 } $$
$$ = \frac{1}{4}\int_0^\infty  \exp{\big(-au+t^2/u\big)} \Gamma(0,\frac{(t+2u)^2}{u} ) \, du $$  where the incomplete gamma function can be connected to the exponential function Ei function through
$$\Gamma(0,x) = -\text{Ei(-x)}. $$
The OP's problem is solved by setting $a=-1.$
The proof is not difficult. Use the Gauss integral for $1/(a+4xy)^2,$
$$\int_0^\infty u \, e^{-bu} du = 1/b^2, $$
interchange integrals, let Mathematica do the innermost double integral, and simplify.  The equivalence has been checked numerically for several real $a$ and $t.$  Here is where the non-rigor enters.  There are several representations for Ei(-x) that have a leading $e^{-x}$ factor, e.g., Gradshteyn and Ryzhik 8.212.[3-14] and even the asymptotic form has such a leading factor.  If we used such a formula, then
$$  \exp{\big(-au+t^2/u\big)} \Gamma(0,\frac{(t+2u)^2}{u} ) \approx 
\exp{\big(-au+t^2/u  -(t^2/u + 4t+4u)   \big)} (something) $$
where something should be well-behaved.  The $t^2/u$ terms, which diverge at the origin, are thus cancelled.  (This also justifies the above formal manipulations as long as $a>-4.$) We are therefore left with the much easier to analyze
$$ H(a,t) \approx \frac{1}{4}\int_0^\infty  \exp{\big(-au)} \Gamma(0, 4t+4u) \, du :=\hat{H}(a,t) $$
Shift the argument in $\hat{H},$
$$4\hat{H}(a,t)=\exp{(at)}\int_t^\infty   \exp{\big(-au)} \Gamma(0, 4u) \,du =$$
$$=\exp{(at)}\Big( \int_0^\infty   \exp{\big(-au)} \Gamma(0, 4u)\,du
- \int_0^t\exp{\big(-au)} \Gamma(0, 4u) \,du \Big) $$
The first integral is solvable in closed form.  Since $t \to 0,$ $u$ is small and so approximate the exponential with the constant 1. That integral is solvable as well.  We find that
$$4\hat{H}(a,t) \sim \exp{(at)}\Big( \frac{\log{(1+a/4)}}{a} -(\frac{1}{4}(1-e^{-4t}) + t\Gamma(0, 4t)) \Big) $$
$$\sim \exp{(at)}\Big( \frac{\log{(1+a/4)}}{a} -t (1+ \Gamma(0, 4t)) \Big) $$
One can expand $\Gamma(0, 4t) \sim -(\gamma + \log{4} +\log{t})$ to see that the next to leading order has a $t \log{t}$ dependence.
With $a=-1,$ I did some numerical comparisons with the 2nd equation of the answer vs. the approximation of the preceding formula for $\hat{H}(-1,t)$ for $t=0.01, 0.001,$ and $0.0001.$  The respective errors are 2.0%, 0.23%, and 0.024%, respectively.  I don't think further terms are available via this method until the something mentioned in the middle of the answer is quantified.
