Analytical approximation of the definite Chapman integral for arbitrary zenith angles

This study presents an analytical approximation of the definite Chapman integral, applicable to any zenith angle and finite integration limits. The author also presents the asymptotic expression for the definite Chapman integral, which enables an accurate and efficient implementation free of numerical overflows. The maximum relative error in our analytical solution is below 0.5 %.

1 Introduction

The Chapman function, a specific improper integral, has wide application in diverse fields of study (Chapman, 1931, 1953), including but not limited to radiative transfer theory, aeronomy, and atmospheric absorption and scattering in general (Bauer and Lammer, 2004; Brasseur and Solomon, 2005; Schunk and Nagy, 2009; Grieder, 2010; Engel, 2018). It represents the integration of an exponentially varying density along a slanted path within spherical geometry. In computing atmospheric attenuation and scattering over finite distances, the definite form of this integral becomes essential. Several researchers (Green and Barnum, 1963; Fitzmaurice, 1964; Swider and Gardner, 1969; Titheridge, 1988; Kocifaj, 1996; Huestis, 2001) have proposed various analytical approximations of the Chapman function. A comprehensive review and an improvement of these approximations were recently offered by Vasylyev (2021). Nonetheless, a straightforward solution applicable to arbitrary path angles and finite integration limits remains elusive. Our work addresses this gap by offering a comprehensive solution for the definite Chapman integral, ensuring precision over finite distances and aligning with the Chapman function at infinite limits.

Boltzmann's distribution, at a constant temperature T, describes the exponential decrease in air molecule density with altitude h, as follows:

$$\begin{array}{}\text{(1)}& n\left(h\right)=n\left(\mathrm{0}\right)e^{-\frac{mgh}{k_{\mathrm{B}}T}}.\end{array}$$

Here, m denotes the mass of a single molecule, g the gravitational acceleration, k_B the Boltzmann constant and T the absolute temperature. On a planet with radius R, the assumption of constant g is valid only when h≪R.

Considering a planet of radius R, as shown in Fig. 1, starting from point A in the atmosphere at distance D from the center, along a path at angle z from the zenith, the integral of density along the path A–B is proportional to

$$\begin{array}{}\text{(2)}& I=\underset{\mathrm{0}}{\overset{L}{\int }}{e}^{-\left(\sqrt{D^{\mathrm{2}}+l^{\mathrm{2}}+\mathrm{2}\cos zDl}-R\right)/H}\mathrm{d}l,\end{array}$$

where $L=\left|\mathrm{AB}\right|$. The H in the exponent is termed the scale height. Since the integral is performed in the atmosphere, the exponent is always negative, and the integral is well defined and is smaller than L.

Given the average density ρ of the planet, gravitational acceleration near the surface (h≪R) can be approximated as $g=\frac{\mathrm{4}\mathit{\pi }G\mathit{\rho }R}{\mathrm{3}}$, G being the gravitational constant. For an effective molecular mass m, the scale height H can be expressed as

$$\begin{array}{}\text{(3)}& H={\displaystyle \frac{k_{\mathrm{B}}T}{mg}}={\displaystyle \frac{\mathrm{3}{k}_{\mathrm{B}}T}{\mathrm{4}\mathit{\pi }G\mathit{\rho }Rm}}.\end{array}$$

Using a molecular weight W, and substituting standard values for $G=\mathrm{6.67}\times {\mathrm{10}}^{-\mathrm{11}}$ (in MKS units) and the ideal gas constant 8.31 J mol⁻¹ K⁻¹, we arrive at $H\approx \mathrm{2.97}\times {\mathrm{10}}^{\mathrm{13}}\frac{T}{\mathit{\rho }WR}$. For Earth, with $\mathit{\rho }\approx \mathrm{5.51}\times {\mathrm{10}}^{\mathrm{3}}$ kg m⁻³ and $R\approx \mathrm{6.4}\times {\mathrm{10}}^{\mathrm{6}}$ m, at a temperature of 300 K and molecular weight of 30, H is calculated to be approximately 8.5×10³ m. More pertinent to the Chapman integral is the ratio $R/H$, which for Earth is around 700. Generally, the $R/H$ ratio can be estimated as

$$\begin{array}{}\text{(4)}& R/H={\displaystyle \frac{\mathrm{4}\mathit{\pi }\mathit{\rho }{R}^{\mathrm{2}}m}{\mathrm{3}{k}_{\mathrm{B}}T}}\approx {\displaystyle \frac{\mathit{\rho }W{R}^{\mathrm{2}}}{\mathrm{2.97}\times {\mathrm{10}}^{\mathrm{13}}T}}.\end{array}$$

For rocky planets larger than 1000 km in radius and with similar density to Earth, the $R/H$ ratio is typically in the hundreds. This implies a relatively thin atmospheric layer compared to the planet's size, allowing the assumption of constant gravity as used in Eq. (6). We therefore propose changing the integration over travel distance in the atmosphere to one over the change in the radial distance.

Figure 1 Density integration from A to B.

Defining $\mathit{\lambda }=D/H$, $R_d=R/D$, $t=l/D$, and $x=L/D$, we reformulate integral I as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}I(x,z,\mathit{\lambda })& =H\mathit{\lambda }\underset{\mathrm{0}}{\overset{x}{\int }}{e}^{-\mathit{\lambda }\left(\sqrt{\mathrm{1}+t^{\mathrm{2}}+\mathrm{2}t\cos z}-R_d\right)}\mathrm{d}t\\ & =H{e}^{-\mathit{\lambda }(\mathrm{1}-R_d)}\left[\mathit{\lambda }\underset{\mathrm{0}}{\overset{x}{\int }}{e}^{-\mathit{\lambda }\left(\sqrt{\mathrm{1}+t^{\mathrm{2}}+\mathrm{2}t\cos z}-\mathrm{1}\right)}\mathrm{d}t\right].\end{array}\end{array}$$

Observing that $\mathit{\lambda }(\mathrm{1}-R_d)=(D-R)/H$, we define the term in square brackets as the definite Chapman integral, i.e., the Chapman integral with finite integration limits.

$$\begin{array}{}\text{(6)}& \text{Cd}(x,z,\mathit{\lambda })=\mathit{\lambda }\underset{\mathrm{0}}{\overset{x}{\int }}{e}^{-\mathit{\lambda }\left(\sqrt{\mathrm{1}+t^{\mathrm{2}}+\mathrm{2}t\cos z}-\mathrm{1}\right)}\mathrm{d}t\end{array}$$

Specifically, we identify

$$\begin{array}{}\text{(7)}& \text{Cd}(\mathrm{\infty },z,\mathit{\lambda })\equiv \text{Ch}(\mathit{\lambda },z),\end{array}$$

where Ch(λ,z) is the Chapman function as defined in Chapman (1931).

2 Analytical solution of the definite Chapman integral

To perform the integral Cd in Eq. (6), we make the following change of variable:

$$\begin{array}{}\text{(8)}& u\left(t\right)=\sqrt{\mathrm{1}+t^{\mathrm{2}}+\mathrm{2}t\cos z}-\mathrm{1}.\end{array}$$

Restricting z to $[\mathrm{0},\mathit{\pi }/\mathrm{2}]$, there is a one-to-one mapping between t and u. Using the relationship $\frac{\mathrm{d}t}{\mathrm{d}u}=\frac{u+\mathrm{1}}{\sqrt{u^{\mathrm{2}}+\mathrm{2}u+{\cos }^{\mathrm{2}}z}}$, the integral is transformed to

$$\begin{array}{}\text{(9)}& \text{Cd}(x,z,\mathit{\lambda })=\mathit{\lambda }\underset{\mathrm{0}}{\overset{y=u\left(x\right)}{\int }}{\displaystyle \frac{\mathrm{1}+u}{\sqrt{(\mathrm{1}+u{)}^{\mathrm{2}}-{\sin }^{\mathrm{2}}z}}}{e}^{-\mathit{\lambda }u}\mathrm{d}u.\end{array}$$

Geometrically, referring to Fig. 1, the upper limit of the integral $y=\left|\mathrm{OB}\right|/D-\mathrm{1}$, denoting the variation in radial distance at the boundary of integration, normalized by the initial distance. Given that the atmospheric thickness of a planet is significantly less than its radius, the upper limit y is much smaller than 1 ($u\le y\ll \mathrm{1}$). With the range of z confined to $[\mathrm{0},\mathit{\pi }/\mathrm{2}]$, sin z is non-negative, this permits the approximation $\sqrt{\mathrm{1}+u+\sin z}\approx \sqrt{\mathrm{1}+\sin z}$. With this simplification, the above integral can be approximated as follows:

$$\begin{array}{}\text{(10)}& \text{Cd}(x,z,\mathit{\lambda })\approx {\displaystyle \frac{\mathit{\lambda }}{\sqrt{\mathrm{1}+\sin z}}}\underset{\mathrm{0}}{\overset{y}{\int }}{\displaystyle \frac{\mathrm{1}+u}{\sqrt{\mathrm{1}+u-\sin z}}}{e}^{-\mathit{\lambda }u}\mathrm{d}u.\end{array}$$

Since λ is large, the main contribution to the integral comes from small u values. Moreover, the assumption of constant gravity in deriving the exponential drop of density is valid only when the atmosphere depth is much smaller than planet radius. These considerations further justify our approximation.

By another change of variable, $w=\mathrm{1}+u-\sin z$, then integrate by parts, the integral can be analytically expressed using the erfc(t) function (complementary error function), which is defined as $\mathrm{1}-\mathrm{2}/\mathit{\pi }{\int }_{\mathrm{0}}^t\exp (-u^{\mathrm{2}})\mathrm{d}u$.

To simplify the expression of our result, we define function $Y(y,z,\mathit{\lambda })$ for $z\in [\mathrm{0},\mathit{\pi }/\mathrm{2}]$,

$$\begin{array}{}\text{(11)}& \begin{array}{rl}Y\left(y,z,\mathit{\lambda }\right)& \equiv {\displaystyle \frac{-\mathrm{1}}{\sqrt{\mathit{\lambda }\left(\mathrm{1}+\sin z\right)}}}[e^{-\mathit{\lambda }y}\sqrt{\mathit{\lambda }\left(\mathrm{1}+y-\sin z\right)}\\ & +\sqrt{\mathit{\pi }}{e}^{\mathit{\lambda }(\mathrm{1}-\sin z)}\left(\mathit{\lambda }\sin z+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)\\ & \text{erfc}\left(\sqrt{\mathit{\lambda }\left(\mathrm{1}+y-\sin z\right)}\right)],\end{array}\end{array}$$

and define function $\text{Cy}(y,z,\mathit{\lambda })$ as

$$\begin{array}{}\text{(12)}& \text{Cy}(y,z,\mathit{\lambda })=Y(y,z,\mathit{\lambda })-Y(\mathrm{0},z,\mathit{\lambda }).\end{array}$$

The definite Chapman integral is found to be

$$\begin{array}{}\text{(13)}& \text{Cd}(x,z,\mathit{\lambda })=\text{Cy}(y,z,\mathit{\lambda }),\end{array}$$

where y is defined by Eq. (8); i.e., $y=\sqrt{\mathrm{1}+x^{\mathrm{2}}+\mathrm{2}x\cos z}-\mathrm{1}$. Geometrically, $y=D_{\mathrm{B}}/D-\mathrm{1}$, D_B being the distance from the end point to the center of the planet.

To study the behavior of $\text{Cy}(y,z,\mathit{\lambda })$, we examine its first derivative:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{d}\text{Cy}(y,z,\mathit{\lambda })}{\mathrm{d}y}}={\displaystyle \frac{\mathit{\lambda }(\mathrm{1}+y)e^{-\mathit{\lambda }y}}{\sqrt{(\mathrm{1}+\sin z)(\mathrm{1}+y-\sin z)}}}.\end{array}$$

Since $\mathrm{d}Y/\mathrm{d}y$ is always positive, $Y(y,z,\mathit{\lambda })$ increases monotonically with y. Moreover, due to the factor e^−λy, the derivative quickly approaches 0 at λy≫1. This indicates that the integral's primary contribution comes from within a few multiples of the scale height, while the contribution from higher altitudes becomes inconsequential. For instance,with λ=500, $Y(y,z,\mathit{\lambda })$ plateaus around y≈0.02. Consequently, $\text{Cy}(y=\mathrm{10}/\mathit{\lambda },z,\mathit{\lambda })$ serves as an excellent approximation for the Chapman function, despite the latter having the integration limit extended to infinity. Our results (Eqs. 11–13) agree with the approximate formulas tabulated in Vasylyev (2021) when evaluated under appropriate limits.

Our result is an analytical solution for the definite Chapman integral applicable to zenith angles z restricted to the range $[\mathrm{0},\mathit{\pi }/\mathrm{2}]$. In this context, y must be positive. However, our solution can be easily extended to situations where $z>\mathit{\pi }/\mathrm{2}$, involving a decrease in radial distance along the integration path. In the simplest scenario, reversing the start and end points of the integration ensures that the zenith angle is $z<=\mathit{\pi }/\mathrm{2}$ at the starting point. For such cases, it is merely a matter of redefining D and z based on the new starting point.

Figure 2 depicts a more intricate scenario, where the zenith angles at both integration ends exceed $\mathit{\pi }/\mathrm{2}$.

Figure 2 Illustration of the density integral from point A to B, with zenith angles greater than $\mathit{\pi }/\mathrm{2}$ at both endpoints. The integral is divided into two parts at point A^′, enabling the application of the Y function.

To adapt the $Y(y,z,\mathit{\lambda })$ function to the case illustrated in Fig. 2, the approach involves altering the integration's starting point. This is achieved by drawing a perpendicular line from the center of the planet to the line AB and taking the intersection point A^′ as the new starting point. The Y function is then applied to two segments, from A^′ to A and from A^′ to B, both with a zenith angle of $\mathit{\pi }/\mathrm{2}$.

With this change of the starting point, let $D^{\prime }=D\sin z$, ${\mathit{\lambda }}^{\prime }=D\sin z/H$, $y_{\mathrm{1}}=\mathrm{1}/\sin z-\mathrm{1}$ and $y_{\mathrm{2}}=\left|\mathrm{OB}\right|/D^{\prime }-\mathrm{1}$:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}\text{Cd}(x,z,\mathit{\lambda })& =Y(y_{\mathrm{1}},\mathit{\pi }/\mathrm{2},\mathit{\lambda }\sin z)+Y(y_{\mathrm{2}},\mathit{\pi }/\mathrm{2},\mathit{\lambda }\sin z)\\ & -\mathrm{2}Y(\mathrm{0},\mathit{\pi }/\mathrm{2},\mathit{\lambda }\sin z).\end{array}\end{array}$$

Figure 3 Comparison of the analytical result and numerical integration.

Figure 4 Relative error of the analytical result compared to numerical integration.

3 Asymptotic expression

Given that λ is significantly greater than 1, the erfc function values in Eq. (11) rapidly converge to 0 at both limits for most z values (for instance, $\text{erfc}\left(\mathrm{3}\right)\approx \mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{5}}$). Simultaneously, the exponential factor of the second term inside the equation's brackets becomes exceedingly large for most z values. As previously mentioned, the original integral remains well defined and is smaller than the length of integration. Therefore, Eq. (12) is dependent on the near cancellation of erfc values at the integration limits:

$$\begin{array}{}\text{(16)}& \mathrm{\Delta }(y,z,\mathit{\lambda })=\text{erfc}\left(\sqrt{\mathit{\lambda }(\mathrm{1}-\sin z)}\right)-\text{erfc}\left(\sqrt{\mathit{\lambda }(\mathrm{1}+y-\sin z)}\right).\end{array}$$

For high values of λ, attempting a direct numerical calculation using Eqs. (11)–(12) could lead to overflow issues with the exponential term and imprecise results in the Δ term, due to the limitations in floating-point precision. It is crucial to analytically neutralize the positive exponent in the second term of Eq. (11). When $\mathit{\lambda }(\mathrm{1}-\sin z)\gg \mathrm{1}$, by retaining only the principal term in the asymptotic expansion of erfc(x), namely $\exp (-x^{\mathrm{2}})/x\sqrt{\mathit{\pi }}$, we can simplify the Δ expression:

$$\begin{array}{}\text{(17)}& \mathrm{\Delta }(y,z,\mathit{\lambda })={\displaystyle \frac{e^{-\mathit{\lambda }(\mathrm{1}-\sin z)}}{\sqrt{\mathit{\pi }}\sqrt{\mathit{\lambda }(\mathrm{1}-\sin z)}}}\left(\mathrm{1}-e^{-\mathit{\lambda }y}\sqrt{{\displaystyle \frac{\mathrm{1}-\sin z}{y+\mathrm{1}-\sin z}}}\right).\end{array}$$

Using the above result, at large λ(1−sin z), Eq. (13) becomes

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\text{Cd}(x,z,\mathit{\lambda })& \approx {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathit{\lambda }(\mathrm{1}+\sin z)}}}[\sqrt{\mathit{\lambda }(\mathrm{1}-\sin z)}\\ & -e^{-\mathit{\lambda }y}\sqrt{\mathit{\lambda }(y+\mathrm{1}-\sin z)}\\ & +{\displaystyle \frac{\left(\mathit{\lambda }\sin z+\frac{\mathrm{1}}{\mathrm{2}}\right)}{\sqrt{\mathit{\lambda }(\mathrm{1}-\sin z)}}}\left(\mathrm{1}-e^{-\mathit{\lambda }y}\sqrt{{\displaystyle \frac{\mathrm{1}-\sin z}{y+\mathrm{1}-\sin z}}}\right)]\end{array},\end{array}$$

and when λy≫1, the exponentially small terms in $\text{Cd}(x,z,\mathit{\lambda })$ above can be dropped. The formula is reduced to the well-known result in the limiting case.

$$\begin{array}{}\text{(19)}& \text{Cd}(x,z,\mathit{\lambda })\approx {\displaystyle \frac{\mathrm{1}+\frac{\mathrm{1}}{\mathrm{2}\mathit{\lambda }}}{\cos z}}\approx {\displaystyle \frac{\mathrm{1}}{\cos z}}\end{array}$$

It is important to observe that this approximation holds true only when $\mathit{\lambda }(\mathrm{1}-\sin z)\gg \mathrm{1}$, and cos z is non-zero at this limit. This indicates that for small zenith angles, the atmospheric curvature can be disregarded, and the optical depth calculations can be based simply on the length of the slanted path.

4 Numerical evaluation

The sole approximation in our derivation was made in Eq. (10). Our analytical results, spanning Eqs. (11)–(15), are valid for any zenith angle, including $z={\mathrm{90}}^{\circ }$. To evaluate our solution, we compared the analytical results from $Y(y,z,\mathit{\lambda })$ (Eqs. 11–12) with direct numerical integration of the original integral $\text{Cd}(x,z,\mathit{\lambda })$ in Eq. (6), across a range of λ values and zenith angles within $[\mathrm{0},\mathit{\pi }/\mathrm{2}]$. Then we plotted the relative error of our analytical solution, calculated as the discrepancy between the analytical and numerical results, normalized by the numerical integral. The full evaluation is demonstrated in the GitHub repository (Yue, 2024). The key resulting plots are presented in Figs. 3 and 4.

Our numerical evaluations revealed that the maximum relative error in the analytical solution remained under 0.5 % for λ values ranging from 50 to 10 000. Furthermore, the asymptotic approximation in Eq. (18) demonstrates high accuracy, with the maximum relative error of less than 1 % when it is applied at $\sqrt{\mathit{\lambda }(\mathrm{1}-\sin z)}>\mathrm{7.0}$. Even when the asymptotic approximation is switched on at $\sqrt{\mathit{\lambda }(\mathrm{1}-\sin z)}>\mathrm{3.0}$, the relative error stayed below 5 %.

Additionally, we assessed our analytical results at an upper limit of y=0.1 for λ values between 50 and 10 000, juxtaposing them with the numerical values of the Chapman function. The comparisons indicated that they are within 0.5 % of each other.

5 Conclusions

In summary, our study provides a comprehensive analytical solution for the definite Chapman integral, applicable to any zenith angle and realistic λ values. The accuracy of our solution has been rigorously tested against direct numerical integrations, demonstrating a high degree of precision with relative errors consistently below 0.5 %. The solution is notable in its simplicity and versatility. This work paves the way for more efficient and accurate atmospheric effect analyses and related studies.