Abstract
The title reaction is studied using a quasi-classical trajectory method for collision energies between 0.1 meV and 10 eV, considering the vibrational excitation of reactant. A new potential energy surface is developed based on a Neural Network many body correction of a triatomics-in-molecules potential, which significantly improves the accuracy of the potential up to energies of 17 eV, higher than in other previous fits. The effect of the fit accuracy and the non-adiabatic transitions on the dynamics are analysed in detail. The reaction cross section for collision energies above 1 eV increases significantly with the increasing of the vibrational excitation of (), for values up to =6. The total reaction cross section (including the double fragmentation channel) obtained for =6 matches the new experimental results obtained by Savic, Schlemmer and Gerlich [Chem. Phys. Chem. 21 (13), 1429.1435 (2020). doi:10.1002/cphc.v21.13]. The differences among several experimental setups, for collision energies above 1 eV, showing cross sections scattered/dispersed over a rather wide interval, can be explained by the differences in the vibrational excitations obtained in the formation of reactants. On the contrary, for collision energies below 1 eV, the cross section is determined by the long range behaviour of the potential and do not depend strongly on the vibrational state of . In addition in this study, the calculated reaction cross sections are used in a plasma model and compared with previous results. We conclude that the efficiency of the formation of in the plasma is affected by the potential energy surface used.
1. Introduction
Hydrogen, as the most abundant element in the Universe, plays a fundamental role in star formation and the chemical evolution of molecular Universe. Its molecular forms are , and [Citation1]. In evolved galaxies, the formation of is usually attributed to atomic hydrogen recombination on cosmic grains and ices [Citation2–4]. is rapidly formed by cosmic rays or electrons, and it collides with to form in the reaction (1) (1) Once the molecular forms of hydrogen are formed, the chemistry in space starts with the formation of the first hydrides. In cold clouds, the most abundant ion is , which is considered to be the universal protonator [Citation5–8] through the proton hop reaction [Citation1, Citation9] (2) (2) where M is an atom or a molecule. The HM cations are very reactive and trigger many chemical networks giving rise to most of the molecular systems detected in space [Citation5–8, Citation10, Citation11]. In cold environments, M colliders in Equation (Equation2(2) (2) ) deposit on ices, and reacts only with the most abundant molecule, , as (3) (3) a proton exchange reaction, which is constrained by nuclear spin statistic [Citation12–18], and is the responsible for the ortho/para ratio of . When the collider is HD, reaction (Equation3(3) (3) ) is the responsible of deuteration [Citation13, Citation17, Citation19–24]. The deuterated species formed following Equation (Equation2(2) (2) ) produce the observed high relative abundance of deuterated species, estimated as ≈ 10 times higher than the D/H ratio of the galaxy [Citation25–28]. This high deuteration efficiency is attributed to the zero-point energy differences among the deuterated species, very significant at the low temperatures of cold molecular clouds [Citation22, Citation29].
Hydrogen plasma [Citation30–32], apart from their technological applications in industry, medicine and fusion reactors [Citation33, Citation34], can be considered as a prototype for Early Universe models [Citation2, Citation35–38]. In the absence of other species but hydrogen, the molecular species are simply , and , but also exists. The distribution of the ions affects the determination of hydrogen particle flux in the plasma [Citation39, Citation40]. The reaction in Equation (Equation1(1) (1) ) is an important process for modelling plasma at low electron and molecular temperatures ( eV, eV) since the process is the only dominant process for the formation of [Citation41]. Therefore, an accurate cross section is essential for the determination of the density of in plasma models. Moreover, the rovibrationally resolved cross sections are required for collisional-radiative (CR) spectroscopic modelling of molecular hydrogen which can be applied to a fusion detached plasma [Citation42, Citation43]. The isotopic effect on the reaction is also essential for plasma modelling in nuclear fusion tokamak [Citation44].
The reaction in Equation (Equation1(1) (1) ) has been widely studied experimentally [Citation45–54] and theoretically [Citation49, Citation55–60]. In the 1 meV and 1 eV energy interval, an excellent agreement between experimental [Citation54] and theoretical [Citation60] reactive cross sections has been achieved. Recently, this reaction has been studied in two extreme regimes, at ultra cold energies [Citation61–64] and at high collision energies [Citation65] up to 10 eV. At ultra cold energies, many isotopic variants have been studied, and the reactive cross section shows a Langevin-like behaviour so that the measured reactive cross sections, in relative units, shows also excellent agreement with the theoretical simulations [Citation60]. However, recent experiments performed at higher collision energies (1-10 eV) [Citation65] differ considerably from previous theoretical simulations [Citation49, Citation55, Citation60]. Moreover, for collision energies above 1 eV, the new experimental measurements show differences with previous experimental ones [Citation45–48, Citation52, Citation53], all of them scattered in a rather wide interval of cross sections.
There are two main goals in this work. First, we focus on the theoretical simulations of the formation reaction, Equation (Equation1(1) (1) ), to reproduce the new experimental data above 1 eV [Citation65]. Second, we study hydrogen or deuterium plasma models in order to analyse the effects of the calculated cross sections and rates on the or D densities.
This work is organised as follows. In Section 2, we develop new potential energy surfaces (PESs), which include non-adiabatic effects and increase the accuracy of the fit. This new fit uses a Neural Network (NN) method to describe the four-body term, to improve a zero-order description using a Triatomics-in-Molecule treatment, which accurately fits very precise ab initio calculations, over a broader energy interval than previous fits, up to 17 eV. In Section 3, we study the reaction dynamics using a quasi-classical trajectory (QCT) method, including transitions among different electronic states, using the fewest switches method of Tully [Citation66]. The reaction dynamics is studied for collision energies between 1 meV and 10 eV, and for several vibrational states of reactant, as well as for the deuterated reaction, focussing on high energy reactive cross sections recently measured by Savic et al. [Citation65]. In Section 4, we investigate how the calculated reactive cross sections affect the population density of and D in CR models of /D plasma. Finally, in Section 5 some conclusions are extracted.
2. Potential energy surfaces of
In this work several PESs are used, and are listed here to clarify the differences:
PES1: | This PES developed by Sanz-Sanz et al. [Citation67], is the most accurate developed so far for this system. It is built as a sum of a triatomics-in-molecules (TRIM) term, , plus a four body term, . The TRIM term is a generalisation of the Diatomics-in-Molecules (DIM) [Citation68–71] method, in which the electronic Hamiltonian is factorised as a sum of triatomic and diatomic fragments as [Citation67, Citation72] (4) (4) where are the monoelectronic Hamiltonians of fragments and are the bielectronic Hamiltonians (for n−i, o−i electrons) describing the system for the ino nuclei. The TRIM representation consists of a 8×8 matrix, whose elements have and matrix elements of different electronic configurations, in which each hydrogen atom is described by a 1s function, except one corresponding to . In PES1 [Citation67], the triatomic terms included in the TRIM matrix were built as the 3×3 DIM matrix for each fragment plus a three-body term added to describe the ground state of either singlet or triplet symmetry. The non-adiabatic terms in these triatomic fragments are therefore approximated by those obtained in a DIM treatment. PES1 is built adding a four-body correction (MB term) to this TRIM description. This MB term is expressed as a linear combination of four-body polynomials that are symmetric with respect to the permutations of identical nuclei (Permutationally invariant polynomials or PIP). These PIP are built in terms of Rydberg polynomials of the interatomic distances () [Citation73–75]. The set of linear and non linear coefficients, , are optimised to minimise the difference with the ab initio energies obtained with the multi reference configuration interaction (MRCI) method using the aug-cc-pV5Z basis set [Citation76]. In PES1, the same four-body correction term is added to all the diagonal elements of the TRIM matrix. The reactive cross sections calculated on this PES1 shows an excellent agreement with the experimental results for collision energies between 1 meV and 1 eV [Citation60–62, Citation65]. | ||||
PESTRIM8×8 and PESTRIM1×1: | This PES is based on the TRIM model explained above, and is an improvement made in this work. The main difference is that the triatomic is represented by 3×3 diabatic matrices fitting its three lowest singlet and triplet electronic states [Citation77] to MRCI ab initio energies extrapolated to the complete basis set (CBS) limit [Citation78]. This improvement is crucial to study the non-adiabatic dynamics of . Analytical derivatives of the potential energy surfaces and non-adiabatic couplings are calculated based on the Hellmann-Feynmann theorem as described in [Citation60]. The full PESTRIM8×8 potential correspond to the eight adiabatic eigenvalues including the non-adiabatic coupling terms. PESTRIM1×1 only considers the ground adiabatic energy. | ||||
PES-NN: | PESTRIM1×1 lacks high accuracy in the interaction region, i.e. where is formed, while showing excellent agreement for or + rearrangement channels. For this reason a four-body neural network term is added to PESTRIM1×1 to improve the accuracy of the system, as described below. |
2.1. Many body neural network term
Many body potential energy terms are widely used to represent potential energy surface of chemical systems, either in a many body expansion [Citation74, Citation79, Citation80], where one up to N-body terms are summed, or as correction terms of a zero-order description of the potential, described by the TRIM method [Citation67, Citation72] or by a reactive force field matrix [Citation81–83].
In this work a many body neural network (MB-NN) [Citation83] potential energy term has been built as a PIP-NN [Citation84], in which the neural network is fed with a permutational invariant polynomial representation of the molecular geometry. A PIP is constructed by projecting a polynomial of the interatomic distances into the totally symmetric irreducible representation of the desired permutation group. A generator of these PIPs is defined as: (5) (5) where is the number of interatomic distances, is a function of the interatomic distance , is the exponent of the ith monomial for the nth polynomial and is the projector to the totally symmetric irreducible representation of the permutation group. A common choice of in PIP-NN-PES is the decaying exponential .
The set of PIP has to be carefully filtered so that it is purely composed of N-body functions. This means that any of these functions evaluated on a geometry where the N bodies are not closely interacting should be zero. In case these terms included any lower body functions, we have shown that it would add spurious interactions between fragments [Citation85], which specially affects the long range regions.
In Appendix A it is shown that a polynomial in Equation (Equation5(5) (5) ) can be represented as a graph, and that the subset of N-body polynomials corresponds to those whose respective graph is connected, meaning that there exists a path which connects all the vertices (particles) of the graph. In this way we can guarantee that any N-body polynomial, with defined as a decaying exponential or Rydberg function, will tend to zero as any particle or set of particles moves away from the rest. This will automatically make zero a PIP-PES and provide constant descriptor for a NN-PIP-PES, that returns a constant energy out of the N-body region, which can be trained to be as close to zero as possible and that produces no net force, since its derivative with respect to the atom coordinates is zero.
The four-body term developed here is expressed as a feed forward neural net with 11 input neurons and two hidden layers with = 32 and = 77, and sigmoid non linearities : (6) (6) were matrix (with elements ) and vector (with elements ) represent the trainable weights and bias on layer l. The 11 input neurons correspond to four-body PIPs produced by setting a maximum polynomial degree max() = 5 and maximum monomial degree max()=2. All lower body polynomials that would introduce spurious energy contributions in reactants and products channels are filtered following the steps detailed in Appendix A.
The four-body term is trained with NeuralPES, an in-house Python code based on PyTorch [Citation86]. New ab initio points have been used in the fit of this work, of higher accuracy of those used in PES1 [Citation67]. The energies are obtained using a two point extrapolation method to complete basis set (CBS) [Citation78], using the results obtained with the aug-cc-pV5Z and aug-cc-pV6Z basis sets. Around 33000 ab initio points were calculated, including the geometries of Ref. [Citation67] and new ones selected to increase the accuracy of the PES above 2 eV. These last new points were chosen from QCT trajectories at different collision energies (from 1 eV to 10 eV) taken on the TRIMPES1×1 to populate physically accessible configurations at this high energy regions. The complete set of ab initio points is randomly split into a training set (containing 80% of the data) and validation and test sets (with 10% of the points each). The training set consists on 27294 geometries, with energies up to 17 eV over the + asymptote, mostly corresponding to four-body interactions, and elongations to lower body geometries. The training process aims to minimise the root mean squared error between the ab initio and PESTRIM1×1 + energies.
2.2. Analysis of the different PESs
As all the four-body terms described above vanish as the systems tend to reactants or product asymptotes, the long-range interactions are purely described by the triatomic terms of the TRIM model. In all the triatomic fits considered here, the long range terms for + and fragments, for either singlet or triplet states are very precisely described [Citation77, Citation87, Citation88]. This produces highly accurate long-range interaction in the + channel [Citation60]. Allmendinger et al. [Citation61, Citation62] applied a new experimental set up to study + reaction at low temperatures. The measured cross section was then scaled to reproduce the cross section calculated with the PES1 potential [Citation60] at a single collision energy. The excellent agreement between calculated and scaled experimental cross sections between 0.5 and 5 meV demonstrates the good behaviour of the long-range interactions included in PES1. The new PESs introduced in this work, describes the long-range interaction even more accurately, by using improved long-range interactions in the triatomic fragments [Citation77]. The RMS errors for the different potential energy surfaces are presented in Table , in different energy intervals. The improvement of PES-NN over PESTRIM1×1 is clear in all energy ranges due to the enhancement of the channels as presented in Table . PES1 and PES-NN show comparable RMSE for energies below 2 eV, but when ab initio points above 2 eV are added, PES-NN shows a much higher accuracy.
The whole configuration space is divided as follows: when all the interatomic distances are below Å, the system is taken to be in region; if all interatomic distances of a triad of hydrogens are shorter than , the region is taken to be ; if any two pairs of atoms present interatomic distances shorter than , the region is denoted as + ; if only one internuclear distance is , the region is called ; otherwise, if all the interatomic distances are large, the system is in fully dissociated. Following this division, the top left panel of Figure shows the distribution of the data set among the different regions as defined above, as a function of the energy, taking the + reactants asymptote as the zero of energy. Most of the data correspond to and geometries which connect this channel to + and . As can be seen in the top right and bottom panels, the PESTRIM1×1 is highly accurate in the latter two channels, but shows deviations in the region, up to several eV. The effect of the four-body term on PES-NN is decisive for the proper description of region. PES1 was fitted for energies up to ≈ 2 eV, and it presents large energy deviations over this threshold. On the contrary, PES-NN keeps a high precision at high energies, which are the interest of the present work.
In Figure , + approaches for different and angles in their equilibrium geometry are shown for the three potential energy surfaces and compared with ab initio calculations. PESTRIM1×1 tends to predict a larger energy for geometries, while PES1 and PES-NN yield effectively the same description. As the interfragment distance increases towards + the four-body terms in both PES1 and PES-NN go to zero and all that remains are the respective TRIM terms. The improvement of the new PES-NN is better seen when the diatomic fragments are not at the equilibrium geometry, for energies above 2 eV.
The more accurate description of the higher energy regime, energies larger than 2 eV, can be seen in Figure where H approaches to a compressed . The differences between PES1 and PES-NN are more pronounced when bonds are compressed than stretched. This is
3. Reaction dynamics
3.1. Quasi-classical trajectory and surface hoping calculations
The quasi-classical trajectory calculations are performed with the MDwQT code [Citation60, Citation82, Citation89, Citation90]. When considering several coupled adiabatic electronic states, the fewest switches approach of Tully [Citation66] is used as described in [Citation60]. Initial conditions are sampled with the usual Monte Carlo method [Citation91]. The initial conditions for the vibrational modes of and are quantised using the adiabatic switching method [Citation92–94], yielding vibrational energies within 0.3 meV with respect to the exact vibrational levels of each diatomic fragment. The rotational states of and are set to zero in these studies, and the initial distance between the two centre-of-mass is set to 105 bohr. The initial impact parameter, b, is sampled between 0 and B, according to a quadratic distribution on b, where B is determined for each energy according to a capture model [Citation95] as , for a charge-induced-dipole interaction, described by , where α is a constant, with dimensions , which is proportional to the average polarisability of , β, at its equilibrium configuration, as . All trajectories are stopped when any internuclear distance becomes longer than 125 bohr, where they are analysed.
The reactive cross section for each collision energy, E, is calculated as [Citation91] (7) (7) where is the number of trajectories leading to products and is the total number of trajectories with initial impact parameter lower than , the maximum impact parameter for which reaction takes place at energy E. Here we have considered (v=0), while vibrational level varies between 0 and 6. For each () couple and each energy a set of 10 trajectories are run, with energy error lower than 0.01 meV.
The final energy distribution of products is also analysed. We simply evaluate classical energies, without trying to consider the permutation symmetry of identical fermions (for ) or bosons (for D). This is done in three steps. First, the kinetic energy of and H products are calculated and substracted. Second, the rotational angular momentum of products is evaluated, and its rotational energy. By setting the origin of energy at the bottom of the well, the remaining energy corresponds to vibrational energy. Here, we do not attempt to assign the internal vibrational modes, which deserves further development and is led for a future work.
3.2. H + H() collisions in PES1
The new experimental results for the title reaction of Savic et al. [Citation65] differ from previous ones, both experimental and theoretical, above 2 eV. The difference with previous experimental data, which are scattered above 1-2 eV, may be due to different conditions in the generation of the reactants. In low temperature plasma the vibrational temperature of is of the order of 2500 K, so that the population of (v=1) is expected to be lower than 10% [Citation30–32]. On the contrary, is formed by electronic impact or photoionization, which may yield to different vibrational and rotational excitations. Vibrationally excited can partially thermalise, yielding to different initial conditions in different experimental setups. As an example, recent theoretical calculations [Citation96, Citation97] on the charge transfer reaction, and some isotopic variants, have found that the reaction cross section highly depends on the initial vibrational state of the diatomic ion. Following this idea, in this work we have performed QCT calculations of the cross section of the reaction for several initial vibrational states of reactant using the global PES1 of Refs. [Citation60, Citation67], which are shown in Figure .
For exothermic reactions without a barrier, the long-range interaction between the reactants dominates the reaction dynamics. In the case of charge-induced dipole long-range interactions the cross section for exothermic reactions takes the form [Citation95] (8) (8) This is approximately the behaviour of the cross sections for energies below 1 eV for every initial vibrational state . Therefore, we can conclude that in this energy interval reaction dynamics is dominated by long range interactions, independently of the initial vibrational state of ().
However, above 1 eV the reactive cross sections present important differences among the considered, showing that the vibrational excitation has a strong impact on the reactivity. In general, the reactive cross section increases with increasing , which could be simply understood assuming that can break more easily. In the left panels of Figure the two main mechanisms to form are separated as H-hop and proton-hop, corresponding to the fragmentation of or reactants, respectively. For collision energies below 1 eV, the cross section for the proton-hop mechanism is slightly larger. However, above 3 eV the H-hop cross section is larger for =0, and the two mechanisms tend to a rather similar value as increases. This means that the vibrational excitation of does not produce significant increase of the proton-hop mechanism. This can be explained looking at the right panels of Figure , where the maximum impact parameter, , is plotted for each proccess and inital vibrational state. Above 1 eV, increases from = 0 to =6, in nearly an identical quantity for the two mechanisms. We conclude that the increase of the cross section when varying is due to the growing of the subunit, whose right turning point increases from 1.25 to 2.12 Å, for =0 and 6 respectively. Therefore, at energies above 1 eV, when long-range interactions are not able to produce important deviations among the reactants, the size of the two reactants approximately determines the value of the maximum impact parameter. At these higher energies a more complex reaction mechanism occurs, in which may nearly insert in the bond, specially at high excitations.
3.3. D + D collisions in PES1
The reactive cross section for the D(v, j) + D() collisions is shown in Figure for v = j = 0 and = = 0. For energies below 1 eV, the results for D + D closely match those for + . This can be explained in terms of the Langevin model, in which the cross section of Equation (Equation8(8) (8) ) does not depend on the mass of the reactants. This explains why the reaction cross sections for D and are nearly the same.
However, it should be noted that for reactions involving partially deuterated species, the reaction cross section presents larger differences, as those already reported in the theoretical study of Ref. [Citation60]. The shift of the centre-of-mass with respect to the geometric centre of the two diatomic reagents introduces important differences, specially related to the effect of the rotation. In particular, in the homonuclear neutral or D case for j=0, only the charge-electric dipole term affects the long-range interaction, determining the reactive cross section below 1 eV. In this case, the charge-electric quadrupole term vanishes when integrating over the angular coordinate for the isotropic j=0 case (but not for 0).
3.4. Non-adiabatic effects and fit accuracy
The increasing of the vibrational excitation yields to an increase in the reactive cross section. However, this increasing, even for =6 does not match the new experimental data by Savic et al. [Citation65]. Since the cross section as a function of the excitation seems to converge to the value of =6, here after we shall focus on the two limiting cases, = 0 and 6.
In order to investigate the role of electronic transitions among the lower electronic states, in this work we use the PESTRIM8×8, comparing it to the results on the PESTRIM1×1 model. The dynamical results obtained for these two potentials are shown in Figure , and compared with the data obtained with PES1. All the results present a similar behaviour, with small differences in the logarithmic scale used in the figure. All the results converge to the same value at the lowest collision energy considered here, 1 meV, since all the PESs used in this work have similar long range interactions.
At intermediate energies, between 0.01 and 2 eV, the PESTRIM8×8 and PESTRIM1×1 cross sections are slightly different, showing a small effect of non-adiabatic transitions in the dynamics. Curiously, these non-adiabatic effects seem to produce a decrease on the reactive cross section, in contrast to what is needed to match the values at higher energies of the experimental results of Savic et al. [Citation65].
In Figure we also compare with the cross section obtained with the adiabatic PES1, which nearly matches the results obtained with the adiabatic PESTRIM1×1 up to 3 eV. Above 3 eV, the results obtained with PES1 are higher than those for PESTRIM1×1, showing that the four body term may be important. However, for the new PES-NN, of higher accuracy than PES1, the reaction cross section is nearly identical to that of PESTRIM1×1 for collision energies below 1 eV and very close to those obtained with PESTRIM8×8 for energies above 4 eV. From this comparison we may conclude that the better accuracy of the four body term of the new PES-NN does not improve significantly the differences between the simulated and measured [Citation65] reactive cross-sections.
In order to analyse whether there are other excited states not well described by the TRIM approximation, we have performed ab initio calculations of the four lower electronic states, considering a larger electronic basis set, with extra orbitals added to describe the 2s and 2p electronic states of atomic hydrogen. With this larger basis set, we have found that the energies obtained (extrapolated to the complete basis set) differ only a few tenths of meV with the previous ab initio calculations, and that no higher electronic state appears below 10 eV.
3.5. Double fragmentation and reaction mechanism
For the double fragmentation (DF) channel opens at ≈ 2 eV, as it is shown in the left panel of Figure . The opening of this channel occurs approximately for values where the extra energy of =6, 1.67 eV, is added to the D of , 2.65 eV
If we add the DF channel contribution to the production of , as shown in the right panel of Figure the cross section increases considerably, reaching a very good agreement with the new experimental data of Savic et al. [Citation65] above 3 eV.
The main mechanism giving rise to the DF channel consists of three steps, as it is shown in the left panels of Figure . First, a highly vibrationally excited intermediate is formed by insertion of in the elongated , which lives short time. Second, a first H atom is ejected (in the Figure is atom 2 with charge 0), forming a very excited (), which lives from 80 to 160 fs, approximately. In the third step, one of the atoms of (atom 4 in the Figure) dissociates, carrying the positive charge, thus leading to neutral . Since this third step occurs much later, it could explain that experimentally the metastable () would be detected together with more stable products. It is also important to notice that the energy transfer among particles of identical mass is possibly overestimated in a classical treatment, as the one used in this work. If this is the case, it would support the inclussion of the DF cross section as a part of the total reaction cross section. Quantum calculations are needed to solve this problem, but they are rather challengig at the high energy considered.
The right panels of Figure show another DF mechanism: in this case and are produced at 70-80 fs. is vibrationally very excited and dissociates later, at ≈100–110 fs. In this case, there are two degenerate electronic states, each one corresponding to the charge in one of the ejected atoms, at long distances from the fragment. In the ground electronic state, shown in the lower panels, the charge is exchanged between the two identical atoms, because of this degeneracy, and shows the nature of the surface hopping occurring in the products channel when including several electronic states (PESTRIM8×8). The electronic transition occurs among degenerate electronic states describing + and + products, what explains the small effect of including electronic transitions in the reaction dynamics.
In addition, the charge transfer in the entrance channel occurs between and , when the two reactants have the same internuclear distance (see bottom panels of Figure ). In this situation there is a degeneracy between the two lower adiabatic states, as discussed in detail in Ref. [Citation60].
The energy distributions of the products are shown in Figure , and it is nearly quantitatively the same for all the PESs. The energy difference between the two vibrational states of (= 0 and 6) is approximately 1.4 eV, close to the exoergicity. For low collision energies, 0.1-1 meV, the initial vibrational energy of ( = 6) is 1.40 eV higher than that of (= 0), and the vibrational energy of the corresponding products is also higher, but only by ≈ 0.9 eV. Therefore the remaining 0.4 eV are nearly equally distributed between the rotational and translational energies of the products. Rotational energy increases with collision energy, except for the = 6 above 1 eV, where rotational excitation reaches a plateau and seems to start decreasing. The translational energy always increases for = 0, while for = 6 it slightly decreases below 0.1 eV, and then increases again. The vibrational energy of products shows two different behaviours: below ≈ 0.7 eV, the vibrational energy slightly decreases (=0) or remains constant (=6); above this energy the vibrational energy increases sharply.
Such behaviour allows to assign two different reaction mechanisms. Below ≈ 0.7 eV, the impact parameter (in Figure ) is large, supporting a stripping mechanism, in which the long range interactions attract the reactants to each other, originating orbits enhancing the relative angular momentum between the two reactants. Above ≈ 0.7 eV, however, the impact parameter is rather small, and the reaction occurs by an insertion mechanism, in which the is greatly excited vibrationaly, specially as initial vibrational and translational energy of the reactants increases.
The dissociation energy of is 4.34 eV, very close to the value reported by [Citation99] of 4.35 eV, and the average vibrational energy distribution of reaches values in the 4-5 eV interval, values above the dissociation energy explaining why dissociates, leading to the DF channel.
4. Plasma modelling
To model the hydrogen plasma, we shall consider a gas-discharge vessel of cylindrical symmetry, so that only z, parallel to the central axis, and R, the distance from the centre to the walls, will be considered, with the cylinder of infinite length in this case. The gas is initially in the form of neutral , and after the discharge ignition new species are formed, H, , , and electrons. is neglected under the conditions considered here. To model the abundance of these species in the stationary condition we use a model similar to those already described previously [Citation41, Citation100], which is outlined in the Appendix B, including all the processes listed in Table , in which the vibrations of molecular species, , and are not considered.
The plasma model is done for pure hydrogen and pure deuterium gases. For modelling deuterium plasma, the cross sections for electron collisions and radiative transitions of D species are used by those of H species.
The cross section for + and D + D reactive collisions calculated in this work are included in the models presented below. The plasma modellings are also done using the cross section by Janev et al. (2003) [Citation98], which are compared with the present calculations in Figure .
Here, the electron temperature eV and electron density cm are used, which were measured by a Langmuir probe for D plasma in a long cylindrical vessel [Citation39]. The molecular temperature is set to = 0.026 eV (300 K) and the atomic temperature is = 0.052 eV (600 K). At these realistic temperatures, the relevant energies are below = 1 eV, so that the vibrational excitation has no significant effect.
Figure shows the resulting population densities of D and H species depending on the cross section used. We can note, the and population densities vary significantly with the cross section used, while other species are nearly unchanged, as shown in Figure . The D and population density changes are also summarised in Table . It is worth noting that the main depopulating mechanism for density is the diffusion process (given in the last term of Equation EquationA7(A7) (A7) ), while the electron impact processes ( in Table ) only contribute to the depopulation in a small fraction.
The present cross section for D + D is larger than that for + by % at low collision energy, below eV. On the contrary for collision energies above 1 eV, the cross section for deuterium is ≈ 20 lower than that of pure hydrogen reaction, as shown in Figure .b. For + , the cross section by Janev et. al [Citation98] is smaller than the present ones by % below 1 eV, but the difference is enlarged up to 100 % (black triangle) at the collision energy over eV as shown in Figure (b).
The reaction rate coefficient calculated at the molecular temperature, eV, for pure deuterium is ≈ 20 larger than for pure hydrogen, and is also larger than those obtained from Janev et al. [Citation98] by %, as listed in Table . These differences have a rather linear impact on the resulting D and population densities, whose difference varies proportionally to the difference between the rate coefficients listed in Table .
As a result, the use of the present results for D + D and + leads to significantly different D and population densities compared with the widely used cross section for + [Citation98] in this plasma modelling. It should also be noted that the use of the cross section for + instead of that for D + D in the modelling of D plasma can give rise to unreliable population density of D, even though the difference between the two cross sections is small, due to the rate coefficient sensitivity to the cross section at the low collision energy.
When T and n are reduced to 3 eV and 1.5 × 10 (cm), respectively, and the molecular pressure is increased by about 10 times, the amount of X (X=D, H) becomes dominant over those of X ions. These conditions are similar to those reported by Tanarro and Herrero [Citation32], who also find a significant increase of the population. This relative increase of the X density is mostly attributed to that the rate coefficient of X + e collision populating X becomes smaller than the rate coefficient of X + X collision depopulating X.
However, the density of X in this high pressure is less sensitive to the change of the cross section for X + X collision than in the low pressure. This is due to the fact that the increased X + X rate coefficient is accompanied by the decreased X quantity since X + X collision is the main depopulation process for X, which leads to the little change of the X formation. While in the previous case of low pressure, the main depopulating of X does not come from the X + X collision but from X + e collision and the density of X is not affected by the X + X rate coefficient. Thus the density of X in the low pressure case is more sensitive to the change of the cross section for X + X collision than in the high pressure limit having its linear dependency on the rate coefficient mentioned above.
On the other hand, when the molecular temperature is increased from = 0.026 eV (300 K) to = 4.2 eV (50,000 K), the rate coefficient for the present cross section of (v=0,j=0) + (= 0,= 0) differs from that for the cross section by Janev et al. [Citation98] only by about 20 %. However, the cross section of (v=0,j=0) + (=6,=0), shown in Figure , is much larger than that of (v=0,j=0) + (=0,=0) by Janev et al. [Citation98] at collision energies higher than 1 eV. The difference between the rate coefficients using these two cross sections is as much as 10-170 % at the temperature range = 0.026-4.0 eV. The larger differences are found at higher , becoming over 40 % above 1 eV. Hence the high state of can contribute to the population of much more than the =0 state. This was analysed by replacing the (v=0,j=0) + (=0,=0) reactive rate coefficient by that obtained for (v=0,j=0) + (=6,=0) collision. To further analyse the effect of vibrational state, new plasma models need to be developed, increasing considerably in complexity for the quantity of processes included.
5. Conclusion
In this work a detailed study on the + () reactive cross section has been done using a quasi-classical treatment, for collision energies from 1 meV up to 10 eV and for several vibrational states of the reactants and several isotopic variations. To this aim, new potential energy surfaces have been developed, one to include non-adiabatic transitions (PESTRIM8×8) and another to increase the accuracy in the whole energy range up to 17 eV (PES-NN). In all cases, it is found that from 1 meV to ≈ 0.5-1 eV the cross section behaves according to a Langevin law for charge induced-dipole long range interactions. For energies above 1 eV, the simulated cross section decreases fast, below the Langevin limit, for all initial vibrational states of (). However, for E > 1 eV, the reactive cross section exhibits a considerable increase with increasing , and the results seems to converge at =6.
It is found that the reactive cross section for =6, summing the and channels, match very well the recent experimental measurements by Savic et al. [Citation65], and also previous measurements [Citation54, Citation62], describing a broad energy interval from 0.5 meV to 10 eV. Moreover, the fact that for collision energies above 1 eV the measured reactive cross section show different values in different experiments, can be explained by a different vibrational excitation of achieved in each experimental setup.
Experimentally, the products are measured [Citation65]. The fact that in the QCT simulations, the cross section for the double fragmentation channel, , needs to be considered to get an agreement with the new experimental results can be explained by a classical artifact, since QCT method usually overstimate the energy transfer, specially dealing with systems with equal masses.
The new reactive cross sections obtained for + () and D + D have been included in a plasma model together with the widely used cross section [Citation98]. The resulting population densities of and D are proportional to the rate coefficient, which in turn indicates the sensitivity of the population density to the adopted cross section. The new cross sections for vibrational sates () of will be useful for the state-resolved CR modelling in plasma with higher molecular temperature, where higher collision energies (above 0.7 eV) become significant and the differences in the reactivity of the vibrationally excited will become important.
Acknowledgements
We thank Profs. S. Schlemmer and I. Savić for providing us the experimental measurement data and for very fruitfull discussions. Also, we acknowledge Prof. I. Tanarro for a carefull reading of the manuscript and very valuable comments and discussions, and Prof. F. Merkt for providing their experimental data. We acknowledge computing time at Cibeles (CCC-UAM) and Trueno (CSIC).
Disclosure statement
No potential conflict of interest was reported by the author(s).
Additional information
Funding
References
- T. Oka, Chem. Rev. 113 (12), 8738–8761 (2013). doi:10.1021/cr400266w
- S.C.O. Glover, Astrophys. J. 584 (1), 331–338 (2003). doi:10.1086/apj.2003.584.issue-1
- E. Herbst and T.J. Millar, in Low Temperatures and Cold Molecules, edited by Ian W.M. Smith (Imperial College Press, London, 2008), p. 1.
- V. Wakelam, E. Bron, S. Cazaux, F. Dulieu, C. Gry, P. Guillard, E. Habart, L. Hornekær, S. Morisset, G. Nyman, V. Pirronello, S.D. Price, V. Valdivia, G. Vidali and N. Watanabe, Mol. Astrophys. 9, 1–36 (2017). doi:10.1016/j.molap.2017.11.001
- T.J. Millar, A. Bernett and E. Herbst, Astrophys. J. 340, 906 (1989). doi:10.1086/167444
- L. Pagani, M. Salez and P. Wannier, Astron. Astrophys. 258, 479 (1992).
- J. Tennyson, Rep. Prog. Phys. 58 (4), 421–476 (1995). doi:10.1088/0034-4885/58/4/002
- B.J. McCall and T. Oka, Science 287 (5460), 1941–1942 (2000). doi:10.1126/science.287.5460.1941
- T. Oka, Phil. Trans. R. Soc. A. 370 (1978), 4991–5000 (2012). doi:10.1098/rsta.2012.0243
- W.D. Watson, Astrophys. J. 183, L17 (1973). doi:10.1086/181242
- E. Herbst and W. Klemperer, Astrophys. J. 185, 505 (1973). doi:10.1086/152436
- M. Quack, Mol. Phys. 34 (2), 477–504 (1977). doi:10.1080/00268977700101861
- D. Uy, M. Cordonnier and T. Oka, Phys. Rev. Lett. 78 (20), 3844–3847 (1997). doi:10.1103/PhysRevLett.78.3844
- T. Oka, J. Mol. Spectros. 228 (2), 635–639 (2004). doi:10.1016/j.jms.2004.08.015
- D. Gerlich, F. Windisch, P. Hlavenka, R. Plasil and J. Glosik, Phil. Trans. R. Soc. A. 364 (1848), 3007–3034 (2006). doi:10.1098/rsta.2006.1865
- K. Park and J.C. Light, J. Chem. Phys. 127 (22), 224101 (2007). doi:10.1063/1.2805394
- E. Hugo, O. Asvany and S. Schlemmer, J. Chem. Phys. 130 (16), 164302 (2009). doi:10.1063/1.3089422
- K.N. Crabtree, B.A. Tom and B.J. McCall, J. Chem. Phys. 134 (19), 194310 (2011). doi:10.1063/1.3587245
- D.S.N.G. Adams and E. Alge, Astrophys. J. 263, 123 (1982). doi:10.1086/160487
- K. Giles, N.G. Adams and D. Smith, J. Phys. Chem. 96 (19), 7645–7650 (1992). doi:10.1021/j100198a030
- M. Cordonnier, D. Uy, R.M. Dickson, K. Kerr, Y. Zhang and T. Oka, J. Chem. Phys. 113 (8), 3181–3193 (2000). doi:10.1063/1.1285852
- D. Gerlich and S. Schlemmer, Planet. Space Sci. 50 (12–13), 1287–1297 (2002). doi:10.1016/S0032-0633(02)00095-8
- D. Gerlich, E. Herbst and E. Roueff, Planet. Space Sci. 50 (12–13), 1275–1285 (2002). doi:10.1016/S0032-0633(02)00094-6
- K.N. Crabtree, C.A. Kauffman, B.A. Tom, E. Becka, B.A. McGuire and B.J. McCall, J. Chem. Phys. 134 (19), 194311 (2011). doi:10.1063/1.3587246
- T.J. Millar, Planet. Space Sci. 50 (12-13), 1189–1195 (2002). doi:10.1016/S0032-0633(02)00082-X
- L. Pagani, C. Vastel, E. Hugo, V. Kokoouline, C.H. Greene, A. Bacmann, E. Bayer, C. Ceccarelli, R. Peng and S. Schlemmer, Astron. Astrophys. 494 (2), 623–636 (2009). doi:10.1051/0004-6361:200810587
- T.A. Bell, K. Willacy, T.G. Phillips, M. Allen and D. Lis, Astrophys. J. 731 (1), 48 (2011). doi:10.1088/0004-637X/731/1/48
- F. Fontani, A. Palau, P. Caselli, A. Sánchez-Monge, M.J. Butler, J.C. Tan, I. Jiménez-Serra, G. Busquet, S. Leurini and M. Audard, Astron. Astrophys. 529, L7 (2011). doi:10.1051/0004-6361/201116631
- D. Gerlich, J. Chem. Phys. 92 (4), 2377–2388 (1990). doi:10.1063/1.457980
- I. Méndez, F.J. Gordillo-Vázquez, V.J. Herrero and I. Tanarro, J. Phys. Chem. A. 110 (18), 6060–6066 (2006). doi:10.1021/jp057182+
- M. Jiménez-Redondo, E. Carrasco, V.J. Herrero and I. Tanarro, PCCP 13 (20), 9655 (2011). doi:10.1039/c1cp20426b
- I. Tanarro and V.J. Herrero, Plasma Sources Sci. Technol. 20 (2), 024006 (2011). doi:10.1088/0963-0252/20/2/024006
- A. Grill, Cold Plasma in Materials Fabrication (IEEE Press, New York, 1994).
- M.A. Lieberman and A.J. Lichtenberg, Principles of Plasma Discharges and Materials Processing (Wiley, New York, 2005).
- D.W. Savin, P.S. Krstic, Z. Haiman and P.C. Stancil, AstroPhys. J. 606 (2), L167–L170 (2004). doi:10.1086/421108
- C.M. Coppola, S. Longo, M. Capitelli, F. Palla and D. Galli, Astrophys. J. Suppl. Ser. 193 (1), 7 (2011). doi:10.1088/0067-0049/193/1/7
- N. Indriolo and B.J. McCall, Astrophys. J. 745 (1), 91 (2012). doi:10.1088/0004-637X/745/1/91
- C.M. Coppola, D. Galli, F. Palla, S. Longo and J. Chluba, Mon. Not. R. Astron. Soc. 434 (1), 114–122 (2013). doi:10.1093/mnras/stt1007
- K.B. Chai, D.H. Kwon and M. Lee, Plasma Phys. Control. Fusion. 63 (12), 125020 (2021). doi:10.1088/1361-6587/ac2eb1
- R. Perillo, R.C.G.R.A. Akkermans, I.G.J. Classen, S.Q. Korving and M.P. team, Phys. Plasmas 26 (10), 102502 (2019). doi:10.1063/1.5120180
- S. Iordanova, T. Paunska, A. Pashov and S. Lishev, J. Phys. D. Appl. Phys. 52 (1), 015202 (2019). doi:10.1088/1361-6463/aae4ab
- K. Sawada and M. Goto, Atoms 4 (4), 29 (2016). doi:10.3390/atoms4040029
- D. Wünderlich and U. Fantz, Atoms 4 (4), 26 (2016). doi:10.3390/atoms4040026
- K. Verhaegh, et al. Nucl. Fusion. 61 (10), 106014 (2021). doi:10.1088/1741-4326/ac1dc5
- C. GIese and W.B. Maier II, J. Chem. Phys. 39 (3), 739–748 (1963). doi:10.1063/1.1734318
- R.H. Naynaber and S.M. Trujillo, Phys. Rev. 167 (1), 63–66 (1968). doi:10.1103/PhysRev.167.63
- W.R. Gentry, D.J. McClure and C.H. Douglas, Rev. Sci. Instrum. 46 (4), 367–375 (1975). doi:10.1063/1.1134225
- L.T. Specht, K.D. Foster and E.E. Muschlitz, J. Chem. Phys. 63 (4), 1582–1585 (1975). doi:10.1063/1.431481
- J.R. Krenos, K.K. Lehmann, J.C. Tully, P.M. Hierl and G.P. Smith, Chem. Phys. 16 (1), 109–116 (1976). doi:10.1016/0301-0104(76)89028-3
- I. Koyano and K. Tanaka, J. Chem. Phys. 72 (9), 4858–4868 (1980). doi:10.1063/1.439824
- S.L. Anderson, F.A. Houle, D. Gerlich and Y.T. Lee, J. Chem. Phys. 75 (5), 2153–2162 (1981). doi:10.1063/1.442320
- J. Shao and C.Y. Ng, J. Chem. Phys. 84 (8), 4317–4326 (1986). doi:10.1063/1.450053
- J.E. Pollard, L.K. Johnson, D.A. Lichtin and R.B. Cohen, J. Chem. Phys. 95 (7), 4877–4893 (1991). doi:10.1063/1.461704
- T. Glenewinkel-Meyer and D. Gerlich, Isr. J. Chem. 37 (4), 343–352 (1997). doi:10.1002/ijch.v37.4
- C.W. Eaker and G.C. Schatz, J. Phys. Chem. 89 (12), 2612–2620 (1985). doi:10.1021/j100258a036
- J.K. Badenhoop, G.C. Schatz and C.W. Eaker, J. Chem. Phys. 87 (9), 5317–5324 (1987). doi:10.1063/1.453649
- J.R. Stine and J.T. Muckerman, J. Chem. Phys. 68 (1), 185 (1978). doi:10.1063/1.435481
- M. Baer and C.Y. Ng, J. Chem. Phys. 93 (11), 7787–7799 (1990). doi:10.1063/1.459359
- A. Alijah and A.J.C. Varandas, J. Chem. Phys. 129 (3), 034303 (2008). doi:10.1063/1.2953571
- C. Sanz-Sanz, A. Aguado, O. Roncero and F. Naumkin, J. Chem. Phys. 143 (23), 234303 (2015). doi:10.1063/1.4937138
- P. Allmendinger, J. Deiglmayr, K. Höveler, O. Schullian and F. Merkt, J. Chem. Phys. 145 (24), 244316 (2016). doi:10.1063/1.4972130
- P. Allmendinger, J. Deiglmayr, O. Schullian, K. Höveler, J.A. Agner, H. Schmutz and F. Merkt, ChemPhysChem 17 (22), 3596–3608 (2016). doi:10.1002/cphc.v17.22
- K. Höveler, J. Delglmayr and F. Merkt, Mol. Phys. 119 (17–18), e1954708 (2021). doi:10.1080/00268976.2021.1954708
- F. Merkt, K. Höveler and J. Deiglmayr, J. Phys. Chem. Lett. 13 (3), 864–871 (2022). doi:10.1021/acs.jpclett.1c03374
- I. Savić, S. Schlemmer and D. Gerlich, Chem. Phys. Chem. 21 (13), 1429–1435 (2020). doi:10.1002/cphc.v21.13
- J.C. Tully, J. Chem. Phys. 93 (2), 1061–1071 (1990). doi:10.1063/1.459170
- C. Sanz-Sanz, O. Roncero, M. Paniagua and A. Aguado, J. Chem. Phys. 139 (18), 184302 (2013). doi:10.1063/1.4827640
- F.O. Ellison, J. Am. Chem. Soc. 85 (22), 3540–3544 (1963). doi:10.1021/ja00905a002
- F.O. Ellison, N.T. Huff and J.C. Patel, J. Am. Chem. Soc. 85 (22), 3544–3547 (1963). doi:10.1021/ja00905a003
- J.C. Tully, Adv. Chem. Phys. 42, 63 (1980). doi:10.1002/9780470142615.ch2
- P.J. Kuntz, in Atom-Molecule Collision Theory: A Guide for Experimentalists, edited by R.B. Bernstein (Plenum Press, New York, 1979), p. 79.
- A. Aguado, P. Barragan, R. Prosmiti, G. Delgado-Barrio, P. Villarreal and O. Roncero, J. Chem. Phys. 133 (2), 024306 (2010). doi:10.1063/1.3454658
- A. Aguado and M. Paniagua, J. Chem. Phys. 96 (2), 1265–1275 (1992). doi:10.1063/1.462163
- A. Aguado, C. Suarez and M. Paniagua, J. Chem. Phys. 101 (5), 4004–4010 (1994). doi:10.1063/1.467518
- C. Tablero, A. Aguado and M. Paniagua, J. Chem. Phys. 110 (16), 7796–7801 (1999). doi:10.1063/1.478688
- T.H. Dunning Jr, J. Chem. Phys. 90 (2), 1007–1023 (1989). doi:10.1063/1.456153
- A. Aguado, O. Roncero and C. Sanz-Sanz, PCCP 23 (13), 7735–7747 (2021). doi:10.1039/D0CP04100A
- L. Velilla, M. Paniagua and A. Aguado, Int. J. Quantum. Chem. 111, 387 (2010). doi:10.1002/qua.22588
- K. Sorbie and J. Murrell, Mol. Phys. 29 (5), 1387–1407 (1975). doi:10.1080/00268977500101221
- J. Murrell, K. Sorbie and A. Varandas, Mol. Phys. 32 (5), 1359–1372 (1976). doi:10.1080/00268977600102741
- A. Zanchet, P. del Mazo, A. Aguado, O. Roncero, E. Jiménez, A. Canosa, M. Agúndez and J. Cernicharo, PCCP 20 (8), 5415–5426 (2018). doi:10.1039/C7CP05307J
- O. Roncero, A. Zanchet and A. Aguado, Phys. Chem. Chem. Phys. 20 (40), 25951–25958 (2018). doi:10.1039/C8CP04970J
- J. Li, Z. Varga, D.G. Truhlar and H. Guo, J. Chem. Theory. Comput. 16 (8), 4822–4832 (2020). doi:10.1021/acs.jctc.0c00430
- B. Jiang and H. Guo, J. Chem. Phys. 139 (5), 054112 ( (2013). doi:10.1063/1.4817187.
- P. del Mazo-Sevillano, A. Aguado and O. Roncero, J. Chem. Phys. 154 (9), 094305 (2021). doi:10.1063/5.0044009
- A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai and S. Chintala, in Advances in Neural Information Processing Systems 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d' Alché-Buc, E. Fox and R. Garnett (Curran Associates, Inc., Vancouver, Canada, 2019), p. 8024–8035.
- L. Velilla, B. Lepetit, A. Aguado, J. Beswick and M. Paniagua, J. Chem. Phys. 129 (8), 084307 (2008). doi:10.1063/1.2973629
- A. Aguado, O. Roncero, C. Tablero, C. Sanz and M. Paniagua, J. Chem. Phys. 112 (3), 1240–1254 (2000). doi:10.1063/1.480539
- A. Zanchet, O. Roncero and N. Bulut, Phys. Chem. Chem. Phys. 18 (16), 11391–11400 (2016). doi:10.1039/C6CP00604C
- A.J. Ocaña, E. Jiménez, B. Ballesteros, A. Canosa, M. Antiñolo, J. Albadalejo, M. Agúndez, J. Cernicharo, A. Zanchet, P. del Mazo, O. Roncero and A. Aguado, Astrophys. J. 850 (1), 28 (2017). doi:10.3847/1538-4357/aa93d9
- M. Karplus, R.N. Porter and R.D. Sharma, J. Chem. Phys. 43 (9), 3259–3287 (1965). doi:10.1063/1.1697301
- T.P. Grozdanov and E.A. Solov'ev, J. Phys. B. 15 (8), 1195–1204 (1982). doi:10.1088/0022-3700/15/8/012
- C. Qu and J.M. Bowman, J. Phys. Chem. A. 120 (27), 4988–4993 (2016). doi:10.1021/acs.jpca.5b12701
- T. Nagy and G. Lendvay, J. Phys. Chem. Lett. 8 (18), 4621–4626 (2017). doi:10.1021/acs.jpclett.7b01838
- R.D. Levine and R.B. Bernstein, Molecular Reaction Dynamics and Chemical Reactivity (Oxford University Press, Oxford, 1987).
- C. Sanz-Sanz, A. Aguado and O. Roncero, J. Chem. Phys. 154 (10), 104104 (2021). doi:10.1063/5.0044320
- O. Roncero, V. Andrianarijaona, A. Aguado and C. Sanz-Sanz, Mol. Phys. 120 (1-2), e1948125 (2022). doi:10.1080/00268976.2021.1948125
- R.K. Janev, D. Reiter and U. Samm, Forschungszentrum Jülich Report JUEL-4105, 2003.
- I.I. Mizus, A. Alijah, N.F. Zobov, L. Lodi, A.A. Kyuberis, S.N. Yurchenko, J. Tennyson and O.L. Polyansky, MNRAS 468 (2), 1717–1725 (2017). doi:10.1093/mnras/stx502
- I. Koleva, T. Paunska, H. Cshlüter, A. Shirarova and K. Tarnev, Plasma Sources Sci. Technol. 12 (4), 597–607 (2003). doi:10.1088/0963-0252/12/4/311
- R.J. Trudeau, Introduction to Graph Theory (Dover Publications, Dover, 1993).
- S. Even, Graph Algorithms (Cambridge University Press, Cambridge, 2011).
- F. Harary and E.M. Palmer, Graphical Enumeration (Elsevier, Academic Press, 2014).
- A.P. Mani and R.J. Stones, SIAM J. Discrete Math. 30 (2), 1046–1057 (2016). doi:10.1137/15M1024615
- L. St-Onge and M. Moisan, Plasma Chem. Plasma Process. 14 (2), 87–116 (1994). doi:10.1007/BF01465741
- V.E. Golant, A.P. Zhilinski and I.E. Sakharov, Fundamentals of Plasma Physics (Wiley, New York, 1977).
- E.W. McDaniel and E.A. Mason, The Mobility and Diffusion of Ions and Gases (Wiley, New York, 1973).
- J.P. Boeuf, G.J.M. Hagelaar, P. Sarrailh, G. Fubiani and N. Kohen, Plasma Sources Sci. Technol. 20 (1), 015002 (2011). doi:10.1088/0963-0252/20/1/015002
- K.B. Chai and D.H. Kwon, J. Quant. Spectrosc. Radiat. Transfer. 227, 136–144 (2019). doi:10.1016/j.jqsrt.2019.02.015
- W.L. Wiese and J.R. Fuhr, J. Phys. Chem. Ref. Data 38 (3), 565–720 (2009). doi:10.1063/1.3077727
- R.K. Janev, W.D. Langer, K.J. Evans and D.E.J. Post, Elementary Processes in Hydrogen-Helium Plasmas (Springer, Berlin, 1987).
- G.J.M. Hagelaar, G. Fubiani and J.P. Boeuf, Plasma Sources Sci. Technol. 20 (1), 015001 (2011). doi:10.1088/0963-0252/20/1/015001
- W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes in C: the art of scientific computing (Cambridge University Press, New York, 1992).
Appendix A.
N-body permutational invariant polynomials
In this section we provide a definition of a permutational invariant polynomial in terms of graph theory, which we latter use as a way of filtering the N-body polynomials from a general set of PIP.
A graph is a collection of vertices and edges () [Citation101]. In an undirected graph the edges are non-ordered pairs of vertices. An undirected graph is said to be connected if every pair of vertices are joined by a path.
A relation can be established between the polynomials in Equation (Equation5(5) (5) ) and an undirected graph, where the vertices represent the particles and the edges the monomials between them.
For instance, the following polynomial (A1) (A1) can be expressed as the following graph:
This corresponds to a disconnected graph since there are various pairs of vertices which are not reachable, for instance from vertex 2–3. An example of a polynomial whose graph is connected is: (A2) (A2)
A polynomial P is said to be a N-body polynomial if its corresponding undirected graph is connected. Any other polynomial that arises as with being any operation of a permutation group will be a N-body polynomial if P is. Note that the effect of a permutation operation in the graph only affects the order of the vertices and relabel of the edges: (A3) (A3)
Given a graph, there are simple algorithms as Depth-First Search (DFS) [Citation102] to compute whether it is a connected graph or not, which recursively traverses the graph marking each visited vertice. If at the end of the execution all nodes were visited, the graph is connected. There exist a finite number of connected N-vertices undirected graphs which can be evaluated using the recursive formula [Citation103]: (A4) (A4) These numbers are tabulated for N up to 16 in [Citation104]. One should note that the number of connected undirected graphs increases fast, for instance for a set of six vertices there exist 26704 of those, so in practice an upper limit on the number of edges has to be set.
At this point we have the tools to determine the number of N-body polynomials, as well as, given a polynomial, to determine if it is N-body. If the system presents some kind of permutational symmetry, a minimal set of N-body PIP will be generated by projecting the above polynomials onto the totally symmetric irreducible representation of the permutation group. The dimension of the minimal PIP set is necessarily lower or equal to the minimal one of polynomials, as many of them are related by permutation operations. Note that we never mentioned the exponents l of the monomials, since they play no role in the graph construction. Hence, what we have defined up to here is a generator of N-body polynomials or PIPs. Following the general procedure, we are now free to set the desired maximum polynomial degree and produce the combinations of monomial exponents l to generate a finite set of functions.
Appendix B.
Method for the plasma model
The method used to model the plasma is similar to that previously described [Citation41, Citation100], in which each particle's level i density can be solved from a continuity equation (A5) (A5) where D denotes diffusion coefficient and the right hand side is the net particle source and sink term by collisional-radiative (CR) process and particle flux into and out a volume. In steady state of , the density balance equations for a long cylindrical vessel plasma can be expressed as follows. For atomic levels of (A6) (A6) and for species , , and (A7) (A7) The density of , is deduced from the quasi neutrality condition (A8) (A8) The electron density, , is assumed to have a radial distribution close to a Bessel-type profile for the ambipolar-diffusion regime considered here and applying the Bohm criterion as boundary conditions between the plasma and the vessel walls with μ = 2.405 for an infinite cylinder of the effective radius of R [Citation34, Citation100]. The effective R is set as 40 cm for our plasma device.
Diffusion time τ for H atom in the device of radius is given by with the thermal velocity for atomic temperature and the hydrogen mass . The wall recombination coefficient γ is given as the empirical expression [Citation105], with being the temperature of molecule.
Under the ambipolar diffusion assumption [Citation106] the diffusion coefficient for is given by (A9) (A9) where is in Torr and the electronic temperature is in eV. The reduced mobility for and D are 15.9 and 11.2, respectively [Citation107]. The diffusion coefficients for and are given with the relation [Citation105].
The conversion coefficients from H ions into H atom () and molecule () on the wall are taken from [Citation108]. The conversion coefficient from H atom to molecule on the wall is also taken from [Citation108]. , Q and V denote the gas flow rate of in the inlet, the pumping rate in the outlet and the volume of the plasma device, respectively. , Q and V are given by 600 sccm (standard cubic centimeters per minute), 4800 lps (liter per second) and 2.64×10 cm, respectively.
The CR processes and the rate coefficient α notations are listed in Table with the references of the collision cross section and the radiative transition rate. The escaping factor η for radiation trapping in optically thick plasma is given as that for the infinite long cylinder [Citation109].
From the cross section the rate coefficient is obtained as (A10) (A10) as the averaging over the relative velocity distribution. When the colliding particles of the masses and have the Maxwellian energy distribution with temperatures and the rate coefficient α can be expressed as [Citation112] (A11) (A11) where for temperatures in eV. For electron collisions and , with being the electron mass. The Maxwellian rate coefficient for the heavy particle collision + (n = 1) is obtained with and .
The balance equations of including the nonlinear terms of , and are solved by the multidimensional secant Broyden's method [Citation113] setting the initial as the solution of the linear part of Equations (EquationA6(A6) (A6) ) and (EquationA7(A7) (A7) ) for various , , and .