I'm trying to obtain an asymptotic large-$k$ approximation for the integral $$I(k) := e^{-k^2}\int_0^1 \frac{(1 + \xi^2)\Gamma(0, \xi^2 k^2) - 2\Gamma(0, k^2)}{1 - \xi} d\xi$$ where $\Gamma$ is the incomplete gamma function $$\Gamma(0, x) = \int_x^\infty \frac{e^{-t}}{t}\,dt$$ (also expressible in terms of the exponential integral). Specifically, I would like to show that as $k\to\infty$, $I(k)$ becomes negligible with respect to $e^{-k^2}\ln k^2$.
My first thought was to use the series expansion $$\Gamma(0, z) = e^{-z}\biggl(\frac{1}{z} - \frac{1}{z^2} + \cdots\biggr)$$ but the problem is that $\xi^2 k^2$ goes all the way down to zero, so it takes values where this series will converge arbitrarily slowly. It might be possible to show that the region of the integral where this effect is significant is small enough not to matter, but I can't think of a way to do that.
I also thought to do a series expansion of the integrand around $\xi = 1$, or $\xi = 0$, and integrate that, or perhaps to split the integral in half and use a different series for each. But the series coefficients I obtained (using Mathematica) depend on positive powers of $k$, increasing with each order of $\xi$ (or $1 - \xi$), which does not bode well for convergence.
Inspired by another question, I tried transforming the expression into a form that would allow me to apply Watson's lemma, but I couldn't find a transformation that would satisfy the conditions.
So, can anyone offer a way to derive the asymptotic behavior? It seems like it shouldn't be that hard, but I'm too tired to keep thinking about this at the moment.
For what it's worth, numerical evaluation suggests that $$I(k) \sim \frac{C e^{-k^2}}{k}$$ where $C = 1.77245\ldots \overset{?}{=} \sqrt{\pi}$. But I would much prefer to have a symbolic demonstration of this, and not just rely on numerics.
Assume we want to show that: $$ \lim_{k\to +\infty}k\int_{0}^{1}\frac{(1+\xi^2)\Gamma(0,\xi^2 k^2)-2\Gamma(0,k^2)}{1-\xi}\,d\xi = C.\tag{3}$$ That is equivalent to: $$ \lim_{k\to +\infty}\sqrt{k}\int_{0}^{1}\frac{(1+\xi^2)\Gamma(0,\xi^2 k)-2\Gamma(0,k)}{1-\xi}\,d\xi = C,\tag{2}$$ and the previous line is a consequence of: $$ \lim_{k\to +\infty} \sqrt{k}\int_{0}^{1}\frac{(1+\xi^2)e^{-k\xi^2}-2e^{-k}}{(1-\xi)}\,d\xi = C \tag{1}$$ by De l'Hopital theorem and differentiation under the integral sign.
Let we consider $f_k(\xi)=\frac{(1+\xi^2)e^{-k\xi^2}-2e^{-k}}{(1-\xi)}$, for large $k$, over $(0,1)$. It is a smooth and positive function that vanishes towards the right endpoint, that is very well approximated by $$ \widetilde{f}_k(\xi) = (1+\xi) \exp\left((2-k)\xi^2\right)\tag{0}$$ by the stationary phase method or by De l'Hopital theorem again.
So we have that $(1)$ holds with $C=\sqrt{\pi}$, and the same happens for $(2)$ and $(3)$.