Formation of the BeH+ and BeD+ Molecules in Be+ + H/D Collisions Through Radiative Association

Cross sections and rate coefficients for the formation of BeH+ and BeD+ molecules in Be+ + H/D collisions through radiative association are calculated using quantum mechanical perturbation theory and Breit-Wigner theory. The local thermodynamic equilibrium limit of the molecule formation is also studied, since the process is also relevant in environments with high-density and/or strong radiation fields. The obtained rate coefficients may facilitate the kinetic modelling of BeH+/BeD+ production in astrochemical environments as well as the corrosion chemistry of thermonuclear fusion reactors.


INTRODUCTION
Lithium, beryllium and boron are surprisingly rare in the Universe compared to other nuclides of similar mass, having abundances more similar to heavy nuclei in meteorites. The widely accepted explanation for this observation is that these elements cannot be produced in normal stellar nucleosynthesis while also being depleted in hydrogen burning processes, producing helium (Burbidge et al., 1957). These elements can still be formed in spallation reactions of the more abundant elements e.g., carbon, nitrogen and oxygen. Out of these three elements, beryllium is the rarest, with a solar abundance of logϵ Be 1.15 ± 0.20 (with logϵ H 12.00 by definition) (Chmielewski and Müller, 1975).
Even though (or exactly because) lithium, beryllium, and boron have low abundance, these elements can play major roles in observational cosmology. The primordial lithium problem is a famous unsolved problem in modern physics, which states that the current nucleosynthesis model correctly predicts the abundances of deuterium and helium but fails to account for the observed amount of lithium (Fields, 2011). Beryllium could still be formed in extremely metal-deficient stars, and for this reason, the use of beryllium concentrations as a "big bang calorimeter" has been proposed (Pospelov and Pradler, 2011). The presence of BeH molecule in the Sun was first proposed by Wohl (1971). Later Shanmugavel et al. (2008) gave an evidence for the presence of BeH and its isotopologues by observing their A-X transitions in sunspot umbral spectra.
Another important aspect of the beryllium chemistry is that beryllium is one of the main candidates of coating material for fusion reactors (Federici et al., 2001;Paméla et al., 2007). The main issue with the beryllium coating is that Be erodes easily under plasma exposure by releasing Be and Be + into the plasma, where molecules and molecular ions will be formed. In principle, Be, Be + , Be 2+ , BeX, BeX + , BeX 2 and BeX + 2 species can all be important with X being H, D or T. So in order to aid fusion reactor design, the formation and the reactions of these molecules must be understood.
Because of their importance, the BeH and the BeH + molecules have been extensively studied through theoretical calculations in the last decades. Machado and Ornellas (1991) have shown the possible pathways of the BeH + formation from Be + and H. For the electronic states Be + ( 2 S) + H ( 2 S) there are two possible correlations to BeH + , the X 1 Σ + being the lowest energy attractive state, and the a 3 Σ + being a higher energy repulsive state, not leading to a product. The Be + ( 2 P) state has about 4.0 eV excess energy compared to the Be + ( 2 S) state, and the process Be + ( 2 P) + H ( 2 S) has three reactive channels (A 1 Σ + ,b 3 Π and B 1 Π in increasing energy order) and one nonreactive channel c 3 Σ + . Recent ab initio calculations of the BeH + molecule can be found in Refs. (Koput, 2013;Wan et al., 2014), and for BeH a recent review is available (Dattani, 2015). Collisions with electrons can be important in fusion reactor applications, and the BeH + + e − process has been investigated for example in Refs. (Chakrabarti and Tennyson, 2012;Niyonzima et al., 2018;Iacob et al., 2020;Pop et al., 2020;Pop et al., 2021).
The molecule formation in the collision of the Be + and H/D particles is forbidden in a strict sense as the relative translation energy of the reactant atoms cannot be stored in the newly formed vibrational degree of freedom of the diatomic molecule without dissociating back to the atomic components. The product can only be stabilized either by the environment (by thermal collisions with free molecules or by interaction with a solid surface) or by emitting a photon. The last process is called radiative association.
The aim of this present work is to investigate the BeH + and BeD + production by radiative association in absence of spinorbit and nonadiabatic couplings using quantum dynamical methods.

Theory
The cross section of radiative association (RA) can be obtained by quantum mechanical perturbation theory (Babb and Dalgarno, 1995) which yields the Golden rule-like formula where the sum runs over the initial angular momentum (J) and final rotational (j′) and vibrational (v′) quantum numbers. k ini 2 μE c /Z is the wavenumber, where μ is the reduced mass of the colliding fragments, S J,j′ is the so-called Hönl-London factor, λ J,j′,v′,Ec is the wavelength of the emitted photon, E c is the collision energy, and f stat is the statistical weight factor determined by the initial atomic and final molecular electronic states. The transition dipole matrix elements are defined as where D(r) is the dipole moment function, χ ini J,E c (r) is the radial part of the energy normalized scattering wavefunction of the initial state and φ fin j′,v′ (r) is the radial part of the final rovibrational wavefunction, normalized to unity.
The thermal rate constant can be calculated from the cross section according to collision theory as where T is the temperature and k B is the Boltzmann constant. Besides quantum mechanical perturbation theory, we use Breit-Wigner theory to calculate the RA cross sections and rate constants, since it provides a more controlled treatment of the resonances than perturbation theory (Bennett et al., 2003). According to Breit-Wigner theory (Breit and Wigner, 1936;Bain and Bardsley, 1972), if the interference effects between the direct and resonance mechanisms may be neglected, then the total RA cross section and the rate coefficient can be decomposed as where σ res E c ( ) and k res T ( ) are the resonance contributions and σ dir E c ( ) and k dir T ( ) are the direct contributions. Since the direct contribution does not contain resonances by definition, σ dir E c ( ) can be identified as the baseline of the exact cross section which can be obtained either by perturbation theory (using Eq. 1) or it can be estimated by the classical trajectory method (Gustafsson, 2013;Nyman et al., 2015). The resonance contribution to the cross section and the rate coefficient can be obtained as  σ where v and j are vibrational and rotational quantum numbers of the quasi-bound states. ϵ j,v is the energy of the quasi-bound state that corresponds to the peak position of resonances. Γ tun v,j is the width associated with tunneling and Γ rad v,j is the width due to the radiative decay to all lower-lying Frontiers in Astronomy and Space Sciences | www.frontiersin.org July 2021 | Volume 8 | Article 704953 bound and quasi-bound levels j′, v′, of the final electronic state. The resonance contribution to the rate constant in Eq. 7, according to the conventional Breit-Wigner formula, assumes that the quasi-bound states are populated only by tunneling through the potential barrier. This is appropriate in the low gas density and low radiation field limit (Forrey et al., 2018). If there are other modes of populating the quasi-bound states (interaction with the radiation field or third body collisions), an upper limit for the resonance contribution is obtained by replacing Eq. 7 by its local thermal equilibrium (LTE) approximation (Gustafsson and Forrey, 2019) where the sum runs over the j and v quantum numbers of the quasi-bound states.

Computational Details
The electronic structure calculations were carried out with the MOLPRO package (Werner, 2015). The potential energy curve used in the dynamical calculations was calculated with the explicit correlated internally contracted multireference configuration interaction method (icMRCI-F12) with Davidson correction, using the aug-cc-pVQZ-F12 basis set as implemented in MOLPRO. All calculations were carried out in the C 2v symmetry group. The molecular orbitals were constructed using the state averaged complete active space self-consistent field (CASSCF) method with an active space consisting of four electrons on eight orbitals (6a 1 , 1b 1 , 1b 2 , 0a 2 ). Five states with A 1 irreducible representation were included in the state average with equal weights in order to obtain a stable wavefunction at each configuration. For the calculation of the permanent dipole function we employed the full configuration interaction (FCI) method with aug-cc-pVQZ basis set. The dipole was calculated as an expectation value of the dipole operator. Dynamical calculations of radiative association cross sections were carried out according to the quantum mechanical perturbation theory (PT) and Breit-Wigner (BW) method. In the PT treatment the vibrational Hamiltonian was represented with a DVR basis (Colbert and Miller, 1992) for the calculations of bound states. The Numerov method was used for the calculation of the scattering wavefunctions. The total angular momentum ranged from J 0 to 60 and the energy range was E c 10 -5 eV-8.0 eV. An energy resolution of ΔE 10-6 eV was used in the entire energy range so that narrow resonances are reasonably well represented. The resonance parameters, which are needed in the BW method, were computed with the program LEVEL (Le Roy, 2007). The value of the statistical weight factor is f stat 1/4 for reaction (9) and the non-zero Hönl-London factors are S J→J+1 J + 1 and S J→J−1 J (Hansson and Watson, 2005).

RESULTS AND DISCUSSION
There is only one potential energy curve, that of the ground state ( 1 Σ + , see Figure 1), which can provide an available pathway for the molecule formation in the following reaction The corresponding permanent dipole moment in the center of mass fixed frame of reference, and the dipole after isotope substitution, are also displayed in Figure 1. (See Appendix A for details on the relation between the dipole moments of the BeH + and the BeD + molecules). It is clear from Figure 1 that the dipoles correlate correctly to the free Be + + H/D limit as the atomic distance goes to infinity. The ab initio data of the potential and dipole curves are tabulated in Appendix B.
The cross sections for the BeH + molecule calculated with the quantum mechanical perturbation theory are shown in Figure 2A. Due to the attractive potential energy curve and the deep potential well, the efficiency of the RA process is more dominant at lower collision energies and it is still significant up to E c 1.0 eV. Above this energy value, nonreactive processes may dominate in the collisions: bremsstrahlung under 10.0 eV and ionization and electronic excitation of the atoms above (Kramida et al., 2017). The cross section curves show a rich resonance structure, especially above 10 -3 eV where the very narrow resonance peaks are close to each other. Several of these resonances contain contributions from quasi-bound to quasi-bound transitions. Nevertheless, these transitions may have a significant contribution to the radiative association of colliding fragments (Bennett et al., 2003). This unusual way of molecule formation can be efficient, since the decay into the continuum manifold can be slower than the radiative relaxation to the stable ro-vibrational quantum states. The isotope substitution slightly changes the dynamics of the RA (see Figure 2B): the base-line of the cross section is diminished, FIGURE 1 | Potential energy (black) and electric dipole curves (green and blue) for reaction (9) used in the dynamical calculations.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org July 2021 | Volume 8 | Article 704953 according to a known scaling law (Gustafsson, 2013), and the resonance structure is affected as well.
The PT method has some peculiarities in the description of the resonances. The evaluation of the cross section even for a diatomic molecule can be computationally challenging at high energies. Even with a very fine energy resolution we cannot describe all of the narrow resonances, especially above 0.1 eV (see Figure 2A). This means that many of the resonances are missing and others are not perfectly resolved. Owing to this, the corresponding rate constant, computed by numerical integration of (Eq. 3), may not be converged. Furthermore, the PT method does not account for radiative broadening of the quasibound states (Le Roy, 2007;Zamecnikova et al., 2020). Therefore, it will give resonances in the cross section where the total width Γ tun v,j + Γ rad v,j in Eq. 6 is replaced by Γ tun v,j . Even if a rate constant was computed from the perfectly resolved PT cross section it would thus result in the LTE rate constant (Gustafsson and Forrey, 2019).
To bypass the shortcomings of the PT method, we also calculated the cross sections with the Breit-Wigner method (see Figure 2B), which makes possible the direct evaluation of the rate constant, Eq. 7. In our calculations we used the baseline of the PT cross sections as the direct contribution to the BW method (Eqs. 4, 5). Figure 2A displays how the BW method provides a comprehensive and accurate description of the resonances in the entire energy range, especially at high energies where the PT method fails to provide resonances.
The rate constants for the formation of BeH + and BeD + molecules are displayed in Figure 3. As expected from the cross sections ( Figure 2B), the rate constant of the BeD + is smaller with a factor of 2-3, compared to those of BeH + using all of the methods. Due to the insufficient description of the resonances in the PT method, it underestimates the rate constants in the whole temperature range. The magnitude of the rate constants is comparable to formation by RA of other diatomic molecules , but it is still much smaller than rates of common chemical exchange reactions. Based on this, the RA mechanism is likely negligible in a higher density astrochemical environment, where complex chemical networks dominate the molecule production owing to the high diversity and density of chemical species. Nevertheless, RA might be an important reaction channel in the low density, dust-poor regions of the interstellar space, since it is a direct and sometimes the only route to the molecule formation.
Since the production of BeH + and BeD + may be important in a variety of environments, as discussed above, we have also included LTE rate constants in Figure 3. Those are larger than the conventional BW rate constant by up to a factor of seven. This may have to be taken into account in studies of astrophysical environments with high density or, in particular, strong radiation fields. It may also play a role in studies of fusion reactors.  Following the recommendation of (Antipov et al., 2009), the BW cross section is used to obtain the rate constant for the formation of the BeH + and BeD + molecules in the interstellar medium. The extended Arrhenius equation is used to fit the BW rate constants. Owing to the complex structure of the reaction rates, the corresponding curves are split into several temperature ranges which are fitted individually. The fitting parameters are summarized for the formation of BeH+ and BeD + in Table 1 according to different temperature ranges. These rate constants will be made available in the database KIDA (Gustafsson et al., 2012;Wakelam et al., 2012).

CONCLUSION
In this work, we studied the radiative association formation of the BeH + and BeD + molecules. The cross sections and the rate constants are calculated by using quantum mechanical perturbation and the Breit-Wigner theory. The obtained data can be utilized in the modelling of the BeH + abundance in the interstellar medium and for its production in thermonuclear fusion reactors. Furthermore, the resonance contribution of the rate constant under local thermal equilibrium conditions is calculated. Such conditions can be particularly important in fusion reactors, as well as in astrochemical environments with strong radiation fields.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
PS and MG contributed to conception and design of the study. PS and GS performed the calculations. PS and MG wrote the manuscript. All authors contributed to manuscript, read, and approved the submitted version.

ACKNOWLEDGMENTS
This article is based upon work from COST Action CA18212-Molecular Dynamics in the GAS phase (MD-GAS), supported by COST (European Cooperation in Science and Technology). Support from the Kempe Foundation is gratefully acknowledged. MG acknowledges support from the Knut and Alice Wallenberg Foundation. Computational resources provided by the Swedish National Infrastructure for Computing (SNIC) at HPC2N are acknowledged.