Efficient computation of collisional $\ell$-mixing rate coefficients in astrophysical plasmas

TL;DR

This paper introduces analytical formulas and a software tool for efficiently calculating collisional $\, ext{l}$-mixing rate coefficients in hydrogen atom collisions, with validation through astrophysical radio recombination line data.

Contribution

It provides new analytical expressions and a software package for rapid, accurate computation of $\, ext{l}$-mixing rate coefficients in astrophysical plasmas, improving upon previous approximation methods.

Findings

Validated rate coefficients against radio recombination lines.

Highlighted the validity range of various approximation methods.

Provided efficient computational tools for astrophysical plasma modeling.

Abstract

We present analytical expressions for direct evaluation of $\ell$-mixing rate coefficients in proton-excited hydrogen atom collisions and describe a software package for efficient numerical evaluation of the collisional rate coefficients. Comparisons between rate coefficients calculated with various levels of approximation are discussed, highlighting their range of validity. These rate coefficients are benchmarked via radio recombination lines for hydrogen, evaluating the corresponding departure coefficients from local thermal equilibrium.

Tables

Table 1. Table 1: Coefficients A k subscript 𝐴 𝑘 A_{k} and B k subscript 𝐵 𝑘 B_{k} for the expansion of the integral in Eq.( 12 ).

k	$A_k$	$B_k$
0	1	7/2 = 3.5
1	1/5 = 0.2	113/20 = 5.65
2	3/350 = 0.00857142857142857	2857/420 = 6.80238095238095
3	2/14175 = 0.000141093474426808	4793/630 = 7.60793650793651
4	1/873180 = 1.14523924047734$\times {10}^{-6}$	28526/3465 = 8.23261183261183
5	1/189189000 = 5.28571957143386$\times {10}^{-9}$	3151273/360360 = 8.74479131979132
6	1/65675610000 = 1.52263526749123$\times {10}^{-11}$	1102667/120120 = 9.17971195471195
7	1/34192364456250 = 2.92462956540933$\times {10}^{-14}$	2439746/255255 = 9.55807329924977
8	1/25408725942600000 = 3.93565581469558$\times {10}^{-17}$	575762023/58198140 = 9.89313443694249
9	1/25932145697017560000 = 3.85621772946852$\times {10}^{-20}$	42376261/4157010 = 10.1939280877361
10	1/35244143470037502000000 = 2.83735083773604$\times {10}^{-23}$	1000753049/95611230 = 10.4668985954893

Equations

$$k(n,\ell ,{\ell }^{\prime },T,R_c)=2\pi {\int }_0^{\infty }{f}_{MB}(v)vdv{\int }_0^{R_c}{P}_{n\ell \to n{\ell }^{\prime }}(b,v)bdb,$$

$$k(n,\ell ,{\ell }^{\prime },T,R_c)=n^4{a}_0^2{v}_0\sqrt{\frac{8\pi \mu v_0^2}{k_{B}T}}{\int }_0^{\infty }z{P}_{n\ell \to n{\ell }^{\prime }}(z)e^{-\theta z^2/2}dz,$$

$$\theta =n^4\frac{\mu v_0^2}{k_{B}T}\frac{a_0^2}{R_c^2},$$

$$\theta =1.704675\times 10^{-10}{n}^4\frac{n_e}{T^2},$$

$$k(n,\ell ,T,n_e)=\sqrt{\frac{\pi }{2}}{a}_0^2{v}_0\sqrt{\frac{\mu v_0^2}{k_{B}T}}{D}_{n\ell }\left[\frac{3\sqrt{\pi }}{4x^{3/2}}\text{erf}(\eta \sqrt{x})-\frac{3\eta }{2x}{e}^{-{\eta }^{2}x}+\sum _{k=0}^{N_t}{A}_k(B_k-3\gamma -\log (4x))x^k\right]$$

$$A_k=\frac{18\times 4^k}{(k+2)(k+3)(2k+3)k!(2k+1)!}$$

$$B_k=\frac{1}{k+2}+\frac{1}{k+3}+\frac{2}{2k+3}+H_k+2H_{2k+1}$$

$$\tilde{k}(n,\ell ,{\ell }^{\prime },\theta )={\int }_0^{\infty }z{P}_{n\ell \to n{\ell }^{\prime }}(z)e^{-\theta z^2/2}dz.$$

$$P^{(B)}=\frac{1}{2}\left\{\begin{array}[]{ll}1&\mbox{, if $z\leq\sqrt{D_{n\ell}/n^{4}}$}\\ D_{n\ell}/(n^{4}z^{2})&\mbox{, if $z>\sqrt{D_{n\ell}/n^{4}}$}\end{array}\right.$$

$${\tilde{k}}^{(B)}(n,\ell ,\theta )=\frac{1-e^{D_{n\ell }\theta /2n^4}}{2\theta }+\frac{D_{n\ell }}{4n^4}\Gamma (0,D_{n\ell }\theta /2n^4)$$

$${\tilde{k}}^{(PS)}(n,\ell ,\theta )=\frac{D_{n\ell }}{4n^4}\left[1-\gamma -\log (D_{n\ell }\theta /2n^4)\right]$$

$$P^{(SC)}=\left\{\begin{array}[]{ll}z/\sqrt{6D_{n\ell}/n^{4}}&\mbox{, if $z\leq\eta\sqrt{3D_{n\ell}/2n^{4}}$}\\ 3j_{1}(\sqrt{3D_{n\ell}/2n^{4}}/z)^{2}&\mbox{, if $z>\eta\sqrt{3D_{n\ell}/2n^{4}}$}\end{array}\right.$$

$${\tilde{k}}^{(SC)}(n,\ell ,\theta )=\frac{D_{n\ell }}{4n^4}\left[\frac{3\sqrt{\pi }}{4x^{3/2}}\text{mbox}erf(\eta \sqrt{x})-\frac{3\eta }{2x}{e}^{-{\eta }^{2}x}+18{\int }_0^{\infty }{j}_1^2(1/y)y{e}^{-x{y}^2}dy\right]$$

$$I(x)=18{\int }_0^{\infty }{j}_1^2(1/y)y{e}^{-x{y}^2}dy=\frac{1}{2\pi i}{\int }_{c-i\infty }^{c+i\infty }F(1-s)G(s)ds$$

$$F(1-s)=18{\int }_0^{\infty }{j}_1^2(1/y)y^{-s}dy=18\frac{2^{3-s}\cos (\pi s/2)\Gamma (s-1)}{(s-6)(s-4)(s-3)}$$

$$G(s)={\int }_0^{\infty }{e}^{-x{y}^2}{y}^{s-1}dy=\frac{\Gamma (s/2)}{2x^{s/2}}$$

$$I(x)=18{\int }_{c-i\infty }^{c+i\infty }\frac{2^{2-s}\cos (\pi s/2)\Gamma (s-1)}{(s-6)(s-4)(s-3)}\frac{\Gamma (s/2)}{x^{s/2}}ds$$

$$I(x)=\sum _{k=0}^N\text{mbox}Res(s=s_k)+\mathcal{O}(x^{N+1})=\sum _{k=0}^N{A}_k(B_k-3\gamma -2\ln 2-\ln x)x^k+\mathcal{O}(x^{N+1})$$

$$P^{(PS-M)}=\frac{1}{2}\left\{\begin{array}[]{ll}z\;(2P_{1})^{3/2}/\sqrt{D_{n\ell}/{n^{4}}}&\mbox{, if $z\leq\sqrt{D_{n\ell}/(2P_{1}n^{4})}$}\\ D_{n\ell}/(n^{4}z^{2})&\mbox{, if $z>\sqrt{D_{n\ell}/(2P_{1}n^{4})}$}\end{array}\right.$$

$${\tilde{k}}^{(PS-M)}(n,\ell ,\theta )=\frac{D_{n\ell }}{4n^4}\left[\frac{\sqrt{\pi }}{2{\beta }^{3/2}}\text{mbox}erf(\sqrt{\beta })-\frac{1}{\beta }{e}^{-\beta }+\Gamma (0,\beta )\right]$$

Full text

numerical methods and codes, Rydberg atoms; abundances-atomic data

††software: Lmixing (Vrinceanu, 2018)

1 Introduction

Energy-conserving angular momentum-changing $n\ell\rightarrow n\ell^{\prime}$ transitions induced in collisions between Rydberg atoms and low velocity ions are needed for accurate comparison between astrophysical observations and models which 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^{4}$, values as large as a few times $10^{5}$ cm3/s are possible for principal quantum numbers near $n\sim 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. 2015, 2017; Salas et al. 2018), and from hydrogen as a tracer of gas ionized by young stars (HII 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 (Izotov & Thuan, 2010; Chluba & Sunyaev, 2006).

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) (PS64 hereafter). Given that the probability for $\Delta\ell=\pm 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 re-examined in Vrinceanu & Flannery (2001a, b) (VF01a and VF01b), where a non-perturbative closed form solution for the transition probability was found.

A semiclassical (SC) rate coefficient for arbitrary $\ell$-changing collisions was derived in (Vrinceanu, Onofrio & Sadeghpour, 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 $|\Delta\ell|>1$. For the dipole allowed transitions ($|\Delta\ell=1|$), probabilities evaluated with VOS12 SC 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 (Guzman 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, Onofrio & Sadeghpour (2017) (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 statistical equilibrium of HRRL.

2 Rate coefficients for angular momentum changing transitions in proton-Rydberg hydrogen scattering

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

[TABLE]

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

[TABLE]

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\alpha)$. The parameter $\theta$, small for large $R_{c}$, is defined as

[TABLE]

where $\mu$ 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

[TABLE]

with $T$ in $K$ and $n_{e}$ in cm*-3*. Depending on the specific physical situations, other choices for $R_{c}$ are possible as discussed in PS64, changing the $\theta$ 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\rightarrow n\ell^{\prime}}$. The relationships among these approximations reflect the organization of the software package (Vrinceanu, 2018), and are illustrated in the diagram in Fig. 1, and discussed in detail below.

The rate coefficients are calculated by numerically integrating Eq. (2). This integral does not pose difficulties for relatively low principal quantum numbers $n\lesssim 40$ when using the VF01b non-perturbative quantum mechanical $P_{n\ell\rightarrow 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 $\sim\ell\times 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 $\ell$, 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.

VOS12 SC formula works well for moderate to large $|\Delta\ell|$, but loses accuracy for larger $|\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 $|\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^{-2}$, as obtained with perturbation theory (PS64).

For dipole transitions, the code provides rate coefficients for combined $\ell\rightarrow\ell+1$ and $\ell\rightarrow\ell-1$ transitions. In the Born approximation, $P_{n\ell\rightarrow n\ell^{\prime}}$ is assumed to be 1/2 for $R<R_{1}$, after which it decreases as $1/b^{2}$. The PS64 formula is derived from this approximation by adopting the additional assumption that $R_{1}<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}<2.98\times 10^{9}$ cm*-3*/K2. 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 (Salgado et al., 2017; Guzman et al., 2017).

The PS-M approximation introduced in Guzman 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_{1}<R_{c}$. The approximate PS-M rate coefficients shown in table I in Guzman 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 Eq. (2), the VOS17 transition probability leads to a SC rate coefficient for dipole allowed transitions

[TABLE]

where erf is the error function, $D_{n\ell}=6n^{2}(n^{2}-\ell^{2}-\ell-1)$, $\eta=0.277855$ is the solution of the equation $j_{1}(1/z)^{2}=z/6$, with $j_{k}(u)$ a spherical Bessel function of order $k$ and argument $u$, and $x=3D_{n\ell}\theta/(4n^{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

[TABLE]

and

[TABLE]

where $H_{k}=\sum_{j=1}^{k}1/j$ for $k>0$, $H_{0}=0$, with $k$ the $k$-th harmonic number and $\gamma=0.57721$ the Euler constant. A derivation of this formula, and a list of the first eleven $A_{k}$ and $B_{k}$ coefficients, is given in the appendix. Our experience shows that a truncation order $N_{t}\sim 16-20$ is sufficient to provide accurate results.

3 Results

As illustrated in Fig. (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 $\ell$, and the temperature in Kelvin. The optional arguments are the electron number density $n_{e}$ in cm*-3*, the final angular momentum $\ell^{\prime}$, and the method of calculation. The choices for this argument are: quantum (default value) using Eq. (2) with the probability defined by Eq. (3) in VF01b, semiclassical using Eq. (5), Born for the Born approximation Eq. (8), PS-M for the approximation in Eq. (18) originally introduced in Guzman et al. (2017), and P_and_S which implements the PS64 method. If $\ell^{\prime}$ is not provided, the rate for the combined $\ell\rightarrow\ell\pm 1$ transitions is calculated. If $\ell^{\prime}$ is given, one can then choose the exact quantum calculation, or the SC approximation (Eq. (8) in VOS12) with the value for the method parameter classical, provided that $|\Delta\ell|>1$.

Figure 2 shows a comparison between the rate coefficients coefficients for dipole allowed transitions for $\ell=1$ and $\ell=n-2$ angular momenta, calculated by integrating the VF01b quantum formula, the PS64 approximation and the SC approximation in Eq. (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 Fig. 2 required 400 digits of precision and took several hours for $\ell=1$ cases and two days for $\ell=n-2$ cases to complete on a single processor. At around $n=200$, in the $\ell=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 Guzman 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\gtrsim 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 Eq. (5) are considered. In the $\ell=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 Eq. (5): on the one hand the prefactor in front of the integral increases as $~{}n^{4}$, while the integral decreases roughly as $1/\theta$. For comparison, we also include in Figure 2 PS-M rate coefficients which behave more consistently with $n$ than the PS64 model in the $\ell=1$ case, and approach the quantum results at large $n$. However, other than the large $n$ limit for $\ell=1$, PS-M predicts rate coefficients which are consistently smaller than QM and SC values, as particularly evident in the $\ell=n-2$ case.

In Figure 3, the accuracy of PS64 approximation (left panel) and of the SC approximation (right panel), both with respect to the quantum rate coefficients, is 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, i.e. 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 Eq. (5). These cases are extreme, but important for the interpretation of HRRLs which 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 rate coefficients approach negative values. This is also expected as for these parameters $\lambda_{D}\leq n^{2}a_{0}$, indicating that the binary collision assumption may fail, with the $\ell$-mixing instead ruled by many-body interactions. If the opposite case of $\ell=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 Fig. 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); Guzman et al. (2016, 2017); 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 Guzman et al. (2017) that for some extreme cases, the PS-M overestimates the quantum results by a factor of 10.

The HRRL results shown in Fig. 4 were produced with the Salgado et al. (2017a) models using the updated $\ell$-changing collision rate coefficients, Eq. (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$.

We show in Fig. 4 the results for the departure coefficients from thermal population, $b_{n}$, for a homogeneous, one dimensional gas slab with: (i) $T_{e}$ = 100 K and $n_{e}$ = 0.1 cm*-3*, and (ii) $T_{e}$ = 10 K and $n_{e}$ = 100 cm*-3*. 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 $\ell$-changing collision rate coefficients primarily affect $b_{n}$ at low to intermediate $n$ values ($n\lesssim 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 $\ell$ sublevels for which radiative decay is less important. Notably, the $\ell$-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}\beta_{n}$, where the correction factor for stimulated emission $\beta_{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\gtrsim 300$. For lower $n$-values the agreement is less good.

We are currently implementing our new $\ell$-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 $\ell$-changing collision rate coefficients than HRRLs (Salgado et al., 2017).

4 Conclusions

We have introduced an efficient and accurate SC approximation for $\ell$-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 multi-precision 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 $\ell$-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 $\ell$ 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.

Acknowledgments

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 (DV) 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.

5 Appendix

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 Fig. 1 we will discuss the four main sections: the integration of the exact quantum formula, the Born approximation, the classical approximation and the semiclassical approximation. These approaches are derived from Eq. (2) either by keeping the exact, but computationally expensive form for the transition probability $P_{n\ell\rightarrow 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 Eq.(2)

[TABLE]

The exact calculation uses the quantum transition probability derived in (VF01b), that 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 pre-computing 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 Eq. (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\rightarrow\ell^{\prime}$ transitions is approximated as

[TABLE]

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

[TABLE]

where $\Gamma(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_{1}$ parameter ($R_{c}>R_{1}$) for any projectile speed. Since the transition impact parameter increases with the speed of the projectile as $\sim 1/v$, the thermal average Eq. (1) will have a contribution from small speeds for which $R_{1}>R_{c}$. This contribution, neglected in PS64, diminishes as $R_{c}\rightarrow\infty$. Indeed, the PS64 rate is obtained from Eq. (8) in the $\theta\rightarrow 0$ limit

[TABLE]

where $\gamma$ is the Euler constant. This corresponds to Eq.(43) in PS64. The rate coefficients in Equation  (9) become negative for $n$ large enough such that $D_{n\ell}\theta/n^{4}>2e^{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}<2.98\times 10^{9}$, which for small $\ell$ reduces roughly to $n^{4}n_{e}/T^{2}<3\times 10^{9}$, and for large $\ell$ to $n^{3}n_{e}/T^{2}<10^{9}$.

A semiclassical 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

[TABLE]

where $j_{1}(x)$ is the spherical Bessel function of order 1, and $\eta=0.277855$ is the solution of the equation $j_{1}^{2}(1/x)=x/6$. When used in Eq. (6), one gets the semiclassical approximation for the rate coefficient

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

Therefore the integral is

[TABLE]

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,\ldots,-2k,\ldots$ 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,\ldots$ are removable and do not contribute to the expansion. Therefore, after calculating the residues, one obtains

[TABLE]

where $\gamma$ 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 Eq. (16) is accurate within one part in $10^{4}$ for $x\leq 50$.

A modified PS64 approximation was obtained in (Guzman et al., 2017) by replacing the original constant 1/2 for $b<R_{1}$ with a linearly increasing transition probability. When written in terms of the dimensionless parameter $z$, as defined in the discussion following Eq. (2), this approximate transition probability is

[TABLE]

with the free parameter $P_{1}$ equal to the transition probability at the matching impact parameter. The integral Eq. (6) in this case yields the PS-M approximation

[TABLE]

with parameter $\beta=D_{n\ell}\theta/(4P_{1}n^{4})$.