Rafael A.
Jara-Toro
a,
Octavio
Roncero
b and
François
Lique
*a
aUniv Rennes, CNRS, IPR (Institut de Physique de Rennes), UMR 6251, F-35000 Rennes, France. E-mail: francois.lique@univ-rennes.fr
bInstituto de Física Fundamental (IFF-CSIC), C.S.I.C., Serrano 123, 28006 Madrid, Spain
First published on 23rd July 2024
CH+ was the first molecular ion identified in the interstellar medium and is found to be ubiquitous in interstellar clouds. However, its formation and destruction paths are not well understood, especially at low temperatures. A new theoretical approach based on the canonical variational transition state theory was used to study the H + CH+ reactive collisions. Rate coefficients for formation of C+ ions are calculated as a function of temperature. We considered the participation of a direct path and an indirect path in which the reactants should overcome an entropic barrier to form a van der Waals complex or pass through a CH2+ intermediate complex, respectively. We show that the contribution of both pathways to the formation of C+ has to be taken into account. The new reactive rate coefficients for the title reaction, complemented by reactive data for CH+/CH2+ in the H/H2/He mixture, have been used to simulate the corresponding kinetics experimentally measured using an Atomic Beam 22 Pole Trap apparatus at low temperature. A good agreement with the experimental findings was found at 50 K. At a lower temperature, the model overestimates the formation of C+. This shows that secondary reactions are not responsible for the weak C+ production in the experiments at such temperature. Then, we discuss the possible impact of non-adiabatic effects in the study of the H + CH+ reactive collisions and we found that such effects can be responsible for the decrease of the H + CH+ rate coefficients at low temperature. This study offers an explanation for the disagreement between H + CH+ theoretical and experimental rate coefficients which has been going on for 20 years and highlights the need for performing non-adiabatic studies for this simple chemical reaction.
In this context, the destruction of CH+ due to H collisions, a key destruction mechanism of CH+ in the ISM, has received considerable attention on both the experimental and theoretical sides.19,21–32 Despite such high interest, theoretical works failed at reproducing the experimental rate coefficients at low temperatures (below ∼60 K). The experimental studies carried out by Federer et al. (1984, 1985),21,22 Luca et al. (2005),24 and Plasil et al. (2011)25 found similar rate coefficients at temperatures above 60 K. All these experimental studies found a reactive rate coefficient of the order of 1.2 ± 0.5 × 10−9 cm3 s−1, in agreement with the most accurate quantum calculations and in accordance with Langevin capture theory. However, all the theoretical studies, including quantum time independent studies29 normally well suited for such studies, failed at reproducing the experimental measurement at low temperatures.19,26–31,33 Surprisingly, the experimental findings suggest a steep fall off of the reactive rate coefficients with decreasing temperature, even reaching a difference of 2 orders of magnitude with room temperature measurements, in contradiction with what is expected from simple capture models.
The most recent experimental results have been obtained using an Atomic Beam 22 Pole Trap (AB-22PT) apparatus and cover the 10–100 K temperature range.25,34 In the experimental device, H atoms were produced by a radio frequency (RF) discharge from H2 molecules so that both H and H2 coexist and influence the kinetics in the device. Due to the competition between reactions of CH+ with H atoms, H2 molecules and secondary reactions leading finally to CH3+ in the experimental setup, it was found to be necessary to take into account not only the diminution of the CH+ concentration but also the subsequent increase of the C+ concentration and the effective number density of H to obtain a reliable estimation of the rate coefficients for the title reaction.
The purpose of this study is to employ a new theoretical approach that allows understanding the temperature variation of the experimental rate coefficients in the 10–1000 K range and to have a better insight into the H + CH+ reaction mechanism. In particular, we found that correctly modeling the kinetics in the experimental device is important for a full understanding of the experimental measurements at low temperatures. Additionally, we found that non-adiabatic effects can be crucial for the theoretical study of the H + CH+ reaction. Indeed, by carefully modeling the kinetics of CH+ in the trap and by considering the possible impact of non-adiabatic effects, we were able to understand the fall off of the H + CH+ reaction rate coefficients and bring an explanation for the disagreement between previous theoretical and experimental studies of the H + CH+ reaction at low temperatures.
The capture rate coefficients for the CH+ + H reaction path were calculated using the canonical variational transition state theory (CVT).38 This approach was chosen because the reaction bottleneck occurs at a configuration on these paths where the free energy reaches a maximum at every temperature considered in this study. For both reaction paths, we used the Rice–Ramsperger–Kassel–Marcus (RRKM) formalism implemented in the Mesmer program,39 to obtain microcanonical rate coefficients for dissociation of the intermediate complexes as well as for the crossing of the reaction barrier. The canonical rate coefficients for these processes were calculated by multiplying the Maxwell–Boltzmann distribution at a given temperature with the microcanonical rate coefficients at approximately the energy of the reactants for the nascent CH2+* and CH+⋯H intermediates on both paths (a loss of 90 cm−1 was considered due to the conversion of internal energy to kinetic energy). A tunneling effect was calculated using the Eckhart method.40
In parallel, quasi-classical trajectory (QCT) calculations have been performed with the classical molecular dynamics with the quantum transition (MDwQT) program41,42 on the CH2+ global potential energy surface of Stoecklin and Halvick.26 See Section S3 in the ESI† for a more detailed description.
For the direct path, we found that, to reach the van der Waals pre-reactive complex (CH+⋯H), the H has to cross positive values of the free energy in the short range. This activation barrier is induced by the repulsive force due to the collinear approach of both hydrogens and this activation barrier is leading to a decrease of the rate coefficients with decreasing temperatures. For the indirect path, for temperatures below 100 K, the free energy is always below the energy of the reactants so that the capture rate coefficients follow the Langevin capture rate coefficients.
In Fig. 2, the thermal rate coefficients for the two direct (kd) and indirect (ki) processes are plotted as a function of temperature. In Fig. 2, we also include the global calculated rate coefficients (kg) resulting from the sum of both processes. As one can see that our thermal rate coefficients resulting from the sum of both processes are in good agreement with the data obtained from accurate quantum scattering methods,29 except at very low temperature where our data are lower than the former one, the deviation being up to a factor ∼2 at 10 K. The deviation between quantum and present rate coefficients can be explained by the relatively simple approach used in this work, but also by the high difficulty of converging pure quantum reactive data for such a complex reaction, especially at low temperatures. The good overall agreement between present and state of the art calculations is however demonstrating the relative accuracy of our modeling. For the indirect process that implies the formation of a stable intermediate complex, we found that the rate coefficients are quasi-independent of the temperature up to 300 K. In contrast, for the direct path, the activation barrier generated by repulsive forces is significant at low temperatures and induces the steep fall off on the rate coefficients with a decreasing temperature.
Fig. 2 Calculated and experimental rate coefficients for the title reaction. The blue‡ and red lines correspond to the present theoretical results for both indirect and direct paths, respectively, while the black line represents the global (sum of both processes) rate coefficients. Purple and gold symbols correspond to our QCT calculations. Experimental data (Plasil et al.25 (red dots), Luca et al.24 (green triangles), Gerlich et al.34 (green dot), and Federer et al.21,22 (open squares)) as well as the most accurate theoretical results of Werfelli et al.29 (grey dashed line) are presented. |
The QCT results support the hypothesis of two competitive reactive paths, provided from the present microscopic interpretation. Indeed, the QCT rate coefficients show a near Langevin behavior and are in near quantitative agreement with the quantum results30 presented in Fig. 2 for temperatures above 100 K. The thermal rate coefficients were calculated under such conditions (see purple values in Fig. 2 and Fig. S3, ESI†) and are in good agreement with the Langevin limiting value and with the exact quantum results.29
In contrast, the analysis of the trajectories shows that below 100 K, the impact parameter increases a lot and the duration of the trajectories increases enormously (see Section S3 in the ESI†). Thus, if we consider only trajectories ending in 7 picoseconds (ps) which could be associated with the direct mechanism, the QCT rate coefficients show a decreasing behavior with decreasing temperature in good agreement with the present rate coefficients for the direct process (kd). In this view, the remaining set of trajectories, for more than 7 ps and associated with the indirect process, could give an opportunity for secondary collisions (as those described in the next section) to compete with the formation of either H2 + C+ or H + CH+.
Fig. 3 Chemical reactions possibly occurring in the experimental setup with the RF discharge set to OFF and ON. |
If only H2 molecules are present in the experimental setup (e.g. the RF discharge is OFF), CH+ reacts quickly with them so that the formation of CH2+ (R1), possibly followed by a secondary reaction leading to the CH3+ product (R2), could occur. Reactions leading to more complex products such as the CH3+ + H2 reaction were not considered because of the absence of CH4+ or CH5+ in the experiment.25,34,43–45
When the RF discharge is ON (e.g. the discharge dissociates H2 into H), the previously mentioned reactions are occurring since H2 is still present in the device but, in addition to R1 and R2, R3 to R8 are likely to influence the kinetics of CH+ and as a consequence the ionic concentrations are observed. R3 to R7 are connected to the indirect path for the H + CH+ reaction while R8 corresponds to the H + CH+ direct path.
The first step for the indirect path is the formation of the CH2+* intermediate complex with an excess of internal energy (R3). This complex can (i) be destroyed to form C+ (R4), (ii) lead to CH3+ formation via secondary (R5, R6) reactions and (iii) be stabilized via collisions with the buffer gas (He) used to cool the stored ions in the system (R7). The formation of CH2+ in the experiment is then not only due to the reaction between CH+ and H2 (R1), but also due to the reaction between CH+ and H (R3 + R7).
When examining the experimentally measured concentration at temperatures of 12 and 50 K, it is obvious that the concentration of CH2+ significantly increases when the RF discharge switches from OFF to ON. However, this increase does not seem to be dependent on the temperature. This finding combined with the fact that the relative concentration of C+ is significantly dependent on the temperature (see Table 1) indicates that the formation of the stable CH2+ intermediate complex is not necessarily leading to the formation of C+ + H2 products. Then, it seems obvious that one has to model the full kinetics in the experimental device if one wants to accurately extract the H + CH+ rate coefficients.
Relative concentrations | |||
---|---|---|---|
12 K | 50 K (RF-ON) | 50 K (RF-OFF) | |
The concentrations are relative to the initial CH+ concentration. a, b, and c are the relative concentration values at 50 ms, 100 ms, and 150 ms, respectively. Values taken from R. Plasil, et al. (2011)25 and Luca et al. (2006).24 | |||
CH2+ | 0.13a | 0.14a | 0.03a |
0.24b | 0.21b | 0.07b | |
0.32c | 0.30c | 0.09c | |
C+ | 3 × 10−4a | 0.03a | — |
4 × 10−4b | 0.06b | — | |
6 × 10−4c | 0.08c | — |
We solved the kinetic equations with the Kintecus’ program46 in order to simulate the experimental concentrations of C+, CH+, CH2+ and CH3+ resulting from all the possible reactions in the experimental setup (see the ESI† for more details about the simulation).
For the rate coefficients corresponding to R1, R2, R5 and R6, we used the data provided by Plasil et al.25 and Luca et al.24 For R7, we used the collisional limit of the rigid sphere theory. More information about the values used for R1 to R8 is found in Table 2.
Rate coefficients | |
---|---|
k R1, kR2 and kR5 were taken from ref. 1. A range of rate coefficients was used for kR6 without any effect on the global concentrations. For kR7, the collisional limit for the rigid sphere theory was used. | |
k R1 | 1.2 × 10−9 cm3 s−1 |
k R2 and kR5 | 1.6 × 10−9 cm3 s−1 |
k R6 | 1.0 × 10−9–1.0 × 10−11 cm3 s−1 |
k R7 | 3 × 10−10 cm3 s−1 |
First, we investigate the experimental conditions at 50 K with the RF discharge set on OFF (only R1 and R2 occur).24 As shown in Fig. 4a, it is possible to reproduce very well the experimental concentrations as a function of time considering only the CH+ + H2 → CH2+ + H and CH2+ + H2 → CH3+ + H reactions.
Fig. 4 Ionic concentration as a function of the storage time. Upper panels correspond to the experimental results of Gerlich et al. (2001)34 and Plasil et al. (2011).25 Lower panels are the results from our kinetic model based on the reaction scheme of Fig. 1. Dotted and solid lines correspond to the model using set 1 and set 2 of data, respectively. |
Then, we simulate the concentrations of C+, CH+, CH2+ and CH3+ at 50 and 12 K but when the RF discharge is set to ON (both H and H2 are present in the device), two simulations were performed. In the first one, R1, R2 and only the formation of C+ through the H + CH+ reaction using kg rate coefficients were used (set 1, dotted line). In the second one, R1 to R8 were included (set 2, solid line).
We compare the results of our simulation (using sets 1 and 2 of the data) with the results at 50 K reported by Gerlich et al. (2001)34 and Plasil et al. (2011)26 as shown in Fig. 4b. The results of the kinetic simulation using datasets 1 and 2 are in good agreement demonstrating the moderate impact of secondary events in the kinetics (the main difference between datasets 1 and 2 being the inclusion of secondary events). The agreement between the calculated and experimental C+, CH+, CH2+ and CH3+ concentrations is very good. This demonstrates that the value of the overall calculated rate coefficients for the H + CH+ reaction is of good accuracy.
Fig. 4c shows the results of our simulation at 12 K. The global agreement between the calculated and experimental CH+, CH2+ and CH3+ concentrations is good for the two datasets. However, the C+ concentration is overestimated by a factor 3–4 in the simulation performed with the two datasets in the kinetic model. Secondary reactions, despite being slightly more important than at 50 K, cannot explain the important decrease of the C+ concentration in the experimental device. Despite the deviation not being dramatic, our model fails at reproducing the experimental C+ concentration at 12 K and shows that secondary reactions cannot explain the deviation between theoretical and experimental results.
These simulations show that the impact of secondary reactions, despite not fully negligible, cannot explain alone the deviations between the calculated and measured H + CH+ reaction rate coefficients at low temperatures. Indeed, considering the lifetime of the intermediate complex, one would need a He density in the experimental setup higher by at least 3 orders of magnitude than what was considered in order to have a secondary event important. Then, one has to investigate the impact of non-adiabatic effects, neglected in the present study, if one wants to have a better insight into the full H + CH+ reaction mechanism. Such effects can be of high importance, especially at low temperatures, and may impact the efficiency of the title reaction.
In the left panels of Fig. 5, the adiabatic energies of the two lower 2A′ and the lower 2A′′ electronic states of CH2+ at three geometries are presented. The calculations have been performed with the MOLPRO suite of programs46 using an explicitly correlated internally contracted multi-reference configuration interaction (ic-MRCI-F12)48–50 and the cc-pCVTZ-F12 electronic basis set.51 The molecular orbitals are optimized using the state-averaged complete active space self-consistent field (SA-CASSCF) method, with 6 active orbitals and considering 6 electronic states (3A′ and 3A′′).
For the H–C–H configuration (top panels of Fig. 5), there is a clear crossing between the two A′ curves, which gives rise to a shallow barrier at this collinear geometry. As discussed by Stoecklin and Halvick,26 CH+ in its ground electronic state (X1Σ+) has no electron available for bonding and its interaction with the hydrogen atom gives rise to a repulsive curve. However, the CH+ in its first excited electronic state (A1Π) forms a deep well at this geometry. The curves of these two states cross at a CH distance of approximately 2.5 Å, giving rise to the ground electronic state potential with a small barrier.
The top of this barrier at collinear geometry corresponds to a conical intersection (CI) and clearly suggests the possibility of non-adiabatic transitions towards the excited 2A′ state. To illustrate this further, diabatization has been performed using the method described in ref. 52. The diabatic curves are shown in the right panels of Fig. 5. Using such a method, an electronic subspace of three states is considered, formed by doubly degenerate Π electronic states and one Σ electronic state, corresponding to E1 and E0 in the figure. The Σ electronic state is coupled to the two Π electronic states, by the coupling V, but there is no direct coupling between the two Π electronic states. As expected, the E0 is repulsive while the two E1 form a deep well and there is a non-negligible coupling among them (it should be noticed that the calculations are done at an angle of 179.99 degrees, and this is the reason for the non-zero coupling). A very similar situation holds for other angles where the crossings are avoided because the couplings are much stronger.
The barrier originated by this CI was included in the fits made for the ground adiabatic electronic state,26,28 and accurate quantum calculations performed on these PESs29–31 did not show any decrease in the H + CH+ → H2 + C+ reactive rate coefficients. Such behaviour is also supported by the present QCT calculations (see the ESI†), which show a near Langevin behaviour as the quantum calculations,29 and we may then conclude that the cusp introduced by the CI in the ground electronic state is not able to reduce the reactivity.
However, near conical intersections, the flux can bifurcate among the different electronic states involved in the crossing, and we therefore expect some reduction in the reactive rate coefficients when including electronic transitions. It is also interesting to note that the non-adiabatic diagonal terms, proportional to the second derivatives of the electronic functions with respect to nuclear coordinates, typically include an additional barrier, which could yield a partial reflection back, thus reducing the reaction rate coefficients.
The complex non-adiabatic character appearing in the CH2+ complex is further illustrated in the middle panels of Fig. 5, where cuts of the first six electronic states (32A′ and 32A′′) are shown for a T-shaped geometry (in fact 89 degrees to keep the Cs symmetry), for a slightly elongated H2 (at r = 1.51 Bohr, r being the H–H distance). The upper three electronic states correlate with the C(3P) + H2+ (X2Σ) asymptote, leading to two states out of the plane and a state in the plane of the molecule. These three states show a slight energy increase of energy between 4 and 6 Bohrs, and at shorter distances they show a sudden energy decrease, due to crossings with other even more excited Rydberg states. As the H2 internuclear distance increases and the C+ insets in the middle, the well appearing in the 22A′ state, at R = 2 Bohr (R is the distance between C+ and the center of the mass of H2) and 1 eV, gets deeper, being responsible for the deep well of the CH2+ ground electronic state. Another well is also apparent at this geometry for the 32A′′ state, which give rises to the second Π electronic state originating from the Renner–Teller effect in CH2+.
This situation was also described in ref. 47, and the CI appears in this case for T-shaped geometries. The crossings are also illustrated in the diabatic representation shown in the right-middle panel of Fig. 5. All these show that below the H + CH+ asymptote there are several electronic states which are strongly coupled. Thus, strong non-adiabatic transitions are expected to occur, which may reduce the reaction probability.
In the H–H–C geometry shown in the top panels of Fig. 5, in the exit channel towards the H2 (X1Σ+) + C+(2P) asymptote, the lower eigenvalues correspond to the doubly degenerate Π electronic states, and a similar situation holds for other bent geometries.
All these results indicate the possibility that non-adiabatic transitions among these electronic states are indeed possible and should have an effect on the dynamical calculations. To show this, we have performed on-the-fly classical trajectories using Tully's method and a lower level of theory, a time-dependent density functional theory (TD-DFT). 100 non-adiabatic dynamics were made (see Fig. 6) starting from the CH2+* molecule (the CH2+* molecule has an energy corresponding to that of the reactants). The non-adiabatic dynamics were carried out with the Newton-X program interfaced with the GAUSSIAN (version 16 Rev. C.01) using the wb97xd/aug-cc-pvdz with the long-range and dispersion corrected level of theory.53 Such a method was selected because of its relative accuracy and especially its numerical efficiency since multiple simulations of the dynamics were required for high confidence in the results. The accuracy of this approach versus that of highly correlated quantum chemistry approaches was tested (see Section S4 in the ESI†) and we found that this level of theory can provide a good representation of the electronic energies and structure of the collisional complex.
Fig. 6 Average populations of the ground (1A′) and first excited (2A′) electronic states for CH2+* as a function of time. |
It can be seen that the population initially in the 12A′ state is transferred very quickly to the 22A′ state, where it remains for a long time. The increase of the density of vibronic states, because of the many electronic states participating in the dynamics, may induce the increase of the life-time of the CH2+* intermediate complex. Under such a situation, the long-lived CH2+* complex could have more time to suffer other relaxing processes (R5–R7 in Fig. 2), thus reducing the C+ channel. Also, the crossings among electronic states of the ionic character could also lead to the increase of diagonal and non-diagonal transition dipole moments which could accelerate to some extent the radiative emission to form CH2+.
Then, it seems crucial to perform a non-adiabatic study of the H + CH+ reaction involving multiple potential energy surfaces and couplings between them as emphasised by the authors of experimental studies on the title system. Indeed, the indirect path considered in most of the theoretical studies involves a Langevin reaction which only occurs under adiabatic conditions. The present study is the first one allowing qualitative and (in part) quantitative analyses of the CH+ + H experiments at low temperatures. It is indeed offering an alternative to the persisting disagreement between the experimental and theoretical studies for the title reaction. It is now obvious that a non-adiabatic study of the reaction has to be performed in order to improve our knowledge of the reaction and possibly for the theory to match the experimental data.
From the astrophysical point of view, the decrease of the rate coefficients at low temperatures is consistent with the relatively high abundance of CH+ in cold astrophysical media. Indeed, the use of much lower rate coefficients for the H + CH+ reaction than those actually used (∼10−9 cm3 s−1, see the KIDA database55) would help in reconciliating observations of CH+ and astrochemical models. We also highlight that new experimental studies at low temperature using an alternative experimental setup such as molecular beams would be very valuable in order to better know and understand the peculiar reaction mechanism of the H + CH+ reaction.
Footnotes |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp01902d |
‡ Dotted blue line is an extrapolation of the rate coefficients at low temperature due to the impossibility to converge the calculations. |
This journal is © the Owner Societies 2024 |