Formation of NaCl by radiative association in interstellar environments

Context. Radiative association is a possible way of sodium chloride (NaCl) formation in interstellar and related environments. Theo- retical studies are essential since laboratory experiments are unavailable and difficult to perform. Aims. The total rate coefficient was calculated for the formation of NaCl by radiative association at 30–750 K. Methods. We included two contributing processes for the total rate-coefficient computation. One of them takes the nonadiabatic cou- pling between the two lowest 1 Σ + states, X 1 Σ + and B 1 Σ + , into account. The other one was calculated conventionally as a single channel and started in the continuum of the A 1 Π state. The individual rate coefficients were calculated from cross sections obtained up to 0.8 eV, which enabled us to calculate the rate coefficients up to 750 K. The cross section was also calculated for a one-state process within the X 1 Σ + state. Results. The nonadiabatic coupling enhances the formation of NaCl by radiative association by two orders of magnitude at about 30 K and by around one order of magnitude at about 750 K. The single-channel process starting in the continuum of the A 1 Π state starts to contribute above around 200 K. The one-state transition model, within the X 1 Σ + state, is not an adequate approximation for collisions in 1 Σ + symmetry. Instead, these collisions are treated in the diabatic representation in the total rate-coefficient calculation. Conclusions. The calculated total rate-coefficient function at 30–750 K can improve the astrochemical reaction networks for the CRL 2688, IRC+10216, and Orion SrcI environments, where NaCl was detected before.


Introduction
Sodium chloride (NaCl) is the first molecule detected as a bearer of a sodium atom in space (Cernicharo & Guélin 1987;Cernicharo et al. 1988). The gas-phase NaCl (together with potassium chloride) has been observed in different space environments, mostly in inner layers of expanding circumstellar envelopes of stars in their final evolution stage (asymptotic giant branch or post asymptotic giant branch stars). It was detected in shells of carbon-rich stars IRC+10216 (Cernicharo & Guélin 1987;Agúndez et al. 2012) and CRL 2688 (Highberger et al. 2003) with abundances relative to H 2 on the order of 10 −10 -10 −8 . Similar fractional abundances of NaCl were observed in the nebulae surrounding oxygen-rich stars VY CMa ), IK Tauri , and QX Pup (Sánchez Contreras et al. 2018). It is assumed that in these environments, the NaCl molecules survive in the gas phase only for a relatively short time before being condensed onto dust grains. Recently, the first observation of the gas-phase NaCl, not related to evolved stars, was made above the protostellar oxygen-rich disk around the massive accreting young star Orion SrcI (Ginsburg et al. 2019).
One of the possible formation mechanisms for NaCl in the interstellar medium and related environments is radiative association (RA). Temperature-dependent RA rate coefficients belong to the input data for various interstellar-environment models and, because they are very difficult to measure experimentally, their accurate calculations are of increased importance. In the recent study Gustafsson (2020), the rate coefficients for radiative association of sodium with chloride in the temperature range 1-30 K were calculated. In that fully quantum mechanical approach, the nonadiabatic coupling between the ionic and atomic electronic states was accounted for, thus resulting in two-channel coupled calculations. Compared to the conventional single-channel calculations, the differences in the calculated cross section (see Gustafsson 2020, Fig. 4), and consequently also in the calculated rate coefficient, were in orders of magnitude. This clearly demonstrates the importance of including the effects of nonadiabatic dynamics in systems with a presence of avoided crossings between the Born-Oppenheimer (BO) potential energy curves.
The origin of significant nonadiabatic effects in NaCl comes from its intricate electronic structure. In the diabatic picture, the ground electronic state of the NaCl molecule is the ionic molecular state Na + Cl − , which dissociates to Na + ( 1 S)+Cl − ( 1 P). On the other hand, the energetically lowest NaCl dissociation channel corresponds to Na( 2 S)+Cl( 2 P). Employing the molecular symmetry group C ∞v , the ionic molecular state belongs to 1 Σ + symmetry, and four molecular states of the respective 1 Σ + , 1 Π, 3 Σ + , and 3 Π symmetries dissociate via the lowest dissociation channel. Consequently, in the adiabatic picture, the BO ground-state X 1 Σ + potential energy curve exhibits an avoided crossing with the BO potential energy curve of the first excited 1 Σ + electronic state, B 1 Σ + . On the BO ground-state potential energy curve, thus, the ionic configuration, which is present in the potential well, changes into the neutral-atom configuration just at this avoided crossing. Within the adiabatic picture, such an exchange of configurations is reflected by a sharp peak in the nonadiabatic coupling element between the ground and excited states at the avoided-crossing position. Within the diabatic picture at the avoided-crossing position, the ionic and atomic states simply cross, and the corresponding nondiagonal potential matrix element is nonzero.
Our fully quantum mechanical study follows that of Gustafsson (2020) and concentrates on the conditions present in some environments, where the NaCl molecules have been observed. For this purpose, it significantly extends the temperature range to 30-750 K, which should cover local temperatures in the CRL 2688 shocks (Highberger et al. 2003), the IRC+10216 expanding shell (Agúndez et al. 2012), and the accretion disk of Orion SrcI (Ginsburg et al. 2019).
Coupled two-channel calculations were performed in order to determine the energy-dependent cross sections including the abovementioned nonadiabatic dynamic effects between the ground and first excited 1 Σ + electronic states (hereby denoted 1 Σ + → 1 Σ + ). Uncoupled two-state and one-state calculations were used to determine the cross sections for processes A 1 Π → X 1 Σ + and X 1 Σ + → X 1 Σ + , respectively. From the former two cross sections, the corresponding rate coefficients were calculated as functions of temperature, and the sum of the 1 Σ + → 1 Σ + and A 1 Π → X 1 Σ + rate coefficients provides our estimate of the total rate coefficient.

Methods
A detailed description of Sects. 2.1-2.3 can be found in Gustafsson (2020).

General theory
The aim is to obtain the total rate coefficient, α tot (T ), for the formation of NaCl via collisions of Na( 2 S) and Cl( 2 P), which is a sum over all rate coefficients for individual transitions in this colliding channel. In Eq. (1), T is temperature, µ is the reduced mass of 23 Na 35 Cl, k B is the Boltzmann constant, and is the RA cross section at collision energy E i (Babb & Dalgarno 1995). Here, ϵ 0 denotes the vacuum permittivity, k 2 = 2µE i /ℏ 2 , ω is the angular frequency of an emitted photon, S J i →J f are the Hönl-London factors, M J i ;v f ,J f (E i ) is the dipole matrix element associated with the dipole moment of the transition, and p i is the statistical weight (the probability of collision in the initial electronic state) expressed as where S and Λ denote the electronic spin quantum number and projection of the orbital electronic angular momentum of the initial molecular state, respectively. For the colliding Na and Cl atoms, S Na and S Cl denote the spin quantum numbers and L Na and L Cl denote the orbital angular momenta.
Nonadiabatic dynamics. The theory was mainly inspired by Tchang-Brillet et al. (1992), and we only summarize it here. In the nonadiabatic picture, we need to solve the coupled Schrödinger equation for the initial state. In Eq. (4), 1 is a two-dimensional unit matrix, where U 1 (R) and U 2 (R) are potential energies of the atomic and ionic electronic states in the diabatic representation, and the coupling between them is denoted U 12 (R) = U 21 (R). We note that Ψ(R, E i , J i ) is a column vector with components ψ 1 (R, E i , J i ) and ψ 2 (R, E i , J i ), and it denotes the radial initial wave function with collision energy E i and the rotational quantum number J i . The component functions ψ 1 (R, E i , J i ) and ψ 2 (R, E i , J i ) are associated with the electronic states U 1 (R) and U 2 (R), respectively, and together describe the initial state. Then, is the dipole matrix element associated with the dipole moment D 2 (R) of the diabatic closed-channel component ψ 2 (R, E i , J i ) of the continuum state, and the radial wave function of the final bound state ψ v f ,J f (R), associated with U 2 (R). The Schrödinger equation, Eq. (4), was solved with the asymptotic boundary condition and with the condition that ψ 2 (R → ∞, E i , J i ) → 0.
Born-Oppenheimer approximation. Conventionally (Zygelman & Dalgarno 1990;Antipov et al. 2009;Zámečníková et al. 2019;Bai et al. 2021), couplings of electronic states in a molecular system are neglected. This negligence is often appropriate or justifiable, but in some cases it is not. For Na( 2 S) + Cl( 2 P) approaching in the 1 Σ + molecular electronic state, it has not been neglected because of the reasons described in Sect. 1.

Ab initio data
All ab initio data were taken from Giese & York (2004) and Zeiri & Balint-Kurti (1983) and are illustrated in Fig. 1. The diabatic potential energy curves U 1 (R) and U 2 (R) are shown in the top-left panel together with the diabatic coupling U 12 (R). The vicinity of 18.37 a 0 is shown as an inset in the top-left panel. The bottom-left panel shows the dipole-moment function D 2 (R) of the ionic state in the diabatic representation. The potential energy curves for the X 1 Σ + and A 1 Π electronic states in the BO approximation are illustrated in the top-right panel. In the bottom-right panel, the dipole-moment function in the X 1 Σ + state and the transition dipole-moment function for the A 1 Π → X 1 Σ + transitions are illustrated. The adiabatic states within the BO approximation, X 1 Σ + , B 1 Σ + , and A 1 Π, are illustrated in Top left: the potential energy curves for diabatic states U 1 (R) and U 2 (R) and their coupling U 12 (R) are illustrated. Bottom left: the dipole moment of the diabatic state 2 is illustrated. Top right: the potential energy curves of the X 1 Σ + and A 1 Π states within the BO approximation. Bottom right: the X 1 Σ + permanent dipole-moment and A 1 Π → X 1 Σ + transition dipole-moment functions. All the data were taken from Giese & York (2004) and Zeiri & Balint-Kurti (1983
The continuum radial wave functions were calculated by the renormalized Numerov method (Numerov 1923). The discrete variable representation was used for the bound-state radial wave function calculation (Colbert & Miller 1992) with the same input parameters as Gustafsson (2020). The trapezoidal method was used for the integration over collision energy in Eq. (1) to obtain the rate coefficients.

Cross sections
The cross sections for formation of NaCl via radiative association were calculated in the collisions of Na( 2 S) and Cl( 2 P). Three different processes were studied: (i) the nonadiabatic dynamics for the 1 Σ + → 1 Σ + process, (ii) the X 1 Σ + → X 1 Σ + process, and (iii) the A 1 Π → X 1 Σ + process. Here, we calculated the cross sections for 0.03 eV ≤ E i ≤ 0.8 eV.
3.1.1. The 1 Σ + → 1 Σ + process The cross section for the 1 Σ + → 1 Σ + process was obtained by including the nonadiabatic coupling as described in Sect. 2.1. The calculations were performed in the same manner as by Gustafsson (2020), but the maximum of the internuclear distance R 2 had to be modified with increasing collision energy E i . While Gustafsson (2020) used 20 a 0 , here it varied from 21 a 0 for 0.03 eV up to 41 a 0 for 0.8 eV. The radial spacing dR stayed the same as in Gustafsson (2020), equal to 0.002 a 0 . The search of appropriate R 2 and confirmation of dR are described in Appendix B.2.
In addition, J i had to be increased with larger E i . The highest rotational quantum number for the final ro-vibrational bound state of NaCl is 489, which was estimated with the discrete variable routine. From S J i →J i ±1 , the maximum initial rotational quantum number J max Cross sections for formation of NaCl by radiative association. Top: the cross section for the 1 Σ + → 1 Σ + process is depicted from 0.03 eV to 0.8 eV with dE i = 10 −5 eV. It takes the nonadiabatic coupling into account. Bottom: the cross sections are depicted from 10 −5 eV to 0.8 eV for three different processes: 1 Σ + → 1 Σ + , A 1 Π → X 1 Σ + , and X 1 Σ + → X 1 Σ + . The green lines represent cross sections up to 0.03 eV taken from Gustafsson (2020). The purple lines represent cross sections calculated here. NAD denotes nonadiabatic dynamics. MG (2020) denotes Gustafsson (2020). A logarithmic scale is used on the x-and y-axes.
Finally, the energy grid was also changed from the set of Gustafsson (2020) because the density of resonances increased with increasing collision energy. The energy grid step was decreased to dE i = 10 −5 eV. The sufficiency of this choice was checked (see Appendix B.3).
The energy-dependent cross section is illustrated in the top panel of Fig. 2. The cross-section function increases with increasing collision energy up to around 0.7 eV. The background contribution then slowly decreases. The cross section has a rich resonance structure with an increasing number of resonances with increasing collision energy. Gustafsson (2020) showed that for low collision energies, one can see the resonances in sort of "beats." These beats become more frequent at larger collision energies. Here, one can vaguely see two or three beats above 0.03 eV because they start to overlap. This structure cannot be seen from around 0.04 eV. In the bottom panel of Fig. 2, the whole cross-section function from 10 −5 eV is shown, where one can see the beats for lower E i .
The low-energy-lying resonances in the bottom panel of Fig. 2 resemble shape resonances. In the BO approximation, shape resonances would stem from the continuum of the X 1 Σ + state. It seems that below 1 meV, the nonadiabatic effects are less important.
The high-energy resonances in the cross section stem from the diabatic potential energy curve U 2 (R), which correlates to the dissociation limit of Na + ( 1 S) + Cl − ( 1 S). The leading term of the effective potential in the long range is the Coulomb potential, which dies out slowly. This enables many Feshbach resonances, related to the continuum component ψ 2 (R, E i , J i ), to exist. The Feshbach resonances are known for showing dips due to the destructive interference. This can be seen clearly in the bottom panel for lower energies, but also in Figs. B.1 and B.2 for higher collision energies.
3.1.2. The X 1 Σ + → X 1 Σ + process The cross section for the X 1 Σ + → X 1 Σ + process was obtained by the conventional approach within the BO approximation. Similarly, as in the previous Sect., they were calculated from 0.03 eV to 0.8 eV with the grid steps 10 −5 eV and 10 −4 eV up to 0.15 eV and from 0.15 eV to 0.8 eV, respectively. The input parameters from Gustafsson (2020) differed slightly: the maximum internuclear distance for the bound-state and continuum wave functions was increased to 21.1 a 0 , and the grid step was increased to 0.024 5 a 0 and 0.002 45 a 0 for the bound-state and continuum-state wave functions, respectively. The maximum initial rotational quantum number was increased to 479 up to 0.06 eV, and 489 up to 0.8 eV.
The cross section is illustrated in the bottom panel of Fig. 2. The values are expectedly low because of the one-electronicstate nature. From about 10 −2 eV, the background contribution sort of "curls." This is known for RA cross sections at larger collision energies, and it has been observed, for example, by Zygelman & Dalgarno (1990) and Augustovičová et al. (2015). The cross-section function decreases with increasing collision energy from about 0.001 eV. Regarding the resonance contribution, many resonances can be seen in this single channel.
3.1.3. The A 1 Π → X 1 Σ + process The cross section for the A 1 Π → X 1 Σ + process was obtained by the conventional approach within the BO approximation. It was calculated for collision energies 0.03-0.8 eV with dE i = 0.001 eV up to 0.1 eV and dE i = 0.01 eV up to 0.8 eV. The maximum internuclear distance was set to 14 a 0 and 10 a 0 up to 0.1 eV and above 0.1 eV, respectively. The grid step was dR = 0.002 a 0 and 0.02 a 0 for continuum-and bound-state wave functions, respectively. The maximum initial rotational quantum number was increased with the increasing collision energy, concretely up to 0.1 eV J max i = 100 and above 0.1 eV J max i = 330. The cross section is illustrated in the bottom panel of Fig. 2. The values increase with increasing collision energy. At 0.03 eV, the values are about two orders of magnitude smaller than those from the nonadiabatic dynamics for the 1 Σ + → 1 Σ + process. At 0.8 eV, the values are less than one order of magnitude smaller than those from the 1 Σ + → 1 Σ + process. The A 1 Π state does not support any resonances because of its repulsive nature.

Rate coefficient
The rate coefficients were calculated from the cross sections for two studied processes: the 1 Σ + → 1 Σ + process treated nonadiabatically and the A 1 Π → X 1 Σ + process obtained conventionally. The corresponding rate coefficients for these two processes were , MG (2020) A 1 Π →X 1 Σ + , this work A 1 Π →X 1 Σ + , MG (2020) Fig. 3. Total rate coefficient for radiative association of Na( 2 S) with Cl( 2 P) obtained as a sum of the rate coefficients of the 1 Σ + → 1 Σ + and A 1 Π → X 1 Σ + processes up to 750 K. The results of individual processes are compared to the work of Gustafsson (2020). The solid red line is on top of the black and dashed red lines. The solid blue line is on top of the dashed blue line. A logarithmic scale is used on the x-and y-axes. obtained from the cross sections calculated here together with those from Gustafsson (2020). The presented rate-coefficient functions were calculated from 10 K to 750 K. The 1 Σ + → 1 Σ + process was preferred over the X 1 Σ + → X 1 Σ + process since the nonadiabatic coupling cannot be neglected in the collisions in 1 Σ + symmetry. This way, we obtained the total rate coefficients for the formation of NaCl in collisions of Na( 2 S) and Cl( 2 P) via RA.
The total rate coefficient is illustrated in Fig. 3 together with its summands. Between 10 K and 30 K, the rate coefficients for the individual processes are compared to the rate coefficients obtained by Gustafsson (2020). The rate coefficients up to 30 K match those from Gustafsson (2020). The rate-coefficient function for the 1 Σ + → 1 Σ + process increases almost linearly in a log-log plot with increasing temperature. While the A 1 Π → X 1 Σ + rate coefficient increases as well, its values are considerably smaller than those for the 1 Σ + → 1 Σ + process. Therefore, the formation of NaCl by RA mainly takes place via nonadiabatic coupling in the 1 Σ + symmetry which should not be neglected in collisions of Na and Cl in their ground electronic states.
The total rate coefficient up to 30 K did not change if the cross sections calculated here (from 0.03 eV) were included (not shown). The A 1 Π → X 1 Σ + process starts to contribute visibly to the total rate coefficient above about 200 K.
Several computations to confirm the convergence of the rate coefficient were performed. Apart from those computations related to the convergence of the cross section, we checked the sufficiency of the energy grid for the integration in Eq. (1). These results are shown in Appendix C.
A rate-coefficient function can be represented by the Kooij function where α, β, and γ are fitting parameters. The fitting parameters for the present total rate coefficient are summarized in Table 1 for three different temperature ranges. The maximum relative error was kept below 0.001. Notes. x(−y) ≡ x × 10 −y .

Conclusion
In this work, we have calculated the total rate coefficient for the formation of NaCl through RA of Na and Cl in their electronic ground states. Temperatures up to 750 K are considered, which makes the data relevant for a number of known astrophysical objects, for example IRC+10216, CRL 2688, and Orion SrcI (Agúndez et al. 2012;Highberger et al. 2003;Ginsburg et al. 2019). The method which accounts for nonadiabatic dynamics was developed by Gustafsson (2020) in a prior work, where temperatures up to 30 K were considered.
Analysis of the cross sections shows that RA for an approach in the 1 Σ + state dominates A 1 Π completely from about 0.5 meV. The total rate coefficient calculated in this work is also dominated by the 1 Σ + → 1 Σ + reaction. An approach on neither of the triplet states is expected to contribute in any significant way since: (1) the corresponding potential energy curves of final triplet states are repulsive, and (2) molecule formation mediated by spin-orbit coupling would, thus, be required. This effect, however, has been shown to be weak (Antipov et al. 2011).
The cross section study also makes it clear that a single channel formation model on the adiabatic X 1 Σ + state is increasingly, and completely, insufficient with higher collision energy. The X 1 Σ + → X 1 Σ + calculation underestimates the cross section by more than one and seven orders of magnitude at 0.1 meV and 0.8 eV, respectively.
The magnitude of the rate coefficient, which is about 10 −16 cm 3 s −1 at 500 K and still increasing with temperature, is higher than for most other RA reactions with atoms or atomic ions in the ground state, for example CO, CN, and CO + (Nyman et al. 2015;Zámečníková et al. 2020). One exception arises in environments where quasibound states of the collision complex are populated by alternative mechanisms, such as by radiation or three-body interactions. Then the local thermal equilibrium (LTE) rate coefficient should be applied (Forrey et al. 2016), instead of that from the conventional non-LTE, zero density limit (NLTE-ZDL), which is typically obtained by computing the resonance contribution with the conventional Breit-Wigner formula by accounting for the radiative broadening of the resonant states (Gustafsson & Forrey 2019;Gustafsson et al. 2012). The rate coefficient in this work was calculated by straightforward integration of the cross section obtained with the perturbation theory formula, Eq. (2), without accounting for the radiative width of the resonant states. However, the resonances that manifested in the 1 Σ + → 1 Σ + cross section appear to be rather short-lived since they are not as tall as for many other systems (e.g., CO, CO + , and CN). Thus, we expect LTE and NLTE-ZDL to give essentially the same rate coefficient for the formation of NaCl studied in this work. As a consequence, the rate coefficient from the present work can be used in environments independent of the level of background radiation and occurrence of three-body interactions.
A5, page 5 of 9 A&A 664, A5 (2022) In future studies, it would be desirable to extend the calculations of cross sections for NaCl formation to energies surpassing the ionic dissociation limit (at about 1.48 eV). At that energy, the channel for chemi-ionization Na( 2 S) + Cl( 2 P) → Na + ( 1 S) + Cl − ( 1 S) would open and compete with RA. An estimate of the corresponding cross section has been obtained for a few energy points by Cooper et al. (1987). Chemi-ionization is implicitly accounted for in the implementation of the coupled Schrödinger equation, Eq. (4). The coulombic boundary condition is chosen for the ionic channel as it opens, and the chemi-ionization and RA cross sections can be obtained simultaneously. It would be interesting and useful to pursue such a calculation in order to obtain a complete picture of the outcomes from collisions between Na and Cl in their ground states. Potential energy curves for the X 1 Σ + , B 1 Σ + , and A 1 Π states of NaCl are illustrated. They were taken from Giese & York (2004) and Zeiri & Balint-Kurti (1983).

Appendix A: Potential energy curves within the Born-Oppenheimer approximation
The potential energy curves of the three lowest electronic singlet states of NaCl are illustrated in Fig. A.1. The potential energy curves were obtained by Giese & York (2004) and Zeiri & Balint-Kurti (1983) within the BO approximation. The X 1 Σ + and B 1 Σ + states undergo an avoided crossing around 18.37 a 0 , the effect of which on RA is studied in this work.

(NAD) cross section
The convergence for the cross-section calculations was studied in: i) the initial rotational quantum number, J i ; ii) the internuclear distance for the radial continuum wave function obtained by the renormalized Numerov method, R 2 ; iii) the step in the Numerov method, dR; and iv) finally in the energy step, dE i , to cover the resonance contribution. In order to achieve the convergence in all aspects, a series of calculations were performed.

B.1. Convergence in J i
Firstly, the convergence in J i was studied. As mentioned in Sect. 3.1, the maximum initial rotational quantum number J max i = 490. However, the convergence in J i was found with lower J i for collision energies lower than 0.5 eV. It is common that the transitions to higher rotational quantum numbers at lower collision energies are negligible. Their significance rises with increasing collision energies.
In order to estimate J max i , we calculated the partial cross sections and summed them from J i = 0 to J i = 20, then from J i = 21 to J i = 40, from J i = 41 to J i = 60, etc. These partially summed partial cross sections were calculated from 0.03 eV to 0.1 eV with the energy step dE i = 10 −4 eV and from 0.1 eV to 0.8 eV with dE i = 10 −3 eV. The lowest J max

B.2. Convergence in R 2 and dR
With larger collision energies, it was needed to increase the maximum internuclear distance, R 2 , up to which the radial continuum wave function of the initial state was calculated. Therefore, we calculated the partial cross sections for J i = 0, 20, 40, 60, . . . , 300 at collision energies from 0.1-0.2 eV and for J i = 0, 20, 40, 60, . . . , 460, 480, 490 at collision energies from 0.5-0.6 eV and from 0.7-0.8 eV with various R 2 . This enabled us to estimate sufficient R 2 for different partial waves at different collision energies. Fig. B.1 illustrates this convergence search. For 0.1 eV -0.2 eV, we show the partial wave J i = 80; for 0.5 eV -0.6 eV, the partial wave J i = 130; and for 0.7 eV -0.8 eV, the partial wave J i = 160. One can see that the partial cross section with J i = 80 between 0.1 eV and 0.2 eV is converged with respect to the choice of R 2 , when the radial continuum wave function is calculated up to 22 a 0 . In the middle panel, one can see that the partial cross section with J i = 130 between 0.5 eV and 0.6 eV converges if R 2 is equal to at least 31 a 0 . Finally, the partial cross section with J i = 160 converges between 0.7 eV and 0.8 eV when R 2 = 40 a 0 . These partial waves belong to the waves that contribute the most at these shown energies. Two things can be observed: 1) For larger E i , larger R 2 was needed. 2) At larger collision energies, the partial waves with larger J i contribute more.
We do not show in Fig. B.1 how we estimated R 2 for 0.1 eV < E i , 0.2 eV < E i < 0.5 eV, and 0.6 eV < E i < 0.7 eV. The idea is that when a partial cross section is converged up to a certain collision energy, we assumed that up to that energy, one can use the same R 2 . We also assumed that for the partial wave between two certain J i , the convergence is achieved with the same R 2 .
When the convergence in J i and R 2 was achieved, we studied the influence of different dR in the renormalized Numerov method for the calculations of radial continuum wave functions. Figure B.2 shows the partial cross sections calculated with different spatial grids: dR = 0.01 a 0 , 0.005 a 0 , 0.002 a 0 , and 0.001 a 0 . One can see that with the larger step, 0.01 a 0 , the resonances' positions are not converged yet. All three smaller steps allowed us to find the resonances' positions at around the same energy; however, their heights differ. The difference is largest for 0.005 a 0 . Some heights are the same with both steps: 0.002 a 0 and 0.001 a 0 . The results confirm that the same step, dR = 0.002 a 0 , as used by Gustafsson (2020) suffices here as well.

B.3. Resonance contribution
Finally, the cross sections were calculated from about 0.03 eV up to 0.8 eV with the energy step 10 −5 eV. From 0.1 eV, we checked after every 0.1 eV whether a larger energy step could suffice for larger collision energies. However, by deleting every other grid point, we would lose some resonances. Therefore, the energy step 10 −5 eV was used in all our calculations presented here in this paper. Fig. B.3 illustrates this procedure and shows that dE i = 2 × 10 −5 eV is not sufficient to include all important resonances at energies 0.789 eV < E i < 0.8 eV.
By the choice of the energy step to 10 −5 eV, we eliminated resonances with the width narrower than this energy step. With a smaller energy step, some of these resonances could have been involved in the rate-coefficient calculation, but there would always be some narrow resonances, narrower than a chosen energy step, that could not be included, and the narrower a resonance is, the less it contributes to the rate coefficient. Assuming that there are not that many of them and that they do not have A5, page 7 of 9 A&A 664, A5 (2022) Convergence search of partial cross sections with respect to the maximum internuclear distance, R 2 , for radial continuum wave functions is illustrated at 0.1 eV to 0.2 eV for J i = 80, at 0.5 eV to 0.6 eV for J i = 130, and at 0.7 eV to 0.8 eV for J i = 160 with the energy spacing dE i = 10 −4 eV. From the bottom to the top panel, the collision energies increase. Different maximum internuclear distances were used for the radial continuum wave functions. The red line represents a nonconverged partial cross section. The black line represents the converged partial cross section. The dashed green line shows that the convergences were achieved at the black line. The black line is on top of the red line, where they coincide, and the dashed green line is on top of the black line. A logarithmic scale is used on the y-axis. J i =160, R 2 =40 a 0 , dR=0.01 a 0 J i =160, R 2 =40 a 0 , dR=0.005 a 0 J i =160, R 2 =40 a 0 , dR=0.002 a 0 J i =160, R 2 =40 a 0 , dR=0.001 a 0 Fig. B.2. Convergence search of partial cross sections with respect to grid density for radial continuum wave functions is illustrated at 0.1 eV to 0.2 eV for J i = 80, at 0.5 eV to 0.6 eV for J i = 130, and at 0.7 eV to 0.8 eV for J i = 160 with the energy spacing dE i = 10 −4 eV. The red and dashed green lines represent partial cross sections calculated with dR = 0.01 a 0 and dR = 0.005 a 0 , respectively. The dashed black line represents a partial cross section calculated with dR = 0.002 a 0 , which was chosen for the total cross-section calculation. The dashed brown line represents a partial cross section calculated with a denser grid, dR = 0.001 a 0 . Where they coincide, the dashed green line is on extreme heights, we neglected them for the purpose of calculating the rate-coefficient function.
There is one more approach how one can treat the resonance contribution and that is to treat the background and resonance contribution separately. One can use the Breit-Wigner theory to obtain the resonance contribution if the positions, rotational quantum numbers, and widths are known. For example, readers can refer to Zámečníková et al. (2019) for more details.
One can be even more prudent and divide this resonance contribution into narrow and wide categories (Augustovičová et al. 2012;Zámečníková et al. 2018). The wide resonances can be treated similarly as the background contribution, but with the energy step related to the tunneling width. The Breit-Wigner theory can then merely be used for the narrow resonances.
The abovementioned separate treatments are, however, well elaborated within the BO approximation where the nature of resonances is purely additive. Feshbach resonances generally show constructive and destructive interference and are not straightforward to parametrize. Therefore, we decided to choose a grid dense enough to obtain the "best" cross section and corresponding rate coefficient.

B.4. Note
We note that the background contribution in this work was already roughly converged with the parameters used by Gustafsson (2020), who obtained the cross section up to 30 meV, concretely R 2 = 20 a 0 , dR = 0.002 a 0 , and with a varying dE i up to 2 × 10 −5 eV. It is important to note that J max i was expected to change with increasing collision energy. We could see this when we calculated the cross section with the same parameters with dE i = 10 −3 eV (not shown here).