The C(³P) + CH₂(³B₁) reaction: accurate electronic structure calculations and kinetics for astrochemical modeling

Abstract

In this work, we present a theoretical study of the C + CH₂ → H + CCH reaction, which is expected to play a key role in the formation of carbon-chain species in molecular clouds. We employ high-level ab initio electronic structure methods, including the HEAT-345Q composite scheme and multi-reference configuration interaction, to characterize the potential energy surface. Reaction kinetics are investigated through both conventional transition state theory for tight transition states and variational reaction coordinate transition state theory (VRC-TST) for barrierless pathways. The master equation technique is used to account for competing complex formation, isomerization, and dissociation processes, particularly under collisionless conditions relevant to astrochemistry. Rate constants for product formation are calculated over a wide temperature range (10–300 K), and fitted to extended Arrhenius expressions. The results reveal that both singlet and triplet pathways contribute significantly to the overall reaction rate, with the triplet surface dominating under typical ISM conditions.

---

1 Introduction

To date, more than 300 molecules have been detected in the interstellar medium (ISM), and the number is continuously growing.^1,2 These molecules range from the simplest molecules that exist, such as methylidene (CH), the first molecule detected in the ISM in 1937,³ to complex molecules such as fullerenes,⁴ and cyanonaphtalene.⁵ The ISM consists of diverse regions, where molecular abundances vary depending on local conditions.

To model the chemistry happening in different regions of the ISM, reaction networks such as Nautilus,⁶ RATE22⁷ and UCLCHEM⁸ have been developed. These networks typically incorporate hundreds of species and thousands of reactions, occurring in both the gas phase and on dust grain surfaces. They are used to simulate the evolution of the abundances of different molecules over time and to compare them with observational data. In these models, rate constants for individual reactions are crucial, but due to the extreme conditions in the ISM, these rate constants are not easily measured. Indeed, in dense molecular clouds, temperatures can be as low as 10 K with pressures down to about 10⁻¹⁵ bar.⁹ These rate constants are tabulated in databases, such as the KInetic Database for Astrochemistry (KIDA).¹⁰ The lack of data under conditions similar to those in the ISM results in many of the rate constants having significant uncertainties, with values being extrapolated from similar reactions or experimental data from largely different conditions. Because these reaction networks encompass too many reactions to each be exhaustively studied, sensitivity analysis are used to determine the reactions that will have significant impacts on the abundances of the species of interest.^11,12

On the experimental side, the CRESU technique¹³ allows for the direct measurements of rate constants over a wide range of temperatures, but reactions between two radical species, such as atomic carbon (C) and methylene (CH₂), are difficult to study experimentally. Additionally, the CRESU technique¹³ is not always able to probe the products formed and can then not determine the branching ratios of the reaction. In these cases, computational chemistry is a precious tool for computing rate constants that will be included in reaction networks.

Among the numerous astrochemical reactions, the reaction C + CH₂ → C₂H + H has not been studied before but has been shown to be a key reaction in astrochemical networks.¹² Atomic carbon is expected to be one of the most abundant species in molecular clouds, especially in the early and intermediate stages of cloud evolution, with an abundance possibly reaching up to 10⁻⁴ with respect to H₂.^14–18 CH₂ has been detected in various molecular clouds^19–21 and is expected to be formed mainly through the CH₃⁺ + e⁻ → CH₂ + H reaction.²² It is also expected to play a role in the chemistry of Titan as a photodissociation product of methane.^23,24 The C + CH₂ reaction produces hydrogen atoms (H) and the ethynyl radical (C₂H), the latter being involved in the growth of polycyclic aromatic hydrocarbons (PAHs)²⁵ and was first detected in 1974.²⁶ Acetylene (HCCH), which is an intermediate in the reaction pathway, was detected in 1989²⁷ and is an important molecule in extraterrestrial chemistry.²⁸ Vinylidene (CCH₂), another intermediate, remains undetected,²⁹ in part due to the short lifetime of its singlet ground state.³⁰

In this work, we present high-accuracy electronic structure calculations for the C + CH₂ reaction. We employ the high-accuracy extrapolated ab initio thermochemistry (HEAT) scheme,³¹ along with the multi-reference configuration interaction (MRCI) method,³² to obtain accurate information regarding the stationary points on the potential energy surface (PES) of the system. These electronic structure calculations allow us to characterize intermediates and transition states, which are necessary to take into account when studying the reaction kinetics. To investigate the reaction kinetics, we apply a statistical rate theory based master equation method, allowing us to obtain rate constants for the C + CH₂ → H + C₂H reaction.

2 Methodology

In typical chemical kinetics calculations, it is necessary to obtain accurate information on the stationary points of the PES. These can be used to compute microcanonical rate constants that are then input in the master equation (ME) to obtain phenomenological rate constants. This section will describe each step of these processes.

2.1 Electronic structure calculations

Kinetic modeling requires highly accurate evaluation of reaction barriers and sometimes also reaction energies to obtain reliable results. Therefore, we employed the renowned HEAT scheme in its 345-Q form, hereafter referred to as HEAT345-Q.³¹ This scheme has been shown to give reaction energies with errors lower than 1 kJ mol⁻¹. The first step in HEAT345-Q is to perform a geometry optimization with coupled-cluster methods including single and double excitations and including triples in a perturbative way (CCSD(T)),^33,34 using the correlation-consistent cc-pVQZ basis set of Dunning et al.^35,36

From the optimized geometry, several energy contributions are included to obtain accurate energies at a reasonable computational cost. In HEAT345-Q, the total energy is defined as a sum of contributions calculated for the obtained geometry, as given in eqn (1):

E = E^∞_HF + ΔE^∞_CCSD(T) + ΔE_CCSDT + ΔE_CCSDTQ + ΔE_rel + ΔE_DBOC + ΔE_ZPE (1)

In eqn (1), E^∞_HF is the complete basis set extrapolated Hartree–Fock energy, ΔE^∞_CCSD(T) is the complete basis set extrapolation of the correlation energy obtained from a CCSD(T) calculation. Further electron correlation terms (ΔE_CCSDT and ΔE_CCSDTQ) are obtained through the use of frozen core calculations and smaller basis sets compared to ΔE^∞_CCSD(T). ΔE_rel includes relativistic corrections such as the Cowan–Griffin term and the two-electron Darwin energy.^37–39 ΔE_DBOC is the diagonal Born–Oppenheimer correction (DBOC), which is added to account for the Born–Oppenheimer approximation.⁴⁰ ΔE_ZPE is finally the ZPE term corrected for anharmonicity employing VPT2, obtained at the same level of theory as the geometry optimization (i.e. CCSD(T)/cc-pVQZ).⁴¹ A more thorough description of the HEAT scheme and each term involved in eqn (1), along with the basis set used in the calculation of each term can be found in the paper by Bomble et al.³¹

Geometries, relativistic corrections, DBOC and anharmonic-ZPE have been obtained using the CFOUR code.⁴² CCSD and CCSD(T) calculations have been performed with the MOLPRO 2024.3 computational chemistry package.^43–45 Calculations including complete treatment of triples and further (i.e. CCSDT and CCSDTQ) have been carried out using the MRCC code.^46,47

Due to their single-reference nature, coupled cluster methods are expected to fail for the processes of bond-forming and bond-breaking. As will be detailed in Section 2.2, our kinetics calculations require energy and frequency calculations along the bond dissociation or formation. Therefore, additional calculations are necessary to describe the breaking and forming of bonds. To this effect, geometry optimizations and energy calculations were carried out at the MRCI³² level using a multiconfiguration self-consistent field (MCSCF)⁴⁸ wavefunction as the reference. We used a full valence complete active space, i.e. 10 electrons in 10 orbitals. Davidson corrections (+Q)⁴⁹ were applied to the MRCI energies to improve the accuracy of the total electronic energies calculated. The correlation-consistent aug-cc-pVTZ basis set of Dunning et al.^35,36 was used. This level of calculation is denoted as MRCI+Q/aug-cc-pVTZ and has been used for each bond scan vide infra. The multireference calculations were performed using MOLPRO 2024.3.^43–45 For each stationary point, a harmonic frequency calculation was performed to confirm its nature. The minima were identified by the absence of imaginary frequencies, whereas the transition states were confirmed by the presence of exactly one imaginary frequency along the reaction coordinate.

2.2 Kinetics treatment

The C + CH₂ reaction is expected to show two important properties that need to be taken into account when considering the kinetics of the system: first, the possibility for the isomerization of the formed adduct to compete with the hydrogen abstraction that leads to C₂H + H formation. Second, simple bond-forming and bond-breaking processes can be expected to proceed without passing through a saddle point on the PES, making it challenging to determine the kinetic bottleneck.

To accurately describe systems where several chemical reactions compete with collisional (de-)energization processes, the ME^50,51 is used in conjunction with RRKM theory. The ME is a probabilistic model of energy transfer. A two-dimensional master equation is used to take both energy and total angular momentum into account. In this model, the time evolution of the probability of finding species A in an energy state E_i and angular momentum quantum number J_i is given by:

In eqn (2), four terms appear. The first two are concerned with collisional energy transfer. Specifically, the first term is the increase in population of the state E_i, J_i from collisional excitation or deexcitation of molecules of A in all energy and angular moment states E_j, J_j (including E_i, J_i), while the second one is the decrease in population of that same state due to collisions. For the evaluation of these terms, a collisional model has to be chosen. The collisional frequency (ω) used here comes from the Lennard-Jones collision frequency model.⁵² The Lennard-Jones collision parameters for the intermediates CCH₂ and HCCH were assumed to be equal and obtained from the work of Hippler et al.⁵³ Those for the bath gas, H₂, have been obtained from the work of Weaver et al.⁵⁴ The energy transfer model also has to be chosen and is in our case a simple down-exponential model.⁵⁵

The third term of eqn (2) relates to specie B forming A, and vice versa for the fourth term. The reactions are governed by microcanonical rate constants k(E,J) and conservation of the total angular momentum is accounted for.

Microcanonical rate constants k(E,J) can be obtained by RRKM theory. In such work, the rate constant for a unimolecular reaction is given by eqn (3):⁵⁶

E is the total nonfixed energy available to the molecule, J the angular momentum quantum number, L the statistical factor, and h is Planck's constant, N^†(E,J) is the sum of active vibrational and internal rotation states for the transition state, and ρ(E,J) is the density of quantum states of the reactant.

The last term of eqn (2) is a source term accounting for the bimolecular association step. R is the total rate of the entrance flux leading to the formation of the energized adduct, and f is the initial energy distribution function of the energized adduct formed by the association reaction and is given by

with N^†(E, J) here being the rovibrational sum of states of the entrance TS, T the temperature and R the ideal gas constant.

In the present work, some reactions proceed through well defined transition states, which are saddle points on the potential energy surface. On the other hand, some simple bond formation or bond breaking occur without existence of a saddle point on the potential energy surface. The first type of transition states, also called tight transition states can be treated easily. Their sum of states can be computed from the geometry and vibrational frequencies calculated at the saddle point, treating it as a rigid rotor harmonic oscillator. The densum module of the MultiWell package^57,58 was used in such cases, as well as for the computation of the density of states of the different intermediates considered. In contrast, for the second type of transition states, or loose transition states, the kinetic bottleneck is located at the point along the reaction coordinate for which the sum of states is the smallest. A suitable method to treat loose transition states is the variable reaction coordinate transition state theory (VRC-TST) approach.^59–61

In VRC-TST, the vibrational modes, perpendicular to the reaction coordinate in a bond cleavage process, are separated into two groups: conserved modes are those that remain of the same nature along the reaction coordinate, whereas transitional modes undergo a change in nature during the reaction. The former can be identified by their vibrational frequencies remaining almost constant as the interfragment distance is changed, while the latter will see their values tend to zero as the interfragment distance is increased. This is expected, as two separate fragments will have overall less vibrational modes and more rotational modes than the adduct formed by combining these two fragments. These two categories of modes are considered to be uncoupled and their overall sum of states N(E,J) is then given by a convolution as in eqn (5):

In eqn (5), E is the accessible energy and ε is the energy available for the transitional modes. N_v is the sum of states for the conserved vibrational modes, and is computed by a direct-count algorithms, while Ω_Jdε is the number of states for the transitional modes at a given J, which is computed by Monte Carlo integration.⁶⁰

To compute N(E,J), a proper potential model describing the interactions between both fragments involved is necessary. The potential energy of the system was evaluated as a function of the distance between the bonding atoms of each fragment, keeping the relative geometries of the fragments fixed, and was then fitted to a Varshni potential.^62,63 Then, an additional 6-12 Lennard-Jones potential was added between the non-bonding atoms to account for the relative orientation of the fragments. The kinetic bottleneck was then found at the inter-fragment distance which has the smallest value of N(E,J). The Variflex code⁶⁴ was used to compute the sum of states for the loose transition states. Alternatively, the variational RRKM approach implemented in the KTOOLS module of the MultiWell suite^57,58 can be employed to compute the sum of states for the loose transition state channels.

Branching ratios for individual reaction channels were then computed at the E/J resolved level using the pTS code, which solves a 2D master equation under the steady-state approximation,⁶⁵ provided the sums of states and densities of states for transition states and minima are known.

The bimolecular association step rate constant was calculated using bimolecular transition state theory as given in eqn (6):⁵⁶

At lower temperatures, numerical issues arise from the terms present in the ME model, see eqn (4). To circumvent this limitation, we decided to study the system in the collisionless limit. For that case, a computer code was developed to compute E/J resolved rate constants, provided the sums of states for the transition states are known.⁶⁶

Considering an intermediate W_i, one can write the evolution of its concentration with respect to time as:

In eqn (7), N^†_ij is the sum of states for the TS connecting species i and j, and it depends on E and J. The first term of the first sum corresponds to isomerization from the other intermediates W_j to W_i, while the second term corresponds to isomerization from W_i to W_j. The second sum corresponds to forming the product P_k from W_i, and N^†_iR × c corresponds to the formation of W_i from the reactants, with c a constant related to the reactants partition function, but also to the reactants' concentration. The various σ's are the rotational symmetry numbers. Since in the collisionless case the energy levels are uncoupled, the system is solved for each (E,J) pair. The contributions from each pair are summed in order to obtain rate constants for each product channel. Note that all microcanonical rate constants are replaced by sum of states and rotational symmetry numbers, as all densities of states will vanish.⁵⁶ In the steady-state approximation, eqn (7) is equal to 0 and one can rewrite it asor in matrix form:

[Doublestruck A]X = B ⇒ X = [Doublestruck A]⁻¹B (9)

In eqn (9), , , if i is the first intermediate formed by the association reaction, B_i = 0 otherwise. Finally, X is the vector of steady-state concentrations for all intermediates. Then, the product formation rate constant for each product channel l is given by

3 Results and discussion

In this section we will first report results from our electronic structure calculations and discuss them. Thereafter we turn to our kinetics calculations.

3.1 Electronic structure calculations

Unless explicitly stated, the energy values mentioned in this section correspond to the MRCI+Q/aug-cc-pVTZ values. These are given relative to the energy of the separated reactants C(³P) + CH₂(³B₁) and have been corrected for vibrational zero-point energy (ZPE).

The calculated stationary points on the PES associated with the C + CH₂ reaction are shown in Fig. 1, along with reaction enthalpies associated to the formation of intermediates and products at 0 K, using benchmark values from the active thermochemical tables (ATcT).⁶⁷ Each reactant has a triplet spin state as its electronic ground state. This can result in the reaction occurring on a singlet, triplet, or quintet surface. In all three cases, C and CH₂ can form CCH₂ without passing an energy barrier. CCH₂(¹A₁) is with an energy of −662 kJ mol⁻¹ more stable than CCH₂(³B₂) which has an energy of −468 kJ mol⁻¹. CCH₂(⁵A) has an energy of −141 kJ mol⁻¹. On the singlet surface, CCH₂ can undergo H-loss to form H + C₂H without going through a saddle point, or isomerize into HCCH through TS1 which has an energy of −654 kJ mol⁻¹, the latter option being more energetically favorable. HCCH, at an energy of −838 kJ mol⁻¹, can also dissociate into H + C₂H without going through a saddle point. The products H + C₂H have an energy of −298 kJ mol⁻¹, and all of the species involved on the singlet surface are therefore lower in energy than the reactants.

Fig. 1 Stationary points on the PES for the C + CH₂ system. In blue, energies at the HEAT345-Q level. In black, energies at the MRCI+Q/aug-cc-pVTZ level, corrected for harmonic ZPE. In brown, associated reaction enthalpy at 0 K calculated from the ATcT database values.

The individual contributions to the total energy in the HEAT scheme described in Section 2.1 are reported in Table 1 for the singlet surface for reactants, products, and some intermediates. In addition to stationary point calculations, surface scans were performed for the entrance (C + CH₂ → CCH₂(¹A₁)) and exit (CCH₂(¹A₁) → C₂H + H and HCCH(¹Σ_g⁺) → C₂H + H) channels. Since these steps are expected to proceed without an intrinsic energy barrier, a series of single-point energy (SPE) and frequency calculations were carried out along the reaction coordinate to accurately characterize the potential energy surface. This provides essential data to fit the potential energy used as input in VRC-TST calculations.

Table 1 HEAT345-Q contributions for stationary points on the singlet PES. All values are given in Hartrees

	CH2(^3B1)	C(^3P)	CCH2(^1A1)	TS1	HCCH(^1Σg^+)	H(^2S)	C2H(^2Σ^+)
E ^∞HF	−38.935334	−37.688648	−76.800053	−76.780535	−76.855469	−0.499994	−76.167291
ΔE^∞CCSD(T)	−0.213427	−0.156099	−0.462735	−0.478655	−0.480411		−0.446052
ΔECCSDT	−0.000470	−0.000562	−0.000703	−0.000205	0.000225		−0.000181
ΔECCSDTQ	−0.000082	−0.000029	−0.000620	−0.000791	−0.000911		−0.000939
ΔErel	−0.016051	−0.016285	−0.032230	−0.032257	−0.032142	−0.000006	−0.032363
ΔEDBOC	0.002160	0.001660	0.003766	0.003728	0.003673	0.000272	0.003522
ΔEZPE	0.016478		0.023011	0.021042	0.026349		0.013993
E HEAT	−39.146727	−37.859966	−77.269566	−77.267674	−77.338685	−0.499729	−76.629309

---

Previous works on the singlet potential energy surface of acetylene have discussed the presence or absence of an intermediate in the isomerization path between CCH₂(¹A₁) and HCCH(¹Σ_g⁺).^68–71 While we have not been able to find such an intermediate in our calculations, we infer that even if existing it would not be important to include it in our kinetics calculations as it should sit in a shallow well between two small transition states. From comparing energetic data and structure, it appears that our TS1 corresponds to the higher energy one of those two TSs and therefore TS1 should be sufficient to assess the kinetics of the system.

The initial HEAT paper reports energies for both C₂H and CH₂.⁷² For C₂H, our results are in accordance with theirs. A difference of ∼6 kJ mol⁻¹ mainly comes from our inclusion of the 2-electrons Darwin energy in the HEAT345-Q scheme that was not included in the initial HEAT scheme. Removing that contribution from our calculations lowers the energy difference to 0.4 kJ mol⁻¹. A similar observation is made for CH₂. Baraban et al.⁷³ also reported values for HCCH(¹Σ_g⁺) and CCH₂(¹A₁) using the HEAT-456QP scheme.³¹ Our total energies differ by 1.4 kJ mol⁻¹ for HCCH(¹Σ_g⁺) and by 2.2 kJ mol⁻¹ for CCH₂(¹A₁), with each individual term being very close between the calculations reported by Baraban et al. and ours. We can conclude that the results are in good agreement, with energy differences being attributed to the slight differences between the two HEAT schemes.

The triplet surface behaves very similar to the singlet one. CCH₂ can undergo H-loss to form the products H + CCH through TS5 or isomerization to trans-HCCH through TS2, which can itself isomerize to cis-HCCH through TS3. cis-HCCH can then dissociate into the products H + C₂H through TS4. A surface scan was performed for the entrance channel C + CH₂ → CCH₂(³B₂) to fit the energies to a Varshni potential. Similarly as for the singlet surface, this association proceeds without an intrinsic energy barrier.

Few studies have reported results for the triplet PES of the system,^74,75 and we could not find any existing data following the HEAT scheme. While we were able to obtain HEAT data for the minima on the triplet spin surface, our coupled-cluster calculations indicated a multireference nature of these structures through high T1 diagnostic values (∼0.08). The T1 diagnostic is regularly used to assess the reliability of coupled-cluster calculations and for any value higher than 0.02 the results should be treated as unreliable.⁷⁶ The high values obtained for the transition state calculations therefore prevent us from calculating reliable HEAT energies. Nevertheless, HEAT results are reported for the minima in Table 2, and are in general agreement with our MRCI energies as well as with energies reported in the work by Vacek et al.⁷⁴ We however note a significant difference for the ATcT value obtained for the formation of CCH₂(³B₂) and our own results which can not be easily explained.

Table 2 HEAT345-Q contributions for minima on the triplet PES. All values are given in Hartrees

	CCH2(^3B2)	trans-HCCH(^3Bu)	cis-HCCH(^3B2)
E ^∞HF	−76.747189	−76.717941	−76.751199
ΔE^∞CCSD(T)	−0.439916	−0.449692	−0.439837
ΔECCSDT	0.0123715	−0.001046	−0.001237
ΔECCSDTQ	−0.000310	−0.000890	−0.000616
ΔErel	−0.0320931	−0.032212	−0.032138
ΔEDBOC	0.003804	0.004264	0.003849
ΔEZPE	0.030450	0.025690	0.023995
E HEAT	−77.172884	−77.171827	−77.197184

---

Finally, the quintet surface was only briefly studied. While the association of C and CH₂ occurs without an intrinsic energy barrier, both hydrogen abstraction and isomerization to HCCH involve energy barriers of about ∼200 kJ mol⁻¹ so the quintet surface is not expected to play a significant role in the kinetics of the reaction.

All the energies and geometries calculated are available in the SI, as well as the harmonic frequencies used to compute the densities and sums of states necessary for the kinetics calculations. Parameters obtained from the Varshni fits for the loose transition states are also included.

3.2 Kinetics calculations

As shown in Section 3.1, the reaction can occur without energy barrier on both the singlet and triplet surfaces. The bimolecular association rate constants for both the singlet and triplet surfaces are shown in Fig. 2. From this figure, it can be seen that the rate constant on the triplet surface is about twice as large as the one on the singlet surface. While the singlet spin surface is more attractive, the higher rate constant on the triplet spin surface is explained by the difference in electronic degeneracy between the two surfaces. Thus both surfaces may contribute significantly to the overall rate constant. We study the rate constants on the singlet and the triplet surface independently and sum them to obtain the overall results.

Fig. 2 Bimolecular association rate constant as a function of temperature for the C + CH₂ reaction. In blue, for the singlet state. In orange, for the triplet state.

In this section, kinetics calculations reported were performed using MRCI energies. The calculations were performed using both HEAT345-Q and MRCI methods on the singlet surface and no significant difference arose. Studying the system in the collisionless limit allows us to get analytical results for the rate constants for individual channels. We argue that this approximation is relevant, as pressures relevant to astrochemistry are very low, and we will show that at higher temperatures, the unimolecular products are negligible for any pressure lower than atmospheric pressure. Master equation calculations have been performed for all angular momentum quantum numbers up to 80, with a maximum energy of 80 000 cm⁻¹ and energy bins of 10 cm⁻¹.

3.2.1 Singlet surface. On the singlet surface, the reaction C + CH₂ → H + C₂H is expected to proceed following the mechanism:

In the collisionless limit, the last 2 processes relating to collisional (de)excitation are ignored.

The branching ratios obtained from the analytical results in the collisionless limit and for low pressure (10⁻⁶ atm) between 10 K and 300 K are shown in Fig. 3. The two sets of results are in good agreement, with differences under 0.005% in the product branching ratios.

Fig. 3 Branching ratios at P = 10⁻⁶ atm (dotted lines) and in the collisionless case (full lines) on the singlet spin surface. In brown and green: H + C₂H formation from the dissociation of CCH₂, in purple and pink: H + C₂H formation from the dissociation of HCCH, in red and blue: backdissociation to reactants from CCH₂.

For temperatures higher than 150 K, branching ratios for each channel were calculated at various pressures. Results are presented in the SI for pressures of 1, 10³, 10⁶ and 10⁹ atm. No significant changes are observed at pressures lower than 1 atm or higher than 10⁹ atm. At pressures lower than 1 atm, product formation is the major channel, while backdissociation and complex stabilization can be neglected. With increasing pressure, complex formation becomes gradually more important. The formation of HCCH competes with formation of bimolecular products at 10³ atm, and with the formation of CCH₂ at 10⁶ atm. At extremely high pressures, CCH₂ becomes the only relevant product as any CCH₂* adduct formed is thermalized by collisions before having time to react further.

The negligible amount of intermediates stabilized at low pressures along with the strong agreement for branching ratios between the collisionless and low pressure results, shown in Fig. 3 for the temperature range 150–300 K, justify the use of the collisionless approximation to obtain results also at temperatures below 150 K. At pressures of the order of 10⁶ atm, H-loss from CCH₂ has a significant impact on the formation of H + C₂H, which is thus an important channel to include in order to have a proper description of the system. In the collisionless case, we expect that the inclusion of this channel is not important, as the isomerization and subsequent dissociation as seen in Fig. 3 dominates over the backdissociation.

3.2.2 Triplet surface. Based on our calculated stationary points for the triplet surface, the reaction C + CH₂ → C₂H + H is expected to proceed following the mechanism:

In the collisionless limit, the last 3 processes relating to collisional (de)excitation are ignored.

As in the singlet case, at any pressure relevant for astrochemistry, the stabilization of intermediates is negligible and the results in the collisionless case are used for all temperatures. The branching ratios for product formation and backdissociation are shown in Fig. 4. The results are similar to those observed on the singlet spin surface. H-loss from CCH₂(³B₂) is the main product channel, followed by loss from c-HCCH(³B₂). The higher yield from CCH₂(³B₂) is explained by the lower barrier to dissociate than the barrier to isomerize to t-HCCH(³B_u). Dissociation back to the reactants is again negligible.

Fig. 4 Branching ratios at P = 10⁻⁶ atm (dotted lines) and in the collisionless case (full lines) on the triplet surface. In brown and green: H + C₂H formation from the dissociation of CCH₂, in purple and pink: H + C₂H formation from the dissociation of HCCH, in red and blue: backdissociation to reactants from the dissociation of CCH₂.

3.2.3 Overall results and astrochemical modeling. As shown in the previous sections, backdissociation is negligible over the temperature range of interest, and the product formation rate constant is almost equal to the bimolecular association rate constant. Fig. 5 shows that around two-thirds of the product formation comes from the reaction on the triplet spin surface. The rate constant does not show any negative temperature dependence at lower temperatures as can be seen in some radical–radical reactions.^61,77,78

Fig. 5 Calculated product formation rate constants for the CH₂ + C → H + C₂H reaction: singlet case (purple), triplet (orange), overall (blue).

The overall product formation rate constant was fitted to the modified Arrhenius expression

k(T) = α × 10⁻¹⁰(T/300)^β exp(−γ/T) cm³ s⁻¹ (13)

for three different temperature ranges: 10–100 K, 100–200 K, and 200–300 K. Parameters corresponding to the best fits are given in Table 3.

Table 3 Parameters obtained from the best fit for the overall product formation rate constant

Temperature range (K)	α (cm^3 s^−1)	β	γ (K)
10–100	1.890 × 10^−10	−0.0052	23.432
100–200	2.193 × 10^−10	0.303	−2.357
200–300	2.280 × 10^−10	0.244	10.204

---

Finally, the Nautilus simulation code⁶ was used to vizualize the impact of replacing the rate constant that is currently incorporated in the kida.uva.2024 network¹⁰ by the rate constant calculated in this work. Nautilus is a 3-phase (gas, dust grain ice surface and dust grain ice mantle) time dependent chemical model for astrochemistry. There are 800 individual species included in the network that are involved in approximately 9000 reactions. Elements are initially in either their neutral or ionic forms in this model (elements with an ionization potential <13.6 eV are considered to be fully ionized). This model is used to simulate the abundances of atoms and molecules with respect to hydrogen, in molecular clouds as a function of time.

Fig. 6 shows the relative abundances of C₂H and CH₂ with the currently included value of k(T) = 10⁻¹⁰ cm³ s⁻¹ and our calculated value for a simulation at T = 10 K. The effects on the abundances of the species involved in the reaction are rather small, which is explained by the fact that at T = 10 K, the fitted rate constant obtained from the calculations in this work is k(10) = 2.168 × 10⁻¹¹ cm³ s⁻¹ which is less than one order of magnitude different from the value 10⁻¹⁰ cm³ s⁻¹ used originally in the network.

Fig. 6 Abundances of C₂H (red) and CH₂ (green) relative to H₂ with the rate constant currently in the kida.uva.2024 network (full lines) and the rate constant calculated in this work (dotted lines). The Nautilus simulation code was used with T = 10 K.

Further kinetic results obtained in the collisionless case for the singlet and triplet spin surface, as well as the overall results are available in the SI.

4 Conclusion and outlook

In this work, we have conducted a comprehensive theoretical investigation of the C + CH₂ → H + C₂H reaction, focusing on interstellar conditions. By combining high-accuracy ab initio calculations with kinetic modeling approaches, we have characterized the relevant stationary points on the potential energy surfaces for both singlet and triplet reaction pathways and derived reliable rate constants over a broad temperature range.

Notably, this work reports for the first time HEAT-scheme energy values for the minima on the triplet potential energy surface of the C₂H₂ system, contributing new high-level reference data to the field.

The main outcome of this study for astrochemistry is the derivation of three fitted expressions between 10 K and 300 K for the rate constant for the C + CH₂ → H + C₂H reaction.

Overall, this study highlights the value of theoretical chemistry in bridging experimental limitations, particularly for radical–radical reactions. Future work could use these methods to explore pressure-dependent effects in planetary atmospheres or to investigate other key reactions that may not be experimentally unfeasible.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp02644j.

Acknowledgements

Financial support from the Knut and Alice Wallenberg foundation is gratefully acknowledged through grant nr. KAW 2020.0081. M. L. is grateful for financial support from the Light S&T Graduate Program of the University of Bordeaux, the Erasmus+ exchange program, and the Adlerbertska foundation. The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through grant agreement no. 2022-06725.

Notes and references

1. M. Araki, Molecules in space, https://molecules-in.space/, Accessed: 2025-05-23.
2. B. A. McGuire, Astrophys. J., Suppl. Ser., 2022, 259, 30.
3. P. Swings and L. Rosenfeld, Astrophys. J., 1937, 86, 483–486.
4. J. Cami, J. Bernard-Salas, E. Peeters and S. E. Malek, Science, 2010, 329, 1180–1182.
5. B. A. McGuire, R. A. Loomis, A. M. Burkhardt, K. L. K. Lee, C. N. Shingledecker, S. B. Charnley, I. R. Cooke, M. A. Cordiner, E. Herbst and S. Kalenskii, et al. , Science, 2021, 371, 1265–1269.
6. M. Ruaud, V. Wakelam and F. Hersant, Mon. Not. R. Astron. Soc., 2016, 459, 3756–3767.
7. T. Millar, C. Walsh, M. Van de Sande and A. Markwick, Astron. Astrophys., 2024, 682, A109.
8. J. Holdship, S. Viti, I. Jiménez-Serra, A. Makrymallis and F. Priestley, Astron. J., 2017, 154, 38.
9. T. P. Snow and B. J. McCall, Annu. Rev. Astron. Astrophys., 2006, 44, 367–414.
10. V. Wakelam, P. Gratier, J.-C. Loison, K. Hickson, J. Penguen and A. Mechineau, Astron. Astrophys., 2024, 689, A63.
11. A. Saltelli, J. Geophys. Res.:Atmos., 1999, 104, 3789–3793.
12. A. N. Byrne, C. Xue, T. Van Voorhis and B. A. McGuire, Phys. Chem. Chem. Phys., 2024, 26, 26734–26747.
13. B. Rowe, G. Dupeyrat, J. Marquette and P. Gaucherel, J. Chem. Phys., 1984, 80, 4915–4921.
14. T. Phillips and P. Huggins, Astrophys. J., 1981, 251, 533–540.
15. F. Bensch, Astron. Astrophys., 2006, 448, 1043–1060.
16. P. Schilke, J. Keene, J. Le Bourlot, G. Pineau des Forêts and E. Roueff, Astron. Astrophys., 1995, 294, L17–L20.
17. R. Stark, P. R. Wesselius, E. F. van Dishoeck and R. Laureijs, Astron. Astrophys., 1996, 311, 282–290.
18. E. Iglesias, Astrophys. J., 1977, 218, 697–715.
19. J. M. Hollis, P. R. Jewell and F. J. Lovas, Astrophys. J., 1995, 438, 259–264.
20. A. Jacob, K. Menten, Y. Gong, P. Bergman, M. Tiwari, S. Brünken and A. Olofsson, Astron. Astrophys., 2021, 647, A42.
21. E. T. Polehampton, K. Menten, S. Brünken, G. Winnewisser and J.-P. Baluteau, Astron. Astrophys., 2005, 431, 203–213.
22. L. Vejby-Christensen, L. H. Andersen, O. Heber, D. Kella, H. B. Pedersen, H. Schmidt and D. Zajfman, Astrophys. J., 1997, 483, 531.
23. C. Romanzin, M.-C. Gazeau, Y. Bénilan, E. Hébrard, A. Jolly, F. Raulin, S. Boyé-Péronne, S. Douin and D. Gauyacq, Adv. Space Res., 2005, 36, 258–267.
24. A. Heays, A. D. V. Bosman and E. Van Dishoeck, Astron. Astrophys., 2017, 602, A105.
25. E. Reizer, B. Viskolcz and B. Fiser, Chemosphere, 2022, 291, 132793.
26. K. Tucker, M. Kutner and P. Thaddeus, Astrophys. J., 1974, 193, L115–L119 . NASA-supported research.
27. J. Lacy, N. J. Evans, J. Achtermann, D. Bruce, J. Arens and J. Carr, Astrophys. J., 1989, 342, L43–L46.
28. E. O. Pentsak, M. S. Murga and V. P. Ananikov, ACS Earth Space Chem., 2024, 8, 798–856.
29. S. Chandra, A. Kumar, B. Kumthekar and M. Sharma, New Astron., 2011, 16, 152–156.
30. P. Thaddeus, C. A. Gottlieb, R. Mollaaghababa and J. M. Vrtilek, J. Chem. Soc., Faraday Trans., 1993, 89, 2125–2129.
31. Y. J. Bomble, J. Vázquez, M. Kállay, C. Michauk, P. G. Szalay, A. G. Császár, J. Gauss and J. F. Stanton, J. Chem. Phys., 2006, 125, 064108.
32. H.-J. Werner and P. J. Knowles, J. Chem. Phys., 1988, 89, 5803–5814.
33. K. Raghavachari, G. W. Trucks, J. A. Pople and M. Head-Gordon, Chem. Phys. Lett., 1989, 157, 479–483.
34. C. Hampel, K. A. Peterson and H.-J. Werner, Chem. Phys. Lett., 1992, 190, 1–12.
35. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023.
36. R. A. Kendall, T. H. Dunning Jr and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
37. R. D. Cowan and D. C. Griffin, J. Opt. Soc. Am., 1976, 66, 1010–1014.
38. R. L. Martin, J. Phys. Chem., 1983, 87, 750–754.
39. W. Klopper, J. Comput. Chem., 1997, 18, 20–27.
40. M. Born and K. Huang, Dynamical theory of crystal lattices, Oxford University Press, 1996.
41. K. N. Rao, Molecular spectroscopy: modern research, Elsevier, 2012, pp. 115–140.
42. D. A. Matthews, L. Cheng, M. E. Harding, F. Lipparini, S. Stopkowicz, T.-C. Jagau, P. G. Szalay, J. Gauss and J. F. Stanton, J. Chem. Phys., 2020, 152, 214108.
43. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 242–253.
44. H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßelmann, D. Kats, A. Köhn, T. Korona and D. A. Kreplin, et al. , J. Chem. Phys., 2020, 152, 144107.
45. H.-J. Werner and P. J. Knowleset al. MOLPRO, 2024.3, a package of ab initio programs, see https://www.molpro.net/.
46. D. Mester, P. R. Nagy, J. Csoka, L. Gyevi-Nagy, P. B. Szabo, R. A. Horvath, K. Petrov, B. Hegely, B. Ladoczki and G. Samu, et al. , J. Phys. Chem. A, 2025, 129, 2086–2107.
47. M. Kallay, P. R. Nagy, D. Mester, L. Gyevi-Nagy, J. Csoka, P. B. Szabo, Z. Rolik, G. Samu, B. Hegely, B. Ladoczki, K. Petrov, J. Csontos, A. Ganyecz, I. Ladjanszki, L. Szegedy, M. Farkas, P. D. Mezei, R. A. Horvath and B. D. Lörincz, Mrcc, a quantum chemical program suite, see https://www.mrcc.hu.
48. H.-J. Werner and P. J. Knowles, J. Chem. Phys., 1985, 82, 5053–5063.
49. S. R. Langhoff and E. R. Davidson, Int. J. Quantum Chem., 1974, 8, 61–72.
50. J. A. Miller and S. J. Klippenstein, J. Phys. Chem. A, 2006, 110, 10528–10544.
51. S. H. Robertson, Comprehensive Chemical Kinetics, Elsevier, 2019, vol. 43, pp. 299–361.
52. J. Troe, J. Chem. Phys., 1977, 66, 4745–4757.
53. H. Hippler, J. Troe and H. Wendelken, J. Chem. Phys., 1983, 78, 6709–6717.
54. A. B. Weaver and A. A. Alexeenko, J. Phys. Chem. Ref. Data, 2015, 44, 023103.
55. R. Gilbert, K. Luther and J. Troe, Ber. Bunsen-Ges. Phys. Chem., 1983, 87, 169–177.
56. J. Troe, J. Chem. Soc., Faraday Trans., 1994, 90, 2303–2317.
57. J. Barker, T. Nguyen, J. Stanton, C. Aieta, M. Ceotto, F. Gabas, T. Kumar, C. Li, L. Lohr and A. Maranzanaet al., Ann Arbor, Michigan, USA, 2023.
58. J. R. Barker, Int. J. Chem. Kinet., 2001, 33, 232–245.
59. S. J. Klippenstein, J. Chem. Phys., 1992, 96, 367–371.
60. S. J. Klippenstein, J. Phys. Chem., 1994, 98, 11459–11464.
61. Y. Georgievskii and S. J. Klippenstein, J. Chem. Phys., 2003, 118, 5442–5455.
62. Y. Varshni, Physics, 1957, 29, 682.
63. S. Klippenstein, L. Khundkar, A. Zewail and R. Marcus, J. Chem. Phys., 1988, 89, 4761–4770.
64. S. Klippenstein, A. Wagner, R. Dunbar, D. Wardlaw and S. Robertson, Variflex, Argonne National Laboratory, Argonne, IL, 1999.
65. T. L. Nguyen and J. F. Stanton, J. Phys. Chem. A, 2015, 119, 7627–7636.
66. M. Lorin, The code is available online on Github at the following adress: https://github.com/EldwinLorin/collisionless.
67. B. Ruscic, R. E. Pinzon, G. Von Laszewski, D. Kodeboyina, A. Burcat, D. Leahy, D. Montoy and A. F. Wagner, J. Phys.:Conf. Ser., 2005, 561.
68. S. Joseph and A. Varandas, J. Phys. Chem. A, 2010, 114, 13277–13287.
69. R. Schork and H. Köppel, Theor. Chem. Acc., 1998, 100, 204–211.
70. N.-Y. Chang, M.-Y. Shen and C.-H. Yu, J. Chem. Phys., 1997, 106, 3237–3242.
71. S. Boyé-Péronne, D. Gauyacq and J. Liévin, J. Chem. Phys., 2006, 124, 214305.
72. A. Tajti, P. G. Szalay, A. G. Császár, M. Kállay, J. Gauss, E. F. Valeev, B. A. Flowers, J. Vázquez and J. F. Stanton, J. Chem. Phys., 2004, 121, 11599–11613.
73. H. Lee, J. H. Baraban, R. W. Field and J. F. Stanton, J. Phys. Chem. A, 2013, 117, 11679–11683.
74. G. Vacek, J. R. Thomas, B. J. DeLeeuw, Y. Yamaguchi, I. Schaefer and F. Henry, J. Chem. Phys., 1993, 98, 4766–4776.
75. E. Ventura, M. Dallos and H. Lischka, J. Chem. Phys., 2003, 118, 1702–1713.
76. T. J. Lee and P. R. Taylor, Int. J. Quantum Chem., 1989, 36, 199–207.
77. J. A. Miller and S. J. Klippenstein, J. Phys. Chem. A, 2003, 107, 7783–7799.
78. L. B. Harding, Y. Georgievskii and S. J. Klippenstein, J. Phys. Chem. A, 2005, 109, 4646–4656.

---

Footnote

† Contributed equally to this work.

---