Sensitivity analysis of aromatic chemistry to gas-phase kinetics in a dark molecular cloud model

Abstract

The increasingly large number of complex organic molecules detected in the interstellar medium necessitates robust kinetic models that can be relied upon for investigating the involved chemical processes. Such models require rate coefficients for each of the thousands of reactions; the values of these are often estimated or extrapolated, leading to large uncertainties that are rarely quantified. We have performed a global Monte Carlo and a more local one-at-a-time sensitivity analysis on the gas-phase rate coefficients in a 3-phase dark cloud model. Time-dependent sensitivities have been calculated using four metrics to determine key reactions for the overall network as well as for the cyanonaphthalene molecule in particular, an important interstellar species that is severely under-produced by current models. All four metrics find that reactions involving small, reactive species that initiate hydrocarbon growth have large effects on the overall network. Cyanonaphthalene is most sensitive to a number of these reactions as well as ring-formation of the phenyl cation (C₆H₅⁺) and aromatic growth from benzene to naphthalene. Future efforts should prioritize constraining rate coefficients of key reactions and expanding the network surrounding these processes. These results highlight the strength of sensitivity analysis techniques to identify critical processes in complex chemical networks, such as those often used in astrochemical modeling.

---

1 Introduction

To date, over 300 molecules and ions have been detected in the interstellar medium (ISM), with the number rapidly growing.¹ Although such detections confirm the presence of individual molecular species and provide information about their abundances, the chemical processes surrounding these species often remain a mystery. Astrochemical kinetic models are thus typically relied upon to gain information about the underlying chemical pathways.² Such models simulate the chemical evolution of an interstellar region given a number of input parameters regarding physical conditions and chemical processes.^3–6 In larger models this can easily lead to thousands of parameters, all of which must be well-constrained for a robust and reliable model.⁷ However, the accuracy of these parameters depends on the available observational, experimental, and theoretical data, which is limited.^8,9 In particular, astrochemical models are often highly sensitive to the rate coefficients of gas-phase reactions.^8,10

The experimental measurement of rate coefficients in general is a time-consuming and difficult task that is further exacerbated by challenges in measuring these values at low temperature in the gas phase. Techniques developed to meet this challenge include ion cyclotron resonance,¹¹ flowing afterglow,¹² heavy ion storage rings,¹³ and CRESU.¹⁴ This last technique is able to study ion–neutral and neutral–neutral reactions at temperatures as low as 10 K,^15,16 but requires a very large pumping capacity.¹⁴ Theoretical treatments also exist for calculating reaction rate coefficients, however, highly-accurate ab initio potential energy surfaces are often needed, as the presence of small barriers can significantly affect calculated rate coefficients.^15,17–19 This can become further complicated by the presence of open-shell radicals, which participate in a large portion of reactions relevant to the ISM.^19,20 Furthermore, effects such as tunneling and the formation of van der Waals complexes often become significant at low temperatures.^15,21 In the absence of low-temperature rate coefficients, either from experiment or theory, values are typically estimated by inspection of similar reactions or by extrapolation of high-temperature values.²² Such estimations are common in astrochemical models but lead to large uncertainties that are often not quantified.⁸

Sensitivity analysis techniques are frequently applied to complex models to gain an understanding of how uncertainties in input parameters affect output parameters.^23,24 The application of sensitivity analysis methods has been extended to a number of research areas including combustion chemistry²⁵ and atmospheric chemistry,²⁶ as well as various astrochemical environments such as the atmospheres of Titan^8,27–29 and Neptune,³⁰ protoplanetary disks,³¹ hot cores,³² diffuse clouds,³³ and dark molecular clouds.^9,33–36 These studies have primarily focused on using the global sensitivity analysis approach to study simple and abundant molecular species, both comparing abundance uncertainties to observed values and determining key reactions. Over the past decade or more since these initial studies, modeling codes have become more robust and incorporated a wider variety of processes, such as dust grain chemistry.⁴ Likewise, increasing numbers of detected species and studied reactions requires consistent updates to chemical networks.

A large portion of molecular detections have been made in dark clouds – relatively dense regions of space with low temperatures (∼10 K) yet complex chemical inventories. In particular, the first unambiguous detections of polycyclic aromatic hydrocarbons (PAHs) in the ISM have recently been made toward the Taurus Molecular Cloud 1 (TMC-1).³⁷ PAHs are thought to be abundant throughout the interstellar medium due to ubiquitous spectral features in the infrared and optical/near-infrared.³⁸ These detections suggest formation within the cloud through “bottom-up” processes, yet current astrochemical models dramatically fail to reproduce the observed abundances.³⁷ In this work, we present the application of two sensitivity analysis techniques to a modern astrochemical model of the dark cloud TMC-1. We investigate the effects of rate coefficient uncertainties on modeled abundances of aromatic species and identify their key reactions, as well as compare the suitability of the two different techniques for this purpose.

2 Methodology

2.1 Astrochemical model

A number of modeling codes and chemical networks have been utilized for astrochemical modeling of TMC-1 and other astronomical objects. NAUTILUS,⁴ UCLCHEM,⁵ and RATE22⁶ are all publicly available rate-equation codes capable of modeling dark clouds with similar general functionality. The former two codes consider physical and chemical processes on interstellar dust grains, while the latter is gas-phase only. We chose the NAUTILUS modeling code as it is a 3-phase code with a number of features based on experimental and theoretical studies. Dust grains are separated into surface and mantle (bulk ice) phases, with separate reactions for the two phases. A general overview of the model can be seen in Fig. 1. Accretion of gas-phase chemical species onto grain surfaces, desorption of grain species into the gas-phase, and difussion between grain surface and mantle allow for transfer of species between phases. A number of desorption mechanisms are considered including thermal desorption, cosmic-ray stochastic heating, chemical desorption due to reaction exothermicity, and photodesorption via external UV and cosmic-ray induced photons. Dark cloud chemical networks are often based on the KIDA or UMIST databases. We use gas-phase and grain-phase networks based on kida.uva.2014⁷ and Ruaud et al.³⁹ respectively, which have been used extensively in combination with the NAUTILUS code and its predecessors to study dark cloud chemistry.^22,40–42 A number of extensions are included based on former works focused on large carbon-chain and aromatic species that have recently been detected in TMC-1.^37,43–48 We thus expect our results to be generally applicable to modeling of dark cloud chemistry, while the methods described in the following section are suitable for any kinetic model.

Fig. 1 A general schematic of the NAUTILUS 3-phase modeling code. This is a rate-equation model that simulates the chemical evolution of an astronomical region given a chemical network, a number of physical parameters, initial abundances, and various numerical parameters. Species and reactions are separated into gas, grain surface, and grain mantle phases with physical processes allowing for transfer between phases.

‘Typical’ dark cloud conditions were assumed, namely a gas density of 2 × 10⁴ cm⁻³,⁴⁹ a kinetic temperature of 10 K for the gas and dust grains,⁵⁰ a cosmic-ray ionization rate, ζ, of 1.3 × 10⁻¹⁷ s⁻¹,⁵¹ and a visual extinction of 10 mag for external UV photons.⁵² Chemical desorption was set to 1% efficiency, assuming that 1% of molecules formed on grain surfaces have enough energy to desorb into the gas phase.⁵³ A peak grain temperature of 70 K lasting for 1 × 10⁻⁵ s was used for the cosmic-ray heating mechanism as estimated by Hasegawa and Herbst.⁵⁴ The initial elemental abundances used and their references can be seen in Table 1 and represent the low-metal abundances from Graedel et al.⁵⁵ with a lower oxygen abundance and a greater sulfur abundance. It has been found that a C/O ratio of 1.1 significantly improves agreement between modeled and observed abundances, particularly for larger hydrocarbons such as unsaturated carbon chains.^44,47,56 Data from a recent survey of sulfur-bearing molecules in molecular clouds suggests a depletion factor of 20 relative to the solar abundance, larger than the previously used value.⁵⁷ We assume that these elements are all initially in atomic form except for hydrogen, which is primarily molecular. All input files used for the model, including the reaction networks and input parameters, are available in a Zenodo repository: https://zenodo.org/doi/10.5281/zenodo.13257328.

Table 1 Initial atomic/molecular abundances referenced to hydrogen nuclei

Species	n/nH	Ref.	Species	n/nH	Ref.
H2	0.499	Ruaud et al.^4	F	6.68 × 10^−9	Neufeld et al.^58
He	0.09	Wakelam and Herbst^59	Cl	1.00 × 10^−9	Graedel et al.^55
C	1.70 × 10^−4	Jenkins^60	Si	8.00 × 10^−9	Graedel et al.^55
N	6.20 × 10^−5	Jenkins^60	Na	2.00 × 10^−9	Graedel et al.^55
O	1.55 × 10^−4	Xue et al.^44	Mg	7.00 × 10^−9	Graedel et al.^55
S	7.50 × 10^−7	Fuente et al.^57	Fe	3.00 × 10^−9	Graedel et al.^55
P	2.00 × 10^−10	Graedel et al.^55

---

2.2 Monte Carlo analysis

Previous sensitivity analyses of gas-phase rate coefficients in astrochemical kinetic models, including those of dark molecular clouds, have primarily made use of the Monte Carlo (MC) approach.^{8,9,27,29–35,61} In this method, each rate coefficient is selected randomly from an assigned uncertainty distribution, typically a log-normal distribution centered on the nominal value. A large number of iterations are then performed to obtain a statistical spread of abundances.³⁰ Spearman rank correlation coefficients (RCCs) can then be calculated for every rate coefficient-abundance pair. This value is a measure of the monotonic correlation between two variables; essentially how often an increase is observed in both variables simultaneously.⁶² The reactions with large RCCs can be considered as the reactions whose rate coefficient uncertainties most contribute to uncertainties in modeled abundances, and thus the reactions for which further constraint of rate coefficient would be most beneficial.⁸ As a global method, the Monte Carlo approach takes into account the entire range of rate coefficient uncertainty and preserves coupling between rate coefficients.^23,28,63

In practice, determining uncertainty distributions for a large number of rate coefficients is a difficult task, particularly when many of them are based on chemical intuition or extrapolation of high-temperature data. Thorough statistical analysis of uncertainties, requiring expertise in the techniques used to obtain these values, has only been performed for a small fraction of reactions present in astrochemical networks.^8,35 As such, we choose to make a blanket assumption and use uncertainties of a factor of 2.0 for every rate coefficient. This factor of 2.0 is the assumed default for most reactions in astrochemical networks and databases.^30,34,35 We also choose to use log-uniform distributions, as they may be more appropriate when uncertainty is not well quantified.²⁵ 10 000 iterations of the model were thus performed and RCCs for every species-reaction pair were calculated for 13 time-points between 4.5 × 10⁴ and 1.1 × 10⁶ years. Monte Carlo analyses with 5000, 15 000, and 20 000 iterations were also performed to test the convergence.

2.3 One-at-a-time analysis

In contrast to the global Monte Carlo approach, local sensitivity analysis methods may be used to calculate or estimate the local derivative surrounding the nominal value. As such, these techniques take a calculus-based approach to determining sensitivity rather than a statistical approach.^23,64 The simplest approach is the one-at-a-time (OAT) method, also called the finite differences or brute force method.^23,28,65 This consists of running the model a number of times, where in each iteration one rate coefficient is slightly perturbed while the others are kept constant. Sensitivity coefficients are often calculated according to the equation:Here X_i(t) is one of the 20 abundances at time t from the ith factor of change, X₀(t) is the nominal abundance of the given species at time t, F_i is the ith factor of change, and S_i(t) is the resulting sensitivity coefficient.²⁸ While this approach has a number of benefits such as being simple and easy to interpret and providing quantitative information, it is often not applicable to nonlinear models where the local derivative is not representative of the greater uncertainty range.

The application of local sensitivity analysis to astrochemical models is somewhat limited. An OAT approach has been applied to hot-core³² and dark cloud³⁴ networks where key reactions were identified based on sensitivity coefficients using either individual rate coefficient uncertainties or an assumed factor of two as the factor of perturbation. Dobrijevic et al.²⁸ applied a number of sensitivity analysis techniques to a photochemical network of Titan's atmosphere, finding that key reactions identified by OAT and MC techniques can differ significantly when there is non-linearity in the model as a result of large uncertainties in rate coefficients. However, this problem may be bypassed through the use of multiple change factors sampled over a larger range of uncertainty, a variation on the Morris method.^63,66 We perform an OAT analysis using 20 factors of change from 0.5 to 2.0 that are equally spaced in log-scale, assuming a factor of two uncertainty in all rate coefficients as with the MC analysis. For each of the 7418 reactions in our gas-phase network, 20 iterations of the model are performed. In each iteration, the rate coefficient of that reaction is multiplied by one of the 20 change factors, and the resulting modeled abundances of all species are collected at the same 13 time points as the MC analysis. Combined with the nominal, unmodified network this results in 21 abundances for every combination of species, reaction, and time. Three sensitivity metrics are calculated from these sets of 21 abundances. The first, hereafter referred to as SC, uses eqn (1) to determine an absolute sensitivity coefficient for each factor of change and averages them for a mean sensitivity coefficient. The other two are based on standard deviation, either dividing by the mean abundance to obtain a relative standard deviation (RSD) or calculating standard deviation using log-scale abundances (log SD). Directionality of the correlation is then included by fitting a linear regression to the abundances as a function of change factor and multiplying the three sensitivity metrics by the sign of the slope.

3 Results & discussion

3.1 Overall network

The sensitivities of the overall network to individual reactions were calculated as the average absolute sensitivities of all species in the network to each reaction. This procedure was repeated for the four aforementioned metrics and for three time points of interest. Reactions were then ranked by decreasing average sensitivity to determine key reactions. The chosen time points for comparison are 106 500, 314 400, and 541 100 years, and will be referred to as 1 × 10⁵, 3 × 10⁵, and 5 × 10⁵ years respectively for simplicity. Astrochemical models often agree that the age of TMC-1 seems to be between 10⁵ and 10⁶ years based on comparison to observations, however the exact point within this range is not known. These selected time points thus represent a range of suggested chemical ages, with the later time of 5 × 10⁵ generally providing the best agreement between our model and observations of larger hydrocarbons.^44,47,48 Hereafter we refer to the rank of a reaction as its importance, with the highest rank reaction being the most important. For the purpose of comparison, we generally show the top 10 reactions, denoted as key reactions, unless otherwise specified. This is a somewhat arbitrary cutoff as there is no strict definition on what constitutes a “key” reaction but is adequate for identifying predominant reactions and comparing these reactions across different metrics and time points. Full reactions lists sorted by average sensitivity metric can be found in the Zenodo repository: https://zenodo.org/doi/10.5281/zenodo.13257328.

Fig. 2 shows the top 10 reactions for the overall network at 5 × 10⁵ years for all three OAT metrics, ordered by average relative standard deviation. All three metrics agree very well in terms of the reactions identified as key reactions and the relative ordering, with only minor differences in ranking for the bottom two reactions. Additionally, the three metrics agree well in terms of relative magnitude for the sensitivity metrics. The cosmic-ray ionization of molecular hydrogen,

H₂ + CR → H₂⁺ + e⁻, (2)

is consistently the most important reaction by a significant margin, with the average sensitivity of the next most important reaction only being ∼a third of this value. This reaction is the limiting step to formation of H₃⁺, which subsequently donates a proton to C or O to initiate hydrocarbon and oxygen chemistry respectively. This reaction also has a strong positive influence on atomic hydrogen abundance, likely due to the reaction H₂ + H₂⁺ → H + H₃⁺ that quickly follows. As a result, the abundances of many protonated and saturated species increase significantly as the cosmic ray ionization rate of H₂ increases. It is important to note that the uncertainty in this rate, along with all other chemical processes involving cosmic rays, can be thought of as a combination of physical uncertainty in the cosmic-ray flux through this region (established through the input parameter ζ) and uncertainty in the efficiency of the chemical process, including uncertainty in photodissociation/photoionization cross-section.^31,67 As we are altering the rates of individual cosmic-ray processes and not the parameter ζ, we are investigating the latter source of uncertainty. In the case of reaction (2), this arises primarily due to uncertainty in the rate of ionization via secondary electrons, which depends on the gas composition.^68,69

Fig. 2 Heat map of key reactions according to the three OAT metrics at 5 × 10⁵ years for the overall network. All reaction within the top 10 with any of the three metrics are shown. The shading of each cell corresponds to the reaction's rank, with the exact sensitivity metric also given. Reactions ranked #11 or lower are given the lightest shading.

The top 10 reactions for the overall network at 1 × 10⁵, 3 × 10⁵, and 5 × 10⁵ years according to the RSD and RCC metrics are shown in Fig. 3. Comparison of key reactions at different time points provides information about how the chemistry evolves as a function of time. The key reactions identified agree well between the OAT and MC methods at all times despite using different definitions of sensitivity. Thus the reaction rates that are correlated with the abundances of many species are generally the same rates for which changes lead to large deviations in abundances of many species. Previous comparison of a number of sensitivity analysis techniques on a small model of Titan's atmosphere found that OAT and MC methods agree well when uncertainties are small, but as uncertainties grow large the agreement can worsen.²⁸ The authors attribute this to nonlinear behavior that cannot be observed with local methods, and thus suggest that global techniques such as MC should be preferred for uncertainty analysis. The strong agreement between MC and OAT over a larger range of variation may indicate that the use of multiple change factors is able to capture nonlinear behavior, which will be discussed later. It also suggests that coupling between rate coefficients, which is preserved through the MC method but not the OAT method, is not very significant.

Fig. 3 Heat map of key reactions according to the RSD metric (top) and RCC metric (bottom) for the overall network. All reaction within the top 10 at any of the 3 time points with either metric are shown. The shading of each cell corresponds to the reaction's rank, with the exact sensitivity metric also given. Reactions ranked #11 or lower are given the lightest shading.

Over this range of times, reaction (2) is consistently the most important reaction by a large margin. At 1 × 10⁵ years the remaining key reactions are dominated by reactions involving atomic carbon. As identified in a previous sensitivity analysis by Wakelam et al.,³⁴ many of these are “early-stage” reactions at the beginning of hydrocarbon growth. The authors also highlight the reaction of atomic carbon with C₃, attributing the influence of this reaction to the catalytic conversion of C and O into CO via the reaction

O + C₄ → CO + C₃. (3)

As time progresses the magnitude of sensitivities generally decreases, suggesting that the network is approaching a steady-state solution and thus becoming less susceptible to individual changes in reaction rates. In particular the importance of reactions involving atomic carbon drastically decrease as it is converted into more complex molecules, although its reactions with H₂ and C₃ remain important. A number of reactions involving formation and destruction of HCO⁺ become important at later time, likely due to its role in protonating neutral species such as carbon-chain species. The protonation of C₃ specifically is a critical step for forming C₃H_n species as well as larger hydrocarbons. The cosmic-ray ionization of helium and the reaction

N + CN → C + N₂ (4)

notably increase in sensitivity as a function of time. The former process is the limiting step to C⁺ and atomic oxygen formation via the reaction

He⁺ + CO → O + C⁺. (5)

These species react with C_mH and C_mH₂ hydrocarbons that are abundant at later times. Silicon-containing molecules also appear to be highly sensitive to the helium cosmic-ray ionization rate, possibly due to the enhanced production of organo-silicon compounds. Reaction (4) has been previously investigated as a major mechanism for partitioning elemental nitrogen into N₂. We also find that a number of large, more-saturated hydrocarbons such as C₁₀H₈, C₆H₆, and C₉H₅⁺ are highly sensitive to this rate. Additionally, there is a strong negative correlation between this rate and the abundance of atomic carbon at later times, despite this being a significant formation mechanism for atomic carbon. This seems to suggest that this reaction has a role in the partitioning of elemental carbon into its reactive, atomic form that is a building block for large hydrocarbons.

For several of these reactions that have been studied experimentally and theoretically, Wakelam et al.³⁵ compiled this data into datasheets along with recommended rate coefficients and uncertainties. The rate coefficient for the reaction

O + CN → N + CO (6)

has been measured a number of times at room temperature, however the temperature dependence of this reaction is unknown. To our knowledge, the reactions

C + H₂ → CH₂ + hν (7)

C + CH₂ → H + CCH (8)

have not been studied experimentally or theoretically and presently use estimated values. Likewise, the majority of the presented ion–neutral reactions do not have experimentally determined constants, although these rates can often be well-predicted by capture theories.⁷⁰ Further constraining the rates of these critical reactions will be significantly beneficial toward creating more robust models of molecular clouds.

3.2 Cyanonaphthalene

While the determination of key reactions for the overall network provides information about which processes are generally important due to their effects on a large number of species, a more useful application of sensitivity analysis results may be to study the key reactions for certain species of interest. As the number of detected species grow, more and more of which are larger complex organic molecules, the challenge of modeling these species and understanding the critical pathways to their formation grows as well. Here in particular, we turn our interest toward the PAHs 1- and 2-cyanonaphthalene, which we will hereafter refer to as C₁₀H₇CN. These species, along with indene⁷¹ and cyanoindene,⁷² are the only PAHs unambiguously detected in the interstellar medium.³⁷ Previous modeling work has been unable to reproduce the observed abundances by a large margin. However, it is unknown whether the discrepancy arises from missing pathways, underestimated rate coefficients, or both, as well as where in the network these missing pathways or underestimated rates may be.

In Fig. 4, the time-dependent abundances of C₁₀H₇CN for all 10 000 MC iterations are plotted, while Table 2 contains the corresponding statistical information for C₁₀H₇CN and some precursor species. As molecular complexity increases the spread of abundances does as well. This is consistent with the increased rate coefficient uncertainty propagated as the amount of reactions necessary to form these species increases. It is evident that a factor of 2 uncertainty in all reactions is not enough to account for the observed abundance of C₁₀H₇CN, as the MC iteration with the maximum abundance model still under-produces this species by almost two and a half orders of magnitude compared to its TMC-1 observed value. This is not entirely surprising given the large discrepancies between the nominal model and observations; however it does indicate that much larger uncertainties and/or missing reactions are responsible. We thus turn our attention to the reactions identified as key processes for C₁₀H₇CN by sensitivity analysis to understand critical processes and areas of the network for future study.

Fig. 4 Abundances of C₁₀H₇CN as a function of time for each of the 10 000 MC iterations. The mean abundance curve is shown in green, along with the nominal model (no changes to rate coefficients) in blue and the 2σ distribution in yellow. The dotted red line indicates the observed TMC-1 abundance, with the shaded region signifying the 1σ error.

Table 2 Statistical analysis of MC abundances (10 000 iterations) for select species at 5 × 10⁵ years

Species	Mean abundance	+2σ	−2σ	Max abundance	Min abundance	Observed abundance
C10H7CN	1.73 × 10^−15	2.05 × 10^−14	1.46 × 10^−16	1.43 × 10^−13	3.04 × 10^−19	7.35 × 10^−11
C10H8	5.96 × 10^−14	7.54 × 10^−13	4.72 × 10^−15	5.78 × 10^−12	8.52 × 10^−17	—
C6H5	1.04 × 10^−10	5.76 × 10^−10	1.88 × 10^−11	2.08 × 10^−9	9.62 × 10^−13	—
CH2CCH	1.39 × 10^−10	3.78 × 10^−10	5.08 × 10^−11	1.03 × 10^−9	2.03 × 10^−11	1.00 × 10^−8

---

3.2.1 Key reactions and their rate coefficients. The top 10 reactions ranked by sensitivity of C₁₀H₇CN according to the three OAT metrics at 5 × 10⁵ years can be seen in Fig. 5. Fig. 6 displays these reactions at three select time points according to relative standard deviation and rank correlation coefficients. As with the overall network there is strong agreement between all four metrics. We continue to refer to the ranking of a reaction as its importance, in this section specifically the importance to C₁₀H₇CN. The cosmic-ray ionization of H₂ is by far the most important reaction at 1 × 10⁵ and 3 × 10⁵ years. At 5 × 10⁵ years there are a number of sensitivities of similar magnitude and this is no longer the case.

Fig. 5 Heat map of key reactions according to the three OAT metrics at 5 × 10⁵ years for C₁₀H₇CN. All reaction within the top 10 with any of the three metrics are shown. The shading of each cell corresponds to the reaction's rank, with the exact sensitivity metric also given. Reactions ranked #11 or lower are given the lightest shading.

Fig. 6 Heat map of key reactions according to the RSD metric (top) and RCC metric (bottom) for C₁₀H₇CN. All reaction within the top 10 at any of the 3 time points with either metric are shown. The shading of each cell corresponds to the reaction's rank, with the exact sensitivity metric also given. Reactions ranked #11 or lower are given the lightest shading.

These key reactions are a mix of “early-stage” hydrocarbon reactions, many of which dominate at early times and are also important to the overall network, and “late-stage” aromatic reactions. Reaction network diagrams centered on aromatic and C₃H_n chemistry are showcased in Fig. 7 and 8 respectively. The reactions of C with H₂, C₃, and CH₂ have strong positive correlations with C₁₀H₇CN abundance as these reactions drive formation of C₃–C₅ hydrocarbons that are aromatic precursors as can be seen in these figures. Conversely, the cosmic-ray induced photoionization of C and cosmic-ray ionization of He have strong negative correlations with C₁₀H₇CN abundance as the conversion of C to C⁺ and reaction of C₃ with He⁺ disrupt the formation of these species. The reactions involving aromatic species primarily consist of ring-formation of the phenyl cation (C₆H₅⁺) from acyclic hydrocarbons

CH₄ + C₅H₂⁺ → H + C₆H₅⁺, (9)

the formation of phenyl radical (C₆H₅) and the subsequent formation of naphthalene (C₁₀H₈) via the HAVA mechanism^73,74

C₆H₆ + CRP → H + C₆H₅ (10)

C₆H₅ + CH₂CHC₂H → H + C₁₀H₈, (11)

and direct formation of C₁₀H₇CN itself

C₁₀H₈ + CN → H + C₁₀H₇CN. (12)

The direct formation of benzene from the dissociative recombination of C₆H₇⁺ does not appear in Fig. 6 but is ranked 11th according to the RSD metric and 13th according to the RCC metric at 5 × 10⁵ years. This reaction is fast and dominates the production of benzene; however the relative importance of this process along with the competing product channel to form C₆H₂ suggests the branching ratios of these channels appear to have a moderate effect on the abundance of benzene and subsequent species.

Fig. 7 Major chemical reactions leading to formation of aromatic species in our network. A blue arrow indicates a neutral–neutral reaction, an orange arrow radiative association, a green arrow dissociative-recombination, a purple arrow ion–neutral, a red arrow photodissociation, and a brown arrow adsorption onto grains. CRP stands for cosmic-ray induced photons, and the prefix ‘s-’ designates a grain-adsorbed species.

Fig. 8 Major chemical reactions leading from atomic carbon to the formation of CH₄, C₃ hydrocarbons, and six-membered rings. A blue arrow indicates a neutral–neutral reaction, an orange arrow radiative association, and a purple arrow ion–neutral. The black arrow leading from C₆H₅⁺ to C₆H₅ is meant to show that pathways exist from the former to the latter, the exact ones being showcased in Fig. 7.

At later times, the relative importance of reactions involving atomic carbon generally decreases as was the case for the overall network. Of the reactions that rise in importance, three of them,

C + HCO⁺ → CO + CH⁺ (13)

C + H₃⁺ → H₂ + CH⁺ (14)

H₂ + CH₃⁺ → CH₅⁺ + hν (15)

are part of the pathway from atomic carbon to methane as visible in Fig. 8. As shown, methane is a key precursor to the aromatic cation C₆H₅⁺ as well as C₃ hydrocarbons that participate in the formation of aromatic species. The reaction

H + CH → C + H₂ (16)

displays a negative correlation with C₁₀H₇CN abundance due to depletion of CH. In Fig. 8, it can be seen that this radical reacts with isomers of C₃H₄ to form vinylacetylene (CH₂CHC₂H), a key component of the HAVA mechanism. Additionally, the relative importance of the four previously mentioned aromatic reactions increases from 1 × 10⁵ years to 5 × 10⁵ years as simple carbon-containing species are converted into more complex hydrocarbons. Finally, we again note that reaction (4) becomes significantly more sensitive with time, becoming the 5th most important reaction to C₁₀H₇CN at 5 × 10⁵ years. As mentioned previously, we suspect this may be due to this reaction's affect on the partitioning of carbon into larger hydrocarbons that are precursors to aromatic species. Notably, the strong positive correlation between C₁₀H₇CN abundance and the rate of this reaction further suggests that this is the predominant effect. If removal of CN was as or more important we would expect this correlation, and those of many other CN-containing species, to be negative.

Of the four aromatic reactions that the C₁₀H₇CN abundance is highly sensitive to, none of them have been studied experimentally at low temperature. The rate prefactor for reaction (10) comes from the high-temperature network of Harada et al.,⁷⁵ with the 154 nm branching ratios from Kislov et al.⁷⁶ corresponding to the internal UV field expected in dense clouds.⁶⁷ Rate coefficients for cosmic-ray induced photoprocesses are typically calculated as a function of the H₂ cosmic-ray ionization rate, the photodissociation/photoionization cross section of the molecule of interest, the total absorption cross section, and the induced UV spectrum,⁷⁷ although the particular details of this calculation are not provided. More recent calculations by Heays et al.⁶⁷ using updated cross sections yielded general agreement within a factor of 2 as well as some large deviations for specific molecules/atoms. Updated calculations have not yet been performed on C₆H₆. Reaction (11) presently uses a rate coefficient of 2.5 × 10⁻¹⁰ cm³ s⁻¹ and a branching ratio of 100% according to variational transition state theory and Rice–Ramsperger–Kessel–Marcus (RRKM) calculations at 100 K and zero pressure respectively.⁷³ Follow-up RRKM calculations at zero pressure and 60 K included radiative stabilization of intermediates, finding that the inclusion of this mechanism could favor the formation of unimolecular products over the C₁₀H₈ + H channel.⁷⁸ The direct formation of C₁₀H₇CN (reaction (12)) was estimated to occur with a rate near the collision rate as it is a reaction between an unsaturated hydrocarbon and a neutral radical where the difference in ionization energy of the former and electron affinity of the latter is less than 8.75 eV.^37,79 CRESU measurements at 15 K of the analogous reaction between C₆H₆ and CN resulted in a rate coefficient of 5.4 × 10⁻¹⁰ cm³ s⁻¹, indicating that the rate of reaction with C₁₀H₈ could be larger than the currently used value. Finally, reaction (9) originates from Herbst and Leung⁸⁰ with an overall rate coefficient of 1.0 × 10⁻⁹ cm³ s⁻¹ and 80/20 branching ratios for the H/H₂ elimination channels estimated from Anicich and Huntress.⁸¹ In addition to uncertainties in the overall rate, which may be better estimated by capture theory,⁷⁰ it is unknown whether this is a viable mechanism to formation of the cyclic C₆H₅⁺ isomer. This product channel, which requires breaking and forming a large number of bonds, is likely in competition with formation of acyclic products and its status as the dominant channel is debatable.

3.2.2 Completeness of the network. Beyond the determination of individual rate coefficients that should be further constrained, the identification of key reactions provides information about areas of the network that are incomplete. As shown in Fig. 7 and 8, a number of ring-formation mechanisms to form the first six-membered aromatic ring, primarily C₆H₅⁺, exist in the network. These include [4 + 2] reactions proposed by McEwan et al.⁸² as well as a number of [3 + 3] reactions suggested by Herbst and Leung.⁸⁰ The sensitivity analysis results clearly indicate that reaction (9), a [5 + 1] mechanism involving methane, is the most important of these. The next most important ring-forming reaction is

C₂H₂ + C₄H₃⁺ → C₆H₅⁺ + hν. (17)

Despite potential energy surface calculations confirming that this radiative association should be barrierless and efficient,⁸³ it is only rank #31 in terms of importance to C₁₀H₇CN. The remaining ring-formation mechanisms of C₆H₅⁺ are significantly less important with ranks of #161 or lower. Direct formation of C₆H₆⁸⁴

CCH + CH₂CHCHCH₂ → H + C₆H₆ (18)

and of C₆H₅^85,86

CH₂CCH + CH₂CCH → C₆H₅ + H (19)

are found to be even less important to C₁₀H₇CN abundance with ranks of #1116 and #976 respectively. There is a notable drop-off in reaction significance once neutral hydrocarbons with 3 or more hydrogen atoms (C₂H₃, C₂H₄, CH₂CCH, CH₃CCH) or hydrocarbon cations with 4 or more hydrogen atoms (C₄H₄⁺) are involved. The exception is reaction (9) involving CH₄, the simplest saturated hydrocarbon that is efficiently formed from a chain of ion–molecule reactions from CH⁺ to CH₅⁺ followed by proton donation to CO⁸⁷ as shown in Fig. 8. This trend suggests a larger problem in the network regarding limited production of more-saturated hydrocarbons that are needed for ring-formation. A large portion of detected TMC-1 molecules and ions have historically been unsaturated hydrocarbons and their derivatives; however, a number of recent detections have found high abundances of more-saturated species that cannot be reproduced by current astrochemical models.^88–92 In particular, the propargyl radical (CH₂CCH) has been detected in this cloud with a large abundance.⁹³ This species is thought to play a central role in the formation of six-membered aromatic rings from smaller hydrocarbons as indicated in Fig. 8, however the inability of the model to reproduce the observed abundance results in little contribution from CH₂CCH to ring-formation.⁴⁸ Likewise, the abundance of C₁₀H₇CN is not highly sensitive to any reaction directly involving CH₂CCH nor is the observed abundance of CH₂CCH reproduced in any of the Monte Carlo iterations. This indicates that even within a factor of two variation the present pathways to CH₂CCH and other C₃ hydrocarbons from C₁ and C₂ hydrocarbons are inadequate, and other mechanisms should be explored.

Hydrogenation of C₃H_n⁺ species via reaction with H₂ has been found to proceed efficiently from C₃⁺ to C₃H⁺ followed by radiative association to C₃H₃⁺ based on low temperature ion trap experiments.⁹⁴ C₃H and C₃H₂ can subsequently be formed by dissociative recombination of C₃H₃⁺.^95,96 Further hydrogenation to C₃H₄⁺ was found to be slow with rate coefficients on the order of 1.0 × 10⁻¹⁴ cm³ s⁻¹. Calculations on the radiative associations of C₃H₃⁺ and C₃H₅⁺ with H₂ have shown sizeable barriers in the entrance channels.⁹⁷ Similarly, selected ion flow tube experiments at room temperature were not able to detect a reaction between C₃H₃⁺ and H, establishing an upper limit of 3.0 × 10⁻¹² cm³ s⁻¹ for the rate coefficient.⁹⁸ The problem of missing information regarding C₃ and larger hydrocarbons has also been discussed in the context of modeling chemistry in Titan's atmosphere, with a number of radiative association reactions proposed.²⁹ Rate coefficients for these processes were determined based on a semi-empirical fit to previous calculations using capture rates and densities of states. These reactions are able to efficiently form more-saturated hydrocarbons such as C₃H₃ and C₃H₆ but have uncertainties in rate coefficients of a factor of 30 or larger. Grain-surface hydrogenation of C₃ up to C₃H₈ through consecutive additions of H has been found to be quickly produce C₃H_n molecules on dust grains,^22,99 although additional parameters such as binding energies and efficiency of nonthermal desorption become significant for such processes. An alternative pathway involving radiative attachment of an electron followed by associative detachment with H exists for C₄ and larger carbon chains

C_m + e⁻ → C_m⁻ + hν (20)

C_m⁻ + H → C_mH + e⁻ (21)

C_mH + e⁻ → C_mH⁻ + hν (22)

C_mH⁻ + H → C_mH₂ + e⁻. (23)

Rate coefficients for reactions (20) for 4 ≤ m ≤ 9 and (22) for 2 ≤ m ≤ 8 have been determined using statistical calculations¹⁰⁰ and phase-space theory respectively.¹⁰¹ Room temperature rate coefficients have been measured for reactions (21) and (23) for a number of carbon-chain lengths using a flowing afterglow-selected ion flow tube.¹⁰² This mechanism presently does not extend beyond C_mH₂ hydrocarbons, possibly due to competition with dissociative attachment. There is a clear need for future work regarding gas-phase formation of more-saturated hydrocarbons, particularly C₃H_n species with n ≥ 3 that can participate in ring-formation.

It can be seen in Fig. 7 that the remaining aromatic reactions that C₁₀H₇CN is highly sensitive to are the three final steps leading from C₆H₆ to C₁₀H₇CN through C₆H₅ and C₁₀H₈. Given that the two most important reactions to C₁₀H₇CN involve C₆H₅, this radical appears to be the major limiting species toward formation of PAHs in the network. C₆H₅ is predominately formed by cosmic-ray induced photodissociation of C₆H₆, with direct ring-formation from CH₂CCH recombination inefficient as previously discussed. Alternative H-abstraction mechanisms for C₆H₆ such as reaction with H, OH, CH₃, and C₂H₅ all have sizeable barriers.^103–105 Once C₆H₅ is formed, only 3% or less of it reacts with CH₂CHC₂H to form C₁₀H₈. Due to a lack of reaction with interstellar cations, the vast majority of C₆H₅ adsorbs onto dust grains where it is essentially trapped. Inclusion of additional destruction mechanisms would limit the amount of C₆H₅ trapped on dust grains but also the amount available in the gas-phase, thus it may be that mechanisms other than the HAVA mechanism are required to produce PAHs in dark clouds. For example, potential energy surface calculations of the reactions

C₆H₆ + C₄H₃⁺ → H + C₁₀H₈⁺ (24)

C₆H₆⁺ + 2C₂H → C₁₀H₈⁺ (25)

C₆H₅⁺ + 2C₂H₂ → C₁₀H₉⁺ (26)

indicate that C₁₀H₈⁺ and C₁₀H₉⁺ may be formed from barrierless ion–molecule processes involving a six-membered ring,^106,107 although experiments on this latter reaction find a small barrier for the second addition of C₂H₂ and the formation of an adduct with a four-membered ring rather than a protonated naphthalene structure.¹⁰⁸ Likewise, C₁₀H₇⁺ can rapidly undergo radiative association with H₂ to form C₁₀H₉⁺.^109,110 Neutral C₁₀H₈ may then be formed by charge transfer of C₁₀H₈⁺ or dissociative recombination of C₁₀H₉⁺.^111,112 In combustion chemistry, ion-based mechanisms to soot formation involving C₃H₃⁺ have been proposed, with a number of cyclic ions such as C₇H₇⁺, C₉H₇⁺, C₉H₉⁺, C₁₀H₈⁺, and C₁₀H₉⁺ detected as products.^113–115 Considering the wide variety of ion–molecule reactions that may be involved in interstellar aromatic chemistry and the large rate coefficients of these processes at low temperatures, further expansion and refinement of the aromatic network with a focus on these processes should be another priority for future work.

The incompleteness of the network can also be seen by comparing the magnitude of C₁₀H₇CN sensitivities in Fig. 6 to those of the average sensitivities in Fig. 3. The fact that C₁₀H₇CN is significantly more sensitive to rate coefficients at these times compared to the average species may be a consequence of it being the largest and most chemically ‘complex’ species in the network. A long sequence of reactions including carbon-chain growth, ring-formation, and the HAVA mechanism are required for it to be formed. Once C₁₀H₇CN is formed, its destruction mechanisms are limited to reaction with abundant cations (H⁺, He⁺, H₃⁺, HCO⁺, C⁺, and H₃O⁺) or adsorption onto dust grains. The former mechanism is assumed to yield C₆H₅⁺, which can reform C₁₀H₇CN. The result is that any increase to reactive flux along the chain of reactions leading to C₁₀H₇CN accumulates in the abundance of this species, making it a kinetic sink. Other processes such as subsequent HAVA to form larger PAHs, or ethynyl addition are likely competing processes that are not currently included in the network. Extension of the aromatic network should thus encompass such reactions and extend beyond C₁₀H₈ and C₁₀H₇CN.

3.2.3 Model linearity. In Fig. 9, the abundance of C₁₀H₇CN at 5 × 10⁵ years as a function of change in rate coefficient from the OAT analysis is shown for its key reactions according to the RSD metric. As can be seen, there is a near-linear relationship between C₁₀H₇CN abundance and change in rate coefficient for many of these reactions, particularly for the rank 1 reaction (10). It is clear that the relationship between the rate of reaction rank 3 (which turns out to be reaction (2), the cosmic-ray ionization of hydrogen), is highly nonlinear and non-monotonic. Decreasing this rate coefficient below the nominal value leads to a strong exponential decrease in C₁₀H₇CN abundance. Conversely, increasing this rate coefficient above the nominal value results in a stationary point of peak C₁₀H₇CN abundance corresponding to a change factor of ∼1.6, after which C₁₀H₇CN abundance begins decreasing. This is likely due to competition between hydrocarbon growth and destruction, both of which proceed via ion–neutral reactions initiated by H₃⁺ formation. A greater rate of H₂ cosmic-ray ionization allows complex hydrocarbons to be formed more quickly and to a greater degree, after which destruction via reaction with ions such as H₃⁺, HCO⁺, and N₂H⁺ take over. The ability of our OAT analysis to capture this relationships and assign appropriately large sensitivity coefficients shows that the use of multiple change factors enables the acquisition of quantitative, reaction-specific sensitivities over a wide range of variations.

Fig. 9 Abundances of C₁₀H₇CN at 5 × 10⁵ years as a function of factor of change in rate constant for the top 10 reactions according to the RSD metric. Each dot represents an abundance when the rate coefficient of the corresponding reaction has been multiplied by some factor, with the color of the dot representing the degree of change. There are 21 factors of change ranging from 0.5 to 2.

4 Conclusion

We have performed Monte Carlo and one-at-a-time sensitivity analyses of a 3-phase dense cloud astrochemical model containing a number of complex organic molecules, including aromatic species. By using multiple factors of change for the one-at-a-time analysis, a larger parameter space can be explored and nonlinear relationships can be discovered. The resulting sensitivities agree well with those determined statistically via the Monte Carlo approach. We find that reactions involved in the initiation of interstellar chemistry through the formation of small reactive species, particularly those responsible for carbon-chain growth, have the greatest average effect on modeled abundances. Between 1 × 10⁵ and 5 × 10⁵ years, the likely age range of TMC-1, the absolute sensitivities and their relative orderings can change significantly. A general decrease in sensitivity metrics with time suggests the approach of a steady-state solution, with reactions involving atomic carbon becoming relatively less important at later times due to increasing complexity of carbon-containing species. Several of these key reactions still use estimated rate coefficients with large uncertainties, and as such more accurate values would significantly benefit the model. Likewise, the use of sensitivity analysis techniques allows us to reveal interesting behavior in the model even for reactions with well-constrained rate coefficients, such as the nonlinear relationship between the rate of reaction of N with CN and the abundance of atomic carbon.

Application of these results to the PAH cyanonaphthalene reveals that it is notably sensitive to four particular aromatic reactions: (1) the ring-formation of C₆H₅⁺ from C₅H₂⁺ and CH₄, (2) the photo-dissociation of benzene to C₆H₅via the cosmic-ray induced UV field and (3) its subsequent reaction with CH₂CHC₂H to form naphthalene, and (4) direct formation of cyanonaphthalene via a CN-addition H-elimination mechanism. Consequently the major bottlenecks to C₁₀H₇CN formation appear to be formation of both the first and second aromatic ring. A lack of experimental data and robust calculations casts significant uncertainty on the rates of these reactions. It is also unclear whether the reaction between CH₄ and C₅H₂⁺ is a viable pathway to C₆H₅⁺, as well as to what degree radiative stabilization competes with the bimolecular product channel for the formation of C₁₀H₈ from the HAVA mechanism. Due to these uncertainties these reactions appear to be good candidates for future experimental and theoretical kinetic and mechanistic studies. Likewise, the high sensitivity of C₁₀H₇CN to these reactions suggests that the network may be incomplete in these areas. Alternate pathways to the first six-membered ring are hindered by inefficiencies in formation of more-saturated C₃ and C₄ hydrocarbons, possibly due to missing hydrogenation mechanisms such as radiative association with H₂ or dissociative attachment. Additional ion–neutral routes to form the second aromatic ring may be feasible and rapid under TMC-1 conditions considering their general lack of activation barriers and rate coefficients near the collision rate. The identification of key reactions and incomplete regions of chemistry provides valuable insight for future work on improving kinetic models of dark molecular clouds. Additionally, these methods are broadly applicable to other areas of astrochemical modeling, as well as to fields such as combustion chemistry and atmospheric chemistry where kinetic models are also extensively used.

Author contributions

ANB: conceptualization, methodology, investigation, formal analysis, data curation, visualization, and writing – original draft. CX: methodology and software. TVV: conceptualization, supervision, and project administration. BAM: conceptualization, resources, supervision, and project administration. All authors contributed to writing – review and editing.

Data availability

The NAUTILUS modeling code used is an open-source code that can be found at https://forge.oasu.ubordeaux.fr/LAB/astrochem-tools/pnautilus. The raw sensitivity analysis data, resulting sensitivity metrics, model inputs, and a data analysis Jupyter notebook are located in the Zenodo repository: https://doi.org/10.5281/zenodo.13257329.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Dr Valentine Wakelam for use of the NAUTILUS v1.1 code. A. N. B acknowledges support from the National Science Foundation Grant Number 2141064. B. A. M. and C. X. gratefully acknowledge support of National Science Foundation grant AST-2205126.

Notes and references

1. B. A. McGuire, Astrophys. J., Suppl. Ser., 2022, 259, 30.
2. E. Herbst and W. Klemperer, Astrophys. J., 1973, 185, 505–534.
3. R. T. Garrod, Astrophys. J., 2013, 765, 60.
4. M. Ruaud, V. Wakelam and F. Hersant, Mon. Not. R. Astron. Soc., 2016, 459, 3756–3767.
5. J. Holdship, S. Viti, I. Jiménez-Serra, A. Makrymallis and F. Priestley, Astron. J., 2017, 154, 38.
6. T. J. Millar, C. Walsh, M. V. D. Sande and A. J. Markwick, Astron. Astrophys., 2024, 682, A109.
7. V. Wakelam, J.-C. Loison, E. Herbst, B. Pavone, A. Bergeat, K. Béroff, M. Chabot, A. Faure, D. Galli, W. D. Geppert, D. Gerlich, P. Gratier, N. Harada, K. M. Hickson, P. Honvault, S. J. Klippenstein, S. D. L. Picard, G. Nyman, M. Ruaud, S. Schlemmer, I. R. Sims, D. Talbi, J. Tennyson and R. Wester, Astrophys. J., Suppl. Ser., 2015, 217, 20.
8. E. Hébrard, M. Dobrijevic, P. Pernot, N. Carrasco, A. Bergeat, K. M. Hickson, A. Canosa, S. D. Le Picard and I. R. Sims, J. Phys. Chem. A, 2009, 113, 11227–11237.
9. V. Wakelam, E. Herbst, J. L. Bourlot, F. Hersant, F. Selsis and S. Guilloteau, Astron. Astrophys., 2010, 517, A21.
10. M. Agúndez and V. Wakelam, Chem. Rev., 2013, 113, 8710–8737.
11. T. B. McMahon and J. L. Beauchamp, Rev. Sci. Instrum., 1972, 43, 509–512.
12. D. B. Dunkin, F. C. Fehsenfeld, A. L. Schmeltekopf and E. E. Ferguson, J. Chem. Phys., 1968, 49, 1365–1371.
13. W. D. Geppert, R. Thomas, A. Ehlerding, F. Hellberg, F. Österdahl, M. Af Ugglas and M. Larsson, Int. J. Mass Spectrom., 2004, 237, 25–32.
14. B. R. Rowe, G. Dupeyrat, J. B. Marquette and P. Gaucherel, J. Chem. Phys., 1998, 80, 4915.
15. I. R. Cooke and I. R. Sims, ACS Earth Space Chem., 2019, 3, 1109–1134.
16. B. M. Hays, D. Gupta, T. Guillaume, O. Abdelkader Khedaoui, I. R. Cooke, F. Thibault, F. Lique and I. R. Sims, Nat. Chem., 2022, 14, 811–815.
17. D. C. Clary, Annu. Rev. Phys. Chem., 1990, 41, 61–90.
18. I. W. M. Smith, Angew. Chem., Int. Ed., 2006, 45, 2842–2861.
19. J. Ma, H. Guo and R. Dawes, Phys. Chem. Chem. Phys., 2012, 14, 12090–12093.
20. A. N. Morozov and A. M. Mebel, Phys. Chem. Chem. Phys., 2020, 22, 6868–6880.
21. Y. Georgievskii and S. J. Klippenstein, J. Phys. Chem. A, 2007, 111, 3802–3811.
22. J.-C. Loison, M. Agúndez, V. Wakelam, E. Roueff, P. Gratier, N. Marcelino, D. N. Reyes, J. Cernicharo and M. Gerin, Mon. Not. R. Astron. Soc., 2017, 470, 4075–4088.
23. H. Rabitz, M. Kramer and D. Dacol, Annu. Rev. Phys. Chem., 1983, 34, 419–461.
24. A. Saltelli, J. Geophys. Res.: Atmos., 1999, 104, 3789–3793.
25. C. W. Gao, M. Liu and W. H. Green, Int. J. Chem. Kinet., 2020, 52, 266–282.
26. A. M. Thompson and R. W. Stewart, J. Geophys. Res.: Atmos., 1991, 96, 13089–13108.
27. N. Carrasco, O. Dutuit, R. Thissen, M. Banaszkiewicz and P. Pernot, Planet. Space Sci., 2007, 55, 141–157.
28. M. Dobrijevic, E. Hébrard, S. Plessis, N. Carrasco, P. Pernot and M. Bruno-Claeys, Adv. Space Res., 2010, 45, 77–91.
29. E. Hébrard, M. Dobrijevic, J. C. Loison, A. Bergeat, K. M. Hickson and F. Caralp, Astron. Astrophys., 2013, 552, A132.
30. M. Dobrijevic and J. Parisot, Planet. Space Sci., 1998, 46, 491–505.
31. A. I. Vasyunin, D. Semenov, T. Henning, V. Wakelam, E. Herbst and A. M. Sobolev, Astrophys. J., 2008, 672, 629–641.
32. V. Wakelam, F. Selsis, E. Herbst and P. Caselli, Astron. Astrophys., 2005, 444, 883–891.
33. A. I. Vasyunin, A. M. Sobolev, D. S. Wiebe and D. A. Semenov, Astron. Lett., 2004, 30, 566–576.
34. V. Wakelam, J.-C. Loison, E. Herbst, D. Talbi, D. Quan and F. Caralp, Astron. Astrophys., 2009, 495, 513–521.
35. V. Wakelam, I. W. M. Smith, E. Herbst, J. Troe, W. Geppert, H. Linnartz, K. Öberg, E. Roueff, M. Agúndez, P. Pernot, H. M. Cuppen, J. C. Loison and D. Talbi, Space Sci. Rev., 2010, 156, 13–72.
36. J. Heyl, T. Lamberts, S. Viti and J. Holdship, Mon. Not. R. Astron. Soc., 2023, 520, 503–512.
37. 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, S. Kalenskii, M. A. Siebert, E. R. Willis, C. Xue, A. J. Remijan and M. C. McCarthy, Science, 2021, 371, 1265–1269.
38. A. Tielens, Annu. Rev. Astron. Astrophys., 2008, 46, 289–337.
39. M. Ruaud, J. C. Loison, K. M. Hickson, P. Gratier, F. Hersant and V. Wakelam, Mon. Not. R. Astron. Soc., 2015, 447, 4004–4017.
40. T. H. G. Vidal, J.-C. Loison, A. Y. Jaziri, M. Ruaud, P. Gratier and V. Wakelam, Mon. Not. R. Astron. Soc., 2017, 469, 435–447.
41. D. Navarro-Almaida, R. L. Gal, A. Fuente, P. Rivière-Marichalar, V. Wakelam, S. Cazaux, P. Caselli, J. C. Laas, T. Alonso-Albi, J. C. Loison, M. Gerin, C. Kramer, E. Roueff, R. Bachiller, B. Commerçon, R. Friesen, S. García-Burillo, J. R. Goicoechea, B. M. Giuliano, I. Jiménez-Serra, J. M. Kirk, V. Lattanzi, J. Malinen, N. Marcelino, R. Martín-Domènech, G. M. M. Caro, J. Pineda, B. Tercero, S. P. Treviño-Morales, O. Roncero, A. Hacar, M. Tafalla and D. Ward-Thompson, Astron. Astrophys., 2020, 637, A39.
42. S. Doddipatla, G. R. Galimova, H. Wei, A. M. Thomas, C. He, Z. Yang, A. N. Morozov, C. N. Shingledecker, A. M. Mebel and R. I. Kaiser, Sci. Adv., 2021, 7, eabd4044.
43. B. A. McGuire, A. M. Burkhardt, S. Kalenskii, C. N. Shingledecker, A. J. Remijan, E. Herbst and M. C. McCarthy, Science, 2018, 359, 202–205.
44. C. Xue, E. R. Willis, R. A. Loomis, K. L. K. Lee, A. M. Burkhardt, C. N. Shingledecker, S. B. Charnley, M. A. Cordiner, S. Kalenskii, M. C. McCarthy, E. Herbst, A. J. Remijan and B. A. McGuire, Astrophys. J. Lett., 2020, 900, L9.
45. B. A. McGuire, A. M. Burkhardt, R. A. Loomis, C. N. Shingledecker, K. L. K. Lee, S. B. Charnley, M. A. Cordiner, E. Herbst, S. Kalenskii, E. Momjian, E. R. Willis, C. Xue, A. J. Remijan and M. C. McCarthy, Astrophys. J. Lett., 2020, 900, L10.
46. M. C. McCarthy, K. L. K. Lee, R. A. Loomis, A. M. Burkhardt, C. N. Shingledecker, S. B. Charnley, M. A. Cordiner, E. Herbst, S. Kalenskii, E. R. Willis, C. Xue, A. J. Remijan and B. A. McGuire, Nat. Astron., 2021, 5, 176–180.
47. R. A. Loomis, A. M. Burkhardt, C. N. Shingledecker, S. B. Charnley, M. A. Cordiner, E. Herbst, S. Kalenskii, K. L. K. Lee, E. R. Willis, C. Xue, A. J. Remijan, M. C. McCarthy and B. A. McGuire, Nat. Astron., 2021, 5, 188–196.
48. A. N. Byrne, C. Xue, I. R. Cooke, M. C. McCarthy and B. A. McGuire, Astrophys. J., 2023, 957, 88.
49. R. L. Snell, W. D. Langer and M. A. Frerking, Astrophys. J., 1982, 255, 149–159.
50. P. Pratap, J. E. Dickens, R. L. Snell, M. P. Miralles, E. A. Bergin, W. M. Irvine and F. P. Schloerb, Astrophys. J., 1997, 486, 862.
51. L. Spitzer, Jr. and M. G. Tomasko, Astrophys. J., 1968, 152, 971.
52. M. Rodríguez-Baras, A. Fuente, P. Riviére-Marichalar, D. Navarro-Almaida, P. Caselli, M. Gerin, C. Kramer, E. Roueff, V. Wakelam, G. Esplugues, S. García-Burillo, R. L. Gal, S. Spezzano, T. Alonso-Albi, R. Bachiller, S. Cazaux, B. Commercon, J. R. Goicoechea, J. C. Loison, S. P. Treviño-Morales, O. Roncero, I. Jiménez-Serra, J. Laas, A. Hacar, J. Kirk, V. Lattanzi, R. Martín-Doménech, G. Muñoz-Caro, J. E. Pineda, B. Tercero, D. Ward-Thompson, M. Tafalla, N. Marcelino, J. Malinen, R. Friesen and B. M. Giuliano, Astron. Astrophys., 2021, 648, A120.
53. R. T. Garrod, V. Wakelam and E. Herbst, Astron. Astrophys., 2007, 467, 1103–1115.
54. T. I. Hasegawa and E. Herbst, Mon. Not. R. Astron. Soc., 1993, 261, 83–102.
55. T. E. Graedel, W. D. Langer and M. A. Frerking, Astrophys. J., Suppl. Ser., 1982, 48, 321–368.
56. U. Hincelin, V. Wakelam, F. Hersant, S. Guilloteau, J. C. Loison, P. Honvault and J. Troe, Astron. Astrophys., 2011, 530, A61.
57. A. Fuente, P. Rivière-Marichalar, L. Beitia-Antero, P. Caselli, V. Wakelam, G. Esplugues, M. Rodríguez-Baras, D. Navarro-Almaida, M. Gerin, C. Kramer, R. Bachiller, J. R. Goicoechea, I. Jiménez-Serra, J. C. Loison, A. Ivlev, R. Martín-Doménech, S. Spezzano, O. Roncero, G. Muñoz-Caro, S. Cazaux and N. Marcelino, Astron. Astrophys., 2023, 670, A114.
58. D. A. Neufeld, M. G. Wolfire and P. Schilke, Astrophys. J., 2005, 628, 260.
59. V. Wakelam and E. Herbst, Astrophys. J., 2008, 680, 371.
60. E. B. Jenkins, Astrophys. J., 2009, 700, 1299.
61. V. Wakelam, E. Herbst, F. Selsis and G. Massacrier, Astron. Astrophys., 2006, 459, 813–820.
62. C. Spearman, Am. J. Psychol., 1904, 15, 72–101.
63. Global sensitivity analysis: the primer, ed. A. Saltelli, Wiley, Chichester, West Sussex, 2008.
64. A. Saltelli, M. Ratto, S. Tarantola and F. Campolongo, Chem. Rev., 2005, 105, 2811–2828.
65. T. Turányi, Reliab. Eng. Syst. Saf., 1997, 57, 41–48.
66. M. D. Morris, Technometrics, 1991, 33, 161–174.
67. A. N. Heays, A. D. Bosman and E. F. V. Dishoeck, Astron. Astrophys., 2017, 602, A105.
68. A. E. Glassgold and W. D. Langer, Astrophys. J., 1973, 186, 859.
69. A. Dalgarno, M. Yan and W. Liu, Astrophys. J., Suppl. Ser., 1999, 125, 237.
70. T. Su and W. J. Chesnavich, J. Chem. Phys., 1982, 76, 5183–5185.
71. A. M. Burkhardt, K. L. K. Lee, P. B. Changala, C. N. Shingledecker, I. R. Cooke, R. A. Loomis, H. Wei, S. B. Charnley, E. Herbst, M. C. McCarthy and B. A. McGuire, Astrophys. J. Lett., 2021, 913, L18.
72. M. L. Sita, P. B. Changala, C. Xue, A. M. Burkhardt, C. N. Shingledecker, K. L. K. Lee, R. A. Loomis, E. Momjian, M. A. Siebert, D. Gupta, E. Herbst, A. J. Remijan, M. C. McCarthy, I. R. Cooke and B. A. McGuire, Astrophys. J. Lett., 2022, 938, L12.
73. D. S. N. Parker, F. Zhang, Y. S. Kim, R. I. Kaiser, A. Landera, V. V. Kislov, A. M. Mebel and A. G. G. M. Tielens, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 53–58.
74. R. I. Kaiser and N. Hansen, J. Phys. Chem. A, 2021, 125, 3826–3840.
75. N. Harada, E. Herbst and V. Wakelam, Astrophys. J., 2010, 721, 1570.
76. V. V. Kislov, T. L. Nguyen, A. M. Mebel, S. H. Lin and S. C. Smith, J. Chem. Phys., 2004, 120, 7008–7017.
77. R. Gredel, S. Lepp, A. Dalgarno and E. Herbst, Astrophys. J., 1989, 347, 289.
78. A. M. Mebel, A. Landera and R. I. Kaiser, J. Phys. Chem. A, 2017, 121, 901–926.
79. I. W. M. Smith, A. M. Sage, N. M. Donahue, E. Herbst and D. Quan, Faraday Discuss., 2006, 133, 137–156.
80. E. Herbst and C. M. Leung, Astrophys. J., Suppl. Ser., 1989, 69, 271.
81. V. G. Anicich and W. T. Huntress, Jr., Astrophys. J., Suppl. Ser., 1986, 62, 553.
82. M. J. McEwan, G. B. I. Scott, N. G. Adams, L. M. Babcock, R. Terzieva and E. Herbst, Astrophys. J., 1999, 513, 287.
83. R. Peverati, P. P. Bera, T. J. Lee and M. Head-Gordon, Astrophys. J., 2016, 830, 128.
84. B. M. Jones, F. Zhang, R. I. Kaiser, A. Jamal, A. M. Mebel, M. A. Cordiner and S. B. Charnley, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 452–457.
85. J. A. Miller and S. J. Klippenstein, J. Phys. Chem. A, 2003, 107, 7783–7799.
86. L. Zhao, W. Lu, M. Ahmed, M. V. Zagidullin, V. N. Azyazov, A. N. Morozov, A. M. Mebel and R. I. Kaiser, Sci. Adv., 2021, 7, eabf0360.
87. E. Herbst, Int. Rev. Phys. Chem., 2017, 36, 287–331.
88. N. Marcelino, J. Cernicharo, M. Agúndez, E. Roueff, M. Gerin, J. Martín-Pintado, R. Mauersberger and C. Thum, Astrophys. J., 2007, 665, L127.
89. P. Gratier, L. Majumdar, M. Ohishi, E. Roueff, J. C. Loison, K. M. Hickson and V. Wakelam, Astrophys. J., Suppl. Ser., 2016, 225, 25.
90. J. Cernicharo, M. Agúndez, C. Cabezas, N. Marcelino, B. Tercero, J. R. Pardo, J. D. Gallego, F. Tercero, J. A. López-Pérez and P. D. Vicente, Astron. Astrophys., 2021, 647, L2.
91. J. Cernicharo, R. Fuentetaja, C. Cabezas, M. Agúndez, N. Marcelino, B. Tercero, J. R. Pardo and P. D. Vicente, Astron. Astrophys., 2022, 663, L5.
92. I. R. Cooke, C. Xue, P. B. Changala, H. T. Shay, A. N. Byrne, Q. Y. Tang, Z. T. P. Fried, K. L. K. Lee, R. A. Loomis, T. Lamberts, A. Remijan, A. M. Burkhardt, E. Herbst, M. C. McCarthy and B. A. McGuire, Astrophys. J., 2023, 948, 133.
93. M. Agúndez, N. Marcelino, C. Cabezas, R. Fuentetaja, B. Tercero, P. D. Vicente and J. Cernicharo, Astron. Astrophys., 2022, 657, A96.
94. I. Savić and D. Gerlich, Phys. Chem. Chem. Phys., 2005, 7, 1026–1035.
95. G. Angelova, O. Novotny, J. B. A. Mitchell, C. Rebrion-Rowe, J. L. Le Garrec, H. Bluhme, A. Svendsen and L. H. Andersen, Int. J. Mass Spectrom., 2004, 235, 7–13.
96. J. L. McLain, V. Poterya, C. D. Molek, D. M. Jackson, L. M. Babcock and N. G. Adams, J. Phys. Chem. A, 2005, 109, 5119–5123.
97. Z. Lin, D. Talbi, E. Roueff, E. Herbst, N. Wehres, C. A. Cole, Z. Yang, T. P. Snow and V. M. Bierbaum, Astrophys. J., 2013, 765, 80.
98. G. B. I. Scott, D. A. Fairley, C. G. Freeman, M. J. McEwan, N. G. Adams and L. M. Babcock, J. Phys. Chem. A, 1997, 101, 4973–4978.
99. K. M. Hickson, V. Wakelam and J.-C. Loison, Mol. Astrophys., 2016, 3–4, 1–9.
100. R. Terzieva and E. Herbst, Int. J. Mass Spectrom., 2000, 201, 135–142.
101. E. Herbst and Y. Osamura, Astrophys. J., 2008, 679, 1670.
102. C. Barckholtz, T. P. Snow and V. M. Bierbaum, Astrophys. J., 2001, 547, L171.
103. D. L. Baulch, C. T. Bowman, C. J. Cobos, R. A. Cox, Th. Just, J. A. Kerr, M. J. Pilling, D. Stocker, J. Troe, W. Tsang, R. W. Walker and J. Warnatz, J. Phys. Chem. Ref. Data, 2005, 34, 757–1397.
104. H.-X. Zhang, S. I. Ahonkhai and M. H. Back, Can. J. Chem., 1989, 67, 1541–1549.
105. D. Hou and X. You, Phys. Chem. Chem. Phys., 2017, 19, 30772–30780.
106. C. W. Bauschlicher Jr and A. Ricca, Chem. Phys. Lett., 2000, 326, 283–287.
107. P. Ghesquière and D. Talbi, Chem. Phys. Lett., 2013, 564, 11–15.
108. A.-R. Soliman, A. M. Hamid, P. O. Momoh, M. S. El-Shall, D. Taylor, L. Gallagher and S. A. Abrash, J. Phys. Chem. A, 2012, 116, 8925–8933.
109. V. Le Page, Y. Keheyan, V. M. Bierbaum and T. P. Snow, J. Am. Chem. Soc., 1997, 119, 8373–8374.
110. T. P. Snow, V. Le Page, Y. Keheyan and V. M. Bierbaum, Nature, 1998, 391, 259–260.
111. D. K. Bohme, Chem. Rev., 1992, 92, 1487–1508.
112. V. Le Page, Y. Keheyan, T. P. Snow and V. M. Bierbaum, J. Am. Chem. Soc., 1999, 121, 9435–9446.
113. G. Baykut, F. W. Brill and J. R. Eyler, Combust. Sci. Technol., 1986, 45, 233–243.
114. K. C. Smyth, S. G. Lias and P. Ausloos, Combust. Sci. Technol., 1982, 28, 147–154.
115. N. G. Adams, L. D. Mathews and J. David Osborne, Faraday Discuss., 2010, 147, 323–335.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp03229b

---