Efficient Computation of Collisional ℓ-mixing Rate Coefficients in Astrophysical Plasmas

Energy-conserving angular momentum-changing $n{\ell }\to n{\ell }^{\prime}$ transitions induced in collisions between Rydberg atoms and low velocity ions are needed for accurate comparison between astrophysical observations and models that employ atomic theory for temperature and density diagnostics in diffuse atomic clouds, H ii regions, and various nebulae. With ion collision-induced angular momentum mixing rate coefficients scaling as n⁴, values as large as a few times 10⁵ cm³ s⁻¹ are possible for principal quantum numbers near n ∼ 200. Accurate and efficiently calculated rate coefficients are hence necessary to interpret a host of astrophysical processes, such as radio recombination lines from hydrogen (HRRL) and carbon (CRRL) as tracers of the neutral phase of the interstellar medium (ISM; e.g., Oonk et al. 2014, 2017; Salas et al. 2018), and from hydrogen as a tracer of gas ionized by young stars (H ii regions; e.g., Roelfsema & Goss 1992; Anderson et al. 2011). The recombination of hydrogen and helium in the early universe, and the primordial abundance of helium, are also examples of processes affected by collision physics (Chluba & Sunyaev 2006; Izotov & Thuan 2010).

In their pioneering work, Pengelly and Seaton obtained proton–Rydberg hydrogen collisional cross sections for dipole allowed transitions within the Born–Bethe approximation (Pengelly & Seaton 1964, hereafter PS64). Given that the probability for Δℓ = ±1 transition falls off asymptotically as the inverse square of the impact parameter and the cross section hence becomes logarithmically divergent, PS64 invoked a set of cutoff conditions to circumvent this divergence. The arbitrariness implicit in choosing these cutoff conditions was reexamined in Vrinceanu & Flannery (2001a, 2001b, hereafter VF01a and VF01b), where a non-perturbative closed form solution for the transition probability was found.

A semiclassical (SC) rate coefficient for arbitrary ℓ-changing collisions was derived in Vrinceanu et al. 2012 (VOS12), and was shown to be in agreement with both classical trajectory Monte Carlo simulations and numerically integrated quantum rate coefficients for transitions with $| {\rm{\Delta }}{\ell }| \gt 1$. For the dipole allowed transitions $(| {\rm{\Delta }}\ell | =1$), SC probabilities evaluated by VOS12 increase linearly with impact parameter and at a critical value abruptly vanish. This unphysical behavior of the SC transition probability, first noticed in Storey & Sochi (2015), was addressed by the modified PS64 (PS-M) approximation (Guzmán et al. 2016, 2017). A further improved SC probability with the correct asymptotic behavior, yielding a more accurate formula for the dipole transitions, was later derived by Vrinceanu et al. (2017, hereafter VOS17).

In this work, we derive a computationally efficient formula of the rate coefficients for collision-induced dipole transitions, accurate for a broad range of n, temperatures (T), and electron number densities (n_e). This is achieved by using a closed-form expression for the dipole transition rate coefficients, therefore overcoming the need for the explicit calculation of cross sections. In addition, this formulation circumvents an unphysical behavior of the PS64 rate coefficients, which become negative for a range of astrophysically relevant parameters. The resulting rate coefficients are then used to evaluate the departure coefficients from the statistical equilibrium of HRRL.

The rate coefficients for n-conserving, ℓ-changing transitions are obtained by integrating the corresponding transition probability, ${P}_{n{\ell }\to n{\ell }^{\prime} }$, over the impact parameter, b, and thermal distribution of the projectile velocity, v,

$$\begin{eqnarray}&&\begin{array}{l}k(n,{\ell },{\ell }^{\prime} ,T,{R}_{c})\\ \quad =\,2\pi {\int }_{0}^{\infty }{f}_{{MB}}(v)\,{vdv}{\int }_{0}^{{R}_{c}}{P}_{n{\ell }\to n{\ell }^{\prime} }(b,v)\,{bdb},\end{array}\end{eqnarray} \tag{ 1 }$$

with f_MB as the Maxwell–Boltzmann distribution at temperature T. The cutoff distance R_c is required to regularize the divergent integral for dipole allowed transitions ($| {\rm{\Delta }}{\ell }| =1$), when the transition probability decreases too slowly for $b\to \infty$. For all other cases, the integral is finite as ${R}_{c}\to \infty$. The transition probability does not depend on b and v independently, but through the collision parameter $\alpha \,=3{Zn}{\hslash }/(2\,{mvb})$ (VF01a), so that the double integral in Equation (1) is reduced to

$$\begin{eqnarray}&&\begin{array}{l}k(n,{\ell },{\ell }^{\prime} ,T,{R}_{c})={n}^{4}{a}_{0}^{2}{v}_{0}\sqrt{\displaystyle \frac{8\pi \mu {v}_{0}^{2}}{{k}_{B}T}}\\ \quad \times \,{\int }_{0}^{\infty }{{zP}}_{n{\ell }\to n{\ell }^{\prime} }(z)\,{e}^{-\theta {z}^{2}/2}\,{dz},\end{array}\end{eqnarray} \tag{ 2 }$$

where ${a}_{0}=5.29177\times {10}^{-11}$ m is the Bohr radius, ${v}_{0}\,=2.18769\times {10}^{6}$ m s⁻¹ is the atomic unit of velocity, and the integration variable is z = 3/(2n α). The parameter θ, small for large R_c, is defined as

$$\begin{eqnarray}&&\theta ={n}^{4}\,\displaystyle \frac{\mu {v}_{0}^{2}}{{k}_{B}T}\displaystyle \frac{{a}_{0}^{2}}{{R}_{c}^{2}},\end{eqnarray} \tag{ 3 }$$

where μ is the reduced mass of the projectile—target system.

When R_c is chosen to be the Debye length ${\lambda }_{D}=\sqrt{{\epsilon }_{0}{k}_{B}T/({n}_{e}{e}^{2})}$, as in PS64, this parameter becomes

$$\begin{eqnarray}&&\theta =1.704675\times {10}^{-10}{n}^{4}\displaystyle \frac{{n}_{e}}{{T}^{2}},\end{eqnarray} \tag{ 4 }$$

with T in K and n_e in cm⁻³. Depending on the specific physical situations, other choices for R_c are possible as discussed in PS64, changing the θ parameter accordingly.

Equation (2) is the starting point for the calculation of the rate coefficient using various approaches depending on the choice of ${P}_{n{\ell }\to n{\ell }^{\prime} }$. The relationships among these approximations reflect the organization of the software package (Vrinceanu 2018), and are illustrated in the diagram in Figure 1, and discussed in detail below.

Zoom In Zoom Out Reset image size

The rate coefficients are calculated by numerically integrating Equation (2). This integral does not pose difficulties for relatively low principal quantum numbers n ≲ 40 when using the VF01b non-perturbative quantum mechanical ${P}_{n{\ell }\to n{\ell }^{\prime} }$. However, for larger n, the calculation does not converge because of truncation errors and near cancellation of large terms. These problems are addressed in our code by using exact number arithmetic for large factorials and extended floating point precision for large order polynomials. The computational time increases as ∼ℓ × n, to the point that accurate calculations become impractical; for example, the calculation for n = 1000 can take up to 5 minutes. More efficient, but less accurate, approximations are implemented in the code for large n and ℓ, when numerical integration of quantum probability is slow. The method parameter in the code selects the desired procedure, either quantum, or from the menu of approximations discussed below, and should be chosen based on the accepted level in the trade-off of accuracy versus computing cost.

The VOS12 SC formula works well for moderate to large $| {\rm{\Delta }}{\ell }|$, but loses accuracy for larger $| {\rm{\Delta }}{\ell }|$, for which the rate coefficients are small and may be neglected. The VOS12 SC formula has been derived by integrating the classical limit of VF01a transition probability, and has been validated by extensive classical trajectory Monte Carlo simulations. The SC approximation fails for dipole transitions $| {\rm{\Delta }}{\ell }| =1$, because the SC transition probability increases linearly with the impact parameter and drops abruptly after a critical value, instead of decreasing as b⁻², as obtained with perturbation theory (PS64).

For dipole transitions, the code provides rate coefficients for combined ${\ell }\to {\ell }+1$ and ${\ell }\to {\ell }-1$ transitions. In the Born approximation, ${P}_{n{\ell }\to n{\ell }^{\prime} }$ is assumed to be 1/2 for R < R₁, after which it decreases as 1/b². The PS64 formula is derived from this approximation by adopting the additional assumption that R₁ <R_c. This assumption fails for large n_e and low T, limiting the range of applicability of PS64 to ${n}^{2}({n}^{2}-{{\ell }}^{2}-{\ell }-1){n}_{e}/{T}^{2}\lt 2.98\,\times {10}^{9}$ cm${}^{-3}$/K${}^{2}$. Within this range of parameters, PS64 is in reasonable agreement with the non-perturbative quantum results, as demonstrated in the next section. Outside this range, PS64 rate coefficients become negative (Guzmán et al. 2017; Salgado et al. 2017).

The PS-M approximation introduced in Guzmán et al. (2017) replaces the constant 1/2 transition probability in the PS64 model with a linearly increasing one, and obtains the rate coefficient by averaging the resulting cross section over all energies, including those neglected by PS64 for R₁ < R_c. The approximate PS-M rate coefficients shown in table I in Guzmán et al. (2017) are positive even when PS64 are negative, and yield an overall better agreement with the quantum results, although significant deviations are noted in some cases, up to a factor of 10, probably due to the simplicity of the model.

The shortcomings of the VOS12 SC approximation for dipole allowed transitions were addressed in VOS17 by deriving a more accurate SC probability with the correct large b asymptotics. When used in Equation (2), the VOS17 transition probability leads to an SC rate coefficient for dipole allowed transitions

$$\begin{eqnarray}&&\begin{array}{l}k(n,{\ell },T,{n}_{e})=\sqrt{\displaystyle \frac{\pi }{2}}{a}_{0}^{2}{v}_{0}\sqrt{\displaystyle \frac{\mu {v}_{0}^{2}}{{k}_{B}T}}{D}_{n{\ell }}\left[\displaystyle \frac{3\sqrt{\pi }}{4{x}^{3/2}}\,\mathrm{erf}(\eta \sqrt{x})\right.\\ \quad \left.-\displaystyle \frac{3\eta }{2x}{e}^{-{\eta }^{2}x}+\displaystyle \sum _{k=0}^{{N}_{t}}{A}_{k}({B}_{k}-3\gamma -\mathrm{log}(4x)){x}^{k}\right],\end{array}\end{eqnarray} \tag{ 5 }$$

where erf is the error function, ${D}_{n{\ell }}=6{n}^{2}({n}^{2}-{{\ell }}^{2}-{\ell }-1)$, η = 0.277855 is the solution of the equation ${j}_{1}{\left(1/z\right)}^{2}\,=z/6$, with j₁(u) the spherical Bessel function of order 1 and argument u, and $x=3{D}_{n{\ell }}\theta /(4{n}^{4})={n}_{e}{D}_{n{\ell }}/(7.82162\times {10}^{9}{T}^{2})$.

The coefficients in the asymptotic expansion can be calculated up to the truncation order N_t by using

$$\begin{eqnarray*}&&{A}_{k}=\displaystyle \frac{18\times {4}^{k}}{(k+2)(k+3)(2k+3)\,k!\,(2k+1)!}\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{B}_{k}=\displaystyle \frac{1}{k+2}+\displaystyle \frac{1}{k+3}+\displaystyle \frac{2}{2k+3}+{H}_{k}+2{H}_{2k+1},\end{eqnarray*}$$

where ${H}_{k}={\sum }_{j=1}^{k}1/j$ for k > 0, H₀ = 0, is k the kth harmonic number, and γ = 0.57721 the Euler constant. A derivation of this formula, and a list of the first eleven A_k and B_k coefficients, are given in the Appendix. Our experience shows that a truncation order N_t ∼ 16 − 20 is sufficient to provide accurate results.

As illustrated in Figure 1, the Python module Lmixing (Vrinceanu 2018) has a function that calculates the angular momentum mixing rate coefficient. The required arguments are the principal quantum number n, the initial angular momentum ℓ, and the temperature in Kelvin. The optional arguments are the electron number density n_e in cm⁻³, the final angular momentum ${\ell }^{\prime}$, and the method of calculation. The choices for this argument are: quantum (default value) using Equation (2) with the probability defined by Equation (3) in VF01b, semiclassical using Equation (5), Born for the Born approximation Equation (8), PS-M for the approximation in Equation (18) originally introduced in Guzmán et al. (2017), and P_and_S which implements the PS64 method. If ${\ell }^{\prime}$ is not provided, the rate for the combined ${\ell }\to {\ell }\pm 1$ transitions is calculated. If ${\ell }^{\prime}$ is given, one can then choose the exact quantum calculation, or the SC approximation (Equation (8) in VOS12) with the value for the method parameter classical, provided that $| {\rm{\Delta }}{\ell }| \gt 1$.

Figure 2 shows a comparison between the rate coefficients coefficients for dipole allowed transitions for ℓ = 1 and ℓ = n − 2 angular momenta, calculated by integrating the VF01b quantum formula, the PS64 approximation and the SC approximation in Equation (5). The evaluation of the quantum case for extremely large n is slow even when low accuracy results are sufficient, due to the use of multiprecision floating point arithmetic necessary to prevent truncation error and loss in precision. The calculations in Figure 2 required 400 digits of precision and took several hours for ℓ = 1 cases and two days for ℓ = n − 2 cases to complete on a single processor. At around n = 200, in the ℓ = 1 case, the PS64 approximation fails and becomes negative for higher n, as first discussed in Salgado et al. (2017), and in the caption of Table 1 in Guzmán et al. (2017), while the SC approximation is in good agreement with the quantum results over a much larger range of n. For sufficiently low T and high n_e, the SC rate coefficients for n ≳ 500 overestimate the corresponding quantum result, then underestimate the latter at even higher n (not shown in Figure 2), eventually becoming negative. However, unlike the case of PS64, the negative-defined region occurs for progressively higher values of n as more terms in the expansion of Equation (5) are considered. In the ℓ = n − 2 case, the PS64 results are more reliable with respect to the corresponding ones at low angular momentum, and begin to diverge from the exact results for n > 700, becoming negative at around n = 1000. The overall shape of the dependence of the rate coefficient k with n is the result of two competing factors, as seen in Equation (5): on the one hand the prefactor in front of the integral increases as n⁴, while the integral decreases roughly as 1/θ. For comparison, we also include in Figure 2 the PS-M rate coefficients which behave more consistently with n than the PS64 model in the ℓ = 1 case, and approach the quantum results at large n. However, other than the large n limit for ℓ = 1, PS-M predicts rate coefficients that are consistently smaller than QM and SC values, as particularly evident in the ℓ = n − 2 case.

Zoom In Zoom Out Reset image size

In Figure 3, the differences of the PS64 approximation (left panel) and of the SC approximation (right panel), both relative to the quantum rate coefficients, are shown for astrophysically relevant values of n, T, and n_e, and for the ${\ell }=1$ transition. The SC approximation is accurate within 1% over a wide range of T and n, while the accuracy of PS64 is roughly one order of magnitude worse at the same point in the $n-{T}^{2}/{n}_{e}$ plane. The PS64 approximation yields negative rate coefficients in the upper left corner of the diagram, in the low T, high n_e, and high n regime. As already discussed for Figure 2, the SC rate coefficients also become negative in the upper left corner, but in comparison to the PS64 case this occurs for a smaller region with size inversely proportional to the number of terms N_t in Equation (5). These cases are extreme, but important for the interpretation of HRRLs that probe the low electron density, cool ISM (e.g., Salgado et al. 2017 and references therein). The SC approximation obviously becomes less accurate in the region for which the PS64 rate coefficients approach negative values. This is also expected because for these parameters ${\lambda }_{D}\leqslant {n}^{2}{a}_{0}$, indicating that the binary collision assumption may fail, with the ℓ-mixing instead dominated by many-body interactions. If the opposite case of ℓ = n − 2 is considered, the rate coefficients for PS64 and the SC approximations are much closer, in line with what is expected by inspecting the corresponding curves in Figure 2. Figure 3 puts on a more quantitative standing the recent debate on the accuracy of various proposed rate coefficients as reported in Storey & Sochi (2015), Guzmán et al. (2016, 2017), and Williams et al. (2017). The SC approximation is in general more accurate than the above approximations because the VOS17 transition probability agrees better with the quantum results. It was noticed in Guzmán et al. (2017) that for some extreme cases, the PS-M overestimates the quantum results by a factor of 10.

Zoom In Zoom Out Reset image size

The HRRL results shown in Figure 4 were produced with the Salgado et al. (2017) models using the updated ℓ-changing collision rate coefficients, Equation (5) with N_t = 10 and, as a comparison, also using the PS64 rate coefficients. The latter are known to be in good agreement with the quantum mechanical rate coefficients for sufficiently high T, high n_e, and low n.

Zoom In Zoom Out Reset image size

We show in Figure 4 the results for the departure coefficients from a thermal population, b_n, for a homogeneous gas slab with (i) T_e = 100 K and n_e = 0.1 cm⁻³, and (ii) T_e = 10 K and n_e = 100 cm⁻³. The former is a typical cool ISM case and the latter is an an extreme case. Both are exposed to a Galactic power-law radiation field ${T}_{R}\propto {\lambda }^{2.6}$ that is normalized at 100 MHz by T_R,100 = 2000 K (Salgado et al. 2017). We find that the different ℓ-changing collision rate coefficients primarily affect b_n at low to intermediate n values (n ≲ 300) with differences up to a few percent for typical cool ISM conditions. As explained in Salgado et al. (2017), but see also Hummer & Storey (1987), this is because at these intermediate n levels collisions compete with spontaneous decay, effectively storing electrons in high ℓ sublevels for which radiative decay is less important. Notably, the ℓ-changing collision rate coefficients presented here allow us to efficiently calculate b_n values for high n, where the PS64 rate coefficients no longer apply, and which are important to studies of cool, partially ionized ISM (e.g., Oonk et al. 2017; Salas et al. 2018).

The HRRL optical depths are calculated using the product of b_n β_n, where the correction factor for stimulated emission β_n can be seen as the derivative of b_n (Salgado et al. 2017), such that small changes in b_n can lead to somewhat larger changes in HRRL optical depth. These changes are measurable, but require very high signal-to-noise observations across a broad frequency range (i.e., 240–2000 MHz). Most HRRL observations have difficulties achieving such accuracy, and hence the calculated differences will be within current observational uncertainties for typical ISM conditions. We have also compared the HRRL results presented here with those computed using the VOS12 SC rate coefficients (Salgado et al. 2017). We find that for these cases the results agree to within a few percent for n ≳ 300. For lower n-values the agreement is less good.

We are currently implementing our new ℓ-changing rate coefficients for CRRLs, as will be presented in a future work. Although we anticipate that the results may be qualitatively similar to those for HRRLs, from a quantitative standpoint we expect to have stronger dependence of the departure coefficients on the ℓ-changing collision rate coefficients than HRRLs (Salgado et al. 2017).

We have introduced an efficient and accurate SC approximation for ℓ-mixing processes, allowing direct evaluation of collisional rate coefficients for a broad range of n, T, and n_e relevant to astrophysical plasmas. We provide a Python module for evaluating these rate coefficients for easy integration in large-scale radiative-collisional simulation codes. The Python code also implements a multiprecision numerical integration of the quantum rate constants for small to moderate n. The relationship and range of validity and accuracy for various schemes to evaluate ℓ-changing rate coefficients are elucidated. Furthermore, we identify the range of parameters for which PS64 rate coefficients become unphysical. Efficient computer codes for the evaluation of ℓ-mixing rate coefficients with various levels of approximation are beneficial for accurate astrophysical modeling. As noted in the recent release of Cloudy (Ferland et al. 2017), different choices for the rate coefficients can lead to differences up to 10% in the predicted line intensities, which are larger than the precision of current observations. The rate coefficients proposed here are in better agreement with the more rigorous, but computationally less efficient quantum rate coefficients.

This work was supported by the National Science Foundation through a grant to ITAMP at the Harvard-Smithsonian Center for Astrophysics. One of the authors (D. V.) is also grateful for the support received from the National Science Foundation through grants PHY-1831977 and HRD-1829184. J. B. R. O. and P. S. acknowledge financial support from the Dutch Science Organization (NWO) through TOP grant 614.001.351.

Software: Lmixing (Vrinceanu 2018).

This section presents the derivations of the main results and implementation details of the computational module. According to the organization of the code illustrated in Figure 1 we will discuss the four main sections, (i) the integration of the exact quantum formula, (ii) the Born approximation, (iii) the classical approximation, and (iv) the SC approximation. These approaches are derived from Equation (2) either by keeping the exact, but computationally expensive form for the transition probability ${P}_{n{\ell }\to n{\ell }^{\prime} }$, or by replacing it by an approximate expression that provides quicker results with a limited range of validity. The goal is to calculate the dimensionless integral in Equation (2)

$$\begin{eqnarray}&&\tilde{k}(n,{\ell },{\ell }^{\prime} ,\theta )={\int }_{0}^{\infty }{{zP}}_{n{\ell }\to n{\ell }^{\prime} }(z){e}^{-\theta {z}^{2}/2}\,{dz}.\end{eqnarray} \tag{ 6 }$$

The exact calculation uses the quantum transition probability derived in (VF01b), which has practical use limited to small (n < 30) quantum numbers in regular computer arithmetic. The main two reasons for this difficulty are the combinatorial Wigner 6-j symbols involving factorials of large integers that cannot be represented exactly with integer types, and the loss of precision in the calculation of polynomials of large order with alternating terms. Our code takes advantage of the unlimited size integer type in Python by precomputing tables of large factorials, and by using fixed-point representations of real numbers with a prescribed, but arbitrary, number of decimal figures. As a rule of thumb, we find that a calculation for quantum numbers n requires setting the precision at n/2 digits. The numerical integral Equation (6) is calculated with a recurrent Gauss–Lobatto-Kronrod algorithm from Press et al. (1992) adapted for the use of extended precision real numbers.

A Born approximation is obtained from perturbation theory (VOS17) only for dipole allowed transitions. The cumulative probability for both ${\ell }\to {\ell }^{\prime}$ transitions is approximated as

$$\begin{eqnarray}{P}^{(B)}=\displaystyle \frac{1}{2}\left\{\begin{array}{ll}1, & \mathrm{if}\ \ z\leqslant \sqrt{{D}_{n{\ell }}/{n}^{4}}\\ {D}_{n{\ell }}/({n}^{4}{z}^{2}), & \mathrm{if}\ \ z\gt \sqrt{{D}_{n{\ell }}/{n}^{4}}\end{array}\right.,\end{eqnarray} \tag{ 7 }$$

with ${D}_{n{\ell }}=6{n}^{2}({n}^{2}-{{\ell }}^{2}-{\ell }-1)$. When used in Equation (6), this transition probability provides the Born approximation for the rate

$$\begin{eqnarray}&&{\tilde{k}}^{(B)}(n,{\ell },\theta )=\displaystyle \frac{1-{e}^{{D}_{n{\ell }}\theta /2{n}^{4}}}{2\theta }+\displaystyle \frac{{D}_{n{\ell }}}{4{n}^{4}}{\rm{\Gamma }}(0,{D}_{n{\ell }}\theta /2{n}^{4}),\end{eqnarray} \tag{ 8 }$$

where Γ(0, x) is the incomplete gamma function. The PS64 result was derived under the same conditions with the additional tacit assumption that the cutoff impact parameter R_c is greater than the R₁ parameter (R_c > R₁) for any projectile speed. Since the transition impact parameter increases with the speed of the projectile as ∼1/v, the thermal average Equation (1) will have a contribution from small speeds for which R₁ > R_c. This contribution, neglected in PS64, diminishes as ${R}_{c}\to \infty$. Indeed, the PS64 rate is obtained from Equation (8) in the $\theta \to 0$ limit

$$\begin{eqnarray}&&{\tilde{k}}^{({PS})}(n,{\ell },\theta )=\displaystyle \frac{{D}_{n{\ell }}}{4{n}^{4}}\left[1-\gamma -\mathrm{log}({D}_{n{\ell }}\theta /2{n}^{4})\right],\end{eqnarray} \tag{ 9 }$$

where γ is the Euler constant. This corresponds to Equation (43) in PS64. The rate coefficients in Equation (9) become negative for n large enough such that ${D}_{n{\ell }}\theta /{n}^{4}\gt 2{e}^{1-\gamma }$. In other words, the PS64 approximation is limited to cases where ${n}^{2}({n}^{2}-{{\ell }}^{2}-{\ell }-1){n}_{e}/{T}^{2}\lt 2.98\,\times \,{10}^{9}$, which for small ℓ reduces roughly to ${n}^{4}{n}_{e}/{T}^{2}\lt 3\times {10}^{9}$, and for large ℓ to ${n}^{3}{n}_{e}/{T}^{2}\lt {10}^{9}$.

An SC transition probability, compatible with the Born approximation and increasing linearly for small impact parameter, as predicted by the classical approximation, was obtained in (VOS17). In terms of the parameter z, the unresolved transition probability, summed over the ${\ell }^{\prime}$ and slightly modified in order to give a better agreement with the quantum results, is

$$\begin{eqnarray}{P}^{(\mathrm{SC})}=\left\{\begin{array}{ll}z/\sqrt{6{D}_{n{\ell }}/{n}^{4}}, & \mathrm{if}\ \ z\leqslant \eta \sqrt{3{D}_{n{\ell }}/2{n}^{4}}\\ 3{j}_{1}{\left(\sqrt{3{D}_{n{\ell }}/2{n}^{4}}/z\right)}^{2}, & \mathrm{if}\ \ z\gt \eta \sqrt{3{D}_{n{\ell }}/2{n}^{4}}\end{array}\right.,\end{eqnarray} \tag{ 10 }$$

where j₁(x) is the spherical Bessel function of order 1, and η = 0.277855 is the solution of the equation ${j}_{1}^{2}(1/x)=x/6$. When used in Equation (6), one gets the SC approximation for the rate coefficient

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{k}}^{(\mathrm{SC})}(n,{\ell },\theta ) & = & \displaystyle \frac{{D}_{n{\ell }}}{4{n}^{4}}\left[\displaystyle \frac{3\sqrt{\pi }}{4{x}^{3/2}}{\rm{erf}}(\eta \sqrt{x})-\displaystyle \frac{3\eta }{2x}\,{e}^{-{\eta }^{2}x}\right.\\ & & \left.+18{\int }_{0}^{\infty }{j}_{1}^{2}(1/y)\,{{ye}}^{-{{xy}}^{2}}\,{dy}\right],\end{array}\end{eqnarray} \tag{ 11 }$$

where $x=3{D}_{n{\ell }}\theta /4{n}^{4}$, and erf is the error function. The integral can be treated as a Mellin–Barnes integral (Paris & Kaminsky 2001) to obtain an asymptotic expansion in the parameter x, because

$$\begin{eqnarray}\begin{array}{rcl}I(x) & = & {\int }_{0}^{\infty }18{j}_{1}^{2}(1/y)\,{{ye}}^{-{{xy}}^{2}}\,{dy}\\ & = & \displaystyle \frac{1}{2\pi i}{\int }_{c-i\infty }^{c+i\infty }F(1-s)G(s)\,{ds},\end{array}\end{eqnarray} \tag{ 12 }$$

where the integral goes along a line parallel with the imaginary axis with 0 < c < 2, and where F and G are the Mellin transforms

$$\begin{eqnarray}\begin{array}{rcl}F(1-s) & = & 18{\int }_{0}^{\infty }{j}_{1}^{2}(1/y)\,{y}^{-s}\,{dy}\\ & = & 18\,\displaystyle \frac{{2}^{3-s}\cos (\pi s/2){\rm{\Gamma }}(s-1)}{(s-6)(s-4)(s-3)}\end{array}\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&G(s)={\int }_{0}^{\infty }{e}^{-{{xy}}^{2}}\,{y}^{s-1}\,{dy}=\displaystyle \frac{{\rm{\Gamma }}(s/2)}{2{x}^{s/2}}.\end{eqnarray} \tag{ 14 }$$

Therefore, the integral is

$$\begin{eqnarray}&&I(x)=18{\int }_{c-i\infty }^{c+i\infty }\displaystyle \frac{{2}^{2-s}\cos (\pi s/2){\rm{\Gamma }}(s-1)}{(s-6)(s-4)(s-3)}\displaystyle \frac{{\rm{\Gamma }}(s/2)}{{x}^{s/2}}\,{ds}.\end{eqnarray} \tag{ 15 }$$

An expansion for small x can be obtained by completing the integration along another line parallel with the imaginary axis, to create a contour around poles on the negative real axis and using the residue theorem. The poles of the integrand occur at negative even numbers s = 0, − 2, − 4, ..., − 2k, ... and are double. The double nature of the poles is the origin of the logarithm in the expansion of I(x). The poles at s = 1, − 1, − 3, ... are removable and do not contribute to the expansion. Therefore, after calculating the residues, one obtains

$$\begin{eqnarray}\begin{array}{rcl}I(x) & = & \displaystyle \sum _{k=0}^{N}{\rm{Res}}(s={s}_{k})+{ \mathcal O }({x}^{N+1})\\ & = & \displaystyle \sum _{k=0}^{N}{A}_{k}({B}_{k}-3\gamma -2\mathrm{ln}2-\mathrm{ln}x){x}^{k}+{ \mathcal O }({x}^{N+1}),\end{array}\end{eqnarray} \tag{ 16 }$$

where γ is the Euler constant, and the coefficients are ${A}_{k}\,=\,18\times {4}^{k}/(k!(2k+1)!(k+2)(k+3)(2k+3))$ and ${B}_{k}\,=1/(k+2)+1/(k+3)+2/(2k+3)$ $+1+1/2+\ldots +1/k+2+2/3$ $+\ldots +2/(2k+1)$. With the first 11 terms listed in Table 1 the expansion of Equation (16) is accurate within one part in 10⁴ for x ≤ 50.

A modified PS64 approximation was obtained in (Guzmán et al. 2017) by replacing the original constant 1/2 for b < R₁ with a linearly increasing transition probability. When written in terms of the dimensionless parameter z, as defined in the discussion following Equation (2), this approximate transition probability is

$$\begin{eqnarray}{P}^{({PS}-M)}=\displaystyle \frac{1}{2}\left\{\begin{array}{ll}z\,{\left(2{P}_{1}\right)}^{3/2}/\sqrt{{D}_{n{\ell }}/{n}^{4}}, & \mathrm{if}\ \ z\leqslant \sqrt{{D}_{n{\ell }}/(2{P}_{1}{n}^{4})}\\ {D}_{n{\ell }}/({n}^{4}{z}^{2}), & \mathrm{if}\ \ z\gt \sqrt{{D}_{n{\ell }}/(2{P}_{1}{n}^{4})}\end{array}\right.,\end{eqnarray} \tag{ 17 }$$

with the free parameter P₁ equal to the transition probability at the matching impact parameter. The integral Equation (6) in this case yields the PS-M approximation

$$\begin{eqnarray}&&\begin{array}{l}{\tilde{k}}^{(\mathrm{PS}-{\rm{M}})}(n,{\ell },\theta )\\ \quad =\,\displaystyle \frac{{D}_{n{\ell }}}{4{n}^{4}}\left[\displaystyle \frac{\sqrt{\pi }}{2{\beta }^{3/2}}\,{\rm{erf}}(\sqrt{\beta })-\displaystyle \frac{1}{\beta }{e}^{-\beta }+{\rm{\Gamma }}(0,\beta )\right],\end{array}\end{eqnarray} \tag{ 18 }$$

with parameter $\beta ={D}_{n{\ell }}\theta /(4{P}_{1}{n}^{4})$.