# 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

$$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{{\displaystyle #1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\on}[1]{\operatorname{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\sr}[2]{\,\,\,\stackrel{{#1}}{{#2}}\,\,\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$$ \begin{align} \on{I}\pars{a,b,c,j,k} & \equiv\color{#44f} {\int_{0}^{\infty}\int_{0}^{\infty} {\expo{-\pars{x + y}/2}\,\,\,x^{j}\,y^{k} \over ax + by + c} \,\dd x\,\dd y} \\[5mm] & \sr{\substack{\,\,\,\,x/c\ \mapsto\ x\\[0.5mm] \atop y/c\ \mapsto\ y\\[0.5mm] \mbox{}}}{=} c^{j + k + 1}\int_{0}^{\infty}\int_{0}^{\infty} {\expo{-c\pars{x + y}/2}\,\,\,x^{j}\,y^{k} \over ax + by + 1} \,\dd x\,\dd y \\[5mm] & \sr{{\rm as}\ c\ \to\ \infty}{\sim} c^{j + k + 1}\pars{\int_{0}^{\infty}\expo{-cx/2}\,x^{j}\,\dd x} \pars{\int_{0}^{\infty}\expo{-cy/2}\,y^{k}\,\dd y}\tag{1}\label{1} \\[5mm] & = c^{j + k + 1}\ \bracks{\pars{2 \over c}^{j + 1}\Gamma\pars{j + 1}} \bracks{\pars{2 \over c}^{k + 1}\Gamma\pars{k + 1}} \end{align} (\ref{1}) involves the $$\underline{Laplace's\ Method}$$ because the main contribution to the integral, as $$\ds{c \to \infty}$$, occurs in a neighborhood of "$$\ds{x = 0\ \mbox{and}\ y = 0}$$". Therefore, \begin{align} \on{I}\pars{a,b,c,j,k} & \equiv\color{#44f} {\int_{0}^{\infty}\int_{0}^{\infty} {\expo{-\pars{x + y}/2}\,\,\,x^{j}\,y^{k} \over ax + by + c} \,\dd x\,\dd y} \\[5mm] & \sr{{\rm as}\ c\ \to\ \infty}{\sim} \bbx{\color{#44f}{2^{j + k + 2}\,\,\,\Gamma\pars{j + 1}\Gamma\pars{k + 1}} \over c} \end{align}
• @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$$.