# Double integral $\int_0^\infty \int_0^\infty \frac{e^{-(x+y)/2}x^j y^k}{ax+by+c} dx dy$

In an ideal scenario I would like to have an analytic expression for the integral $$\operatorname{I}\left(\,{a,b,c,j,k}\,\right) = \int_{0}^{\infty} \int_{0}^{\infty} \dfrac{\operatorname{e}^{-\left(\,{x + y}\,\right)/2}\,\,x^{j}\, y^{k}}{ax + by + c}\,{\rm d}x\,{\rm d}y$$ for positive real $$a,b,c$$ and non-negative integers $$j,k$$, which can then be evaluated for very large values of $$c$$. The reason I mention the latter circumstance is because while Mathematica ($$13.1$$ for sure) is able to take this integral for small values of $$j$$ and $$k$$, the resulting expression is numerically unstable for large values of $$c$$ (meaning that it's a product of terms exponentially large and small in $$c$$). For example, in the case $$j = k = 0$$ one gets: $$\frac{2 e^{\frac{c}{2 b}} \operatorname{Ei}\left(-\frac{c}{2 b}\right)-2 e^{\frac{c}{2 a}} \operatorname{Ei}\left(-\frac{c}{2 a}\right)}{a-b}$$

Would appreciate any advice on how to:

1. Take the original integral analytically,
2. Efficiently evaluate $$\operatorname{e}^\alpha \operatorname{Ei}(-\alpha)$$ for large values of $$\alpha$$.

Given that Mathematica can take this integral for small $$j$$ and $$k$$, I would guess that there exists a recursive relation.

• The title says "indefinite integral", but the integral in the title is a definite integral. Dec 22, 2022 at 21:26
• @Zacky, eq (A.46) in arXiv:1404.6234. Dec 23, 2022 at 2:34
• @Zacky, if you actually take the integral without expanding the Laguerre's polynomials into monomials, my mind will be blown. Dec 23, 2022 at 2:44
• That stuff looks really complicated, but also interesting as I just started reading about quantum computing recently. I like when I see intergrals appearing from different domains. Dec 23, 2022 at 2:56
• Just notice that $E_i(-z)=-E_1(z)=-\int_z^\infty\frac{e^{-t}}{t}dt\,$ and $\,=-\frac{e^{-z}}{z}-\int_z^\infty\frac{e^{-t}}{t^2}dt\,$ (integrating by part). If you continue the process, you get the asymptotics $$\frac{2 e^{\frac{c}{2 b}} \operatorname{Ei}\left(-\frac{c}{2 b}\right)-2 e^{\frac{c}{2 a}} \operatorname{Ei}\left(-\frac{c}{2 a}\right)}{a-b}\sim\frac{2}{a-b}\sum_{k=1}^\infty k!\Big(\frac{2}{c}\Big)^k\big(a^k-b^k\big)$$ The main term (at $k=1$) is equal to $\frac{4}{c}$ and does not depends on $a,b$ - what is visible at the original integral. Dec 23, 2022 at 5:51

• @TymaGaidash Yes. This is the leading contribution as $\displaystyle c \to \infty$. Next terms are found by expanding $\displaystyle\left(ax + by + 1\right)^{-1}$ in powers of $\displaystyle\left(ax + by\right)$. Dec 26, 2022 at 20:05
For the second question $$e^{\alpha } \,\text{Ei}(-\alpha )=\sum _{n=1}^{\infty } (-1)^n \,\frac{(n-1)!}{\alpha^n}$$ You will not need many terms for a very good approximation.
You could also use the $$[k,k+1]$$ Padé approximant $$P_k$$. For example $$P_2=-\frac{\alpha (3 \alpha +13)}{3 \alpha ^3+16 \alpha ^2+10 \alpha -4}$$ whose error is $$\frac {44}{3 \alpha^6 }$$.
For $$\alpha=100$$, this would give $$e^{100 } \,\text{Ei}(-100)=-\frac{7825}{790249}=-0.009901942299$$ while the exact value is $$-0.009901942287$$.
Using $$P_2$$ the relative error is smaller than $$0.01$$% as soon as $$\alpha \geq 8$$ and smaller than $$0.001$$% as soon as $$\alpha \geq 14$$.