Exponential integrals

The formal solution of the radiative transfer equation can be expressed using exponential integrals. These are mathematical functions defined as:
\begin{equation}
E_n(x) = \int_1^{\infty}\, \exp (-yx)\frac{\mathrm{d}x}{y^n}, \;\;\;\;n = 1, 2, ...
\end{equation}
There is several important properties of these functions:


(i) Derivative of the first order exponential integral is:
\begin{equation}
E_1^{\prime}(x) = -\exp(-x)/x.
\end{equation}

(ii) At $x= 0$, for any $n \ge 2$:
\begin{equation}
E_n(0) = \frac{1}{n - 1},
\end{equation}
while $E_1(0) = +\infty$.

(iii) Theintegrals the following recursion equations:
\begin{eqnarray}
n\,E_{n+1}(x) &=& \exp(-x) - x\,E_n(x),\\
E_{n+1}^{\prime}(x) &=& - E_n(x).
\end{eqnarray}

(iv) For large values of $x$, the asymptotic expansion of $E_n(x)$ is:
\begin{equation}
E_n(x) = \frac{\exp(-x)}{x}\left[1 - \frac{n}{x} + \frac{n(n+1)}{x^2} + \cdots\right].
\end{equation}
For very large values of $x$ only the first term survives and it thus converge to:
\begin{equation}
E_n(x \gg 1) \approx \frac{\exp(-x)}{x}.
\end{equation}

The first four exponential integrals are shown in the figure above (n =1, green;  n = 2, red; n = 3, blue; n = 4, orange). Note their asymptotic behaviour at large x. Note the values they get on the y axis.

Evaluation of exponential integrals

Once $E_1(x)$ is known it is trivial to evaluate the higher-order integrals. However, computing $E_1(x)$ is a bit more complicated. A detailed algorithm is given in chapter 6.3 of Numerical Recipes in Fortran 77: The Art of Scientific Computing (Press et al, 1992, p.215). It is coded in the IDL intrinsic function expint (check EXPINT). Plots in the figure above are rendered using that function.

However, the first exponential integral can be specified by a polynomial fit: The approximative expression for $E_1(x)$ is:
\begin{eqnarray*}
E_1(x) &=& \frac{\exp(-x ) \;(a_0 + a_1\,x + a_2\,x^2 + a_3\,x^3 + a_4\,x^4)}{x\; (b_0 + b_1\,x + b_2\,x^2 + b_3\,x^3 + b_4\,x^4)}\;\;\;\;(x > 1),\\
E_1(x) &=& -\ln (x) + c_0 + c_1\,x + c_2\,x^2 + c_3\,x^3 + c_4\,x^4 + c_5\,x^5\;\;\;\;(x \le 1).
\end{eqnarray*}
The coefficients of the best fit are listed in Abramowitz and Stegun (1964):

	a	b	c
0	0.2677737343	3.9584969228	-0.5772156600
1	8.6347608925	21.0996530827	0.9999919300
2	18.0590169730	25.6329561486	-0.2499105500
3	8.5733287401	9.5733223454	0.0551996800
4	1.0000000000	1.0000000000	-0.0097600400
5			0.0010785700

The approximative solution can be easily coded in IDL, see expint_approx.pro. The difference between $E_1$ computed by the IDL intrinsic figure and by this approximative function is shown in the figure below (the colors are same as before).

References

Many books on radiative transfer has sections on exponential integrals. For example, the detailed derivation of the properties listed above can be found in:
Chandrasekhar, 1950, Radiative Transfer, Appendix I
Kourganoff, 1952, Basic Methods in Transfer Problems, Appendix I