{'fraction_non_alphanumeric': 0.03865114048177361, 'fraction_numerical': 0.024794819867832017, 'mean_word_length': 4.473818553092182, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
HCN production in Titan's Atmosphere: Coupling quantum chemistry and disequilibrium atmospheric modeling

Hydrogen cyanide (HCN) is a critical reactive source of nitrogen for building key biomolecules relevant for the origin of life. Still, many HCN reactions remain uncharacterized by experiments and theory, and the complete picture of HCN production in planetary atmospheres is not fully understood. To improve this situation, we develop a novel technique making use of computational quantum chemistry, experimental data, and atmospheric numerical simulations. First, we use quantum chemistry simulations to explore the entire field of possible reactions for a list of primary species in N2-, CH4-, and H2-dominated atmospheres. In this process, we discover 33 new reactions with no previously known rate coefficients. From here, we develop a consistent reduced atmospheric hybrid chemical network (CRAHCN) containing experimental values when available, and our calculated rate coefficients otherwise. Next, we couple CRAHCN to a 1D chemical kinetic model (ChemKM) to compute the HCN abundance as a function of atmospheric depth on Titan. Our simulated atmospheric HCN profile agrees very well with the Cassini observations. CRAHCN contains 104 reactions however nearly all of the simulated atmospheric HCN profile can be obtained using a scaled down network of only 19 dominant reactions. From here, we form a complete picture of HCN chemistry in Titan's atmosphere, from the dissociation of the main atmospheric species, down to the direct production of HCN along 4 major channels. One of these channels was first discovered and characterized in Pearce et al. (2019) and this work.


INTRODUCTION
Hydrogen cyanide (HCN) is a fundamental molecule in the origins of life. Both nucleobases, the building blocks of DNA/RNA, and amino acids, the building blocks of proteins form in HCN reactions (Oró 1961;Miller 1957). Consequently, a terrestrial atmosphere rich in HCN may be a distinct feature of what we term a biogenic planet, i.e. a planet capable of producing key biomolecules without requiring exogenous sources (e.g. meteorites).
HCN has also been detected in the atmospheres of Pluto and Neptune with concentrations of ∼40 ppm and ∼1 ppb (parts-per-billion), respectively (Lellouch et al. 2017;Marten et al. 1993). Observations have also been used to put an upper bound of 0.1 ppb on the HCN concentration in Uranus's atmosphere (Marten et al. 1993). Finally, HCN has been tentatively detected in the exoplanet atmospheres of 55 Cancri e (Tsiaras et al. 2016) and WASP-63b (MacDonald & Madhusudhan 2017), and may have been present in the early Earth atmosphere prior to the origin of life (Catling & Kasting 2007).
Given the abundance of HCN in Titan's atmosphere, and the availability of the Cassini data, Titan is the perfect testbed for validating theoretical chemical networks for HCN production in atmospheres. In the past, large-scale networks containing 800-3000+ reactions have been paired with 1D chemical kinetic codes to calculate the HCN profile as well as other chemical profiles in Titan's atmosphere (Vuitton et al. 2019;Willacy et al. 2016;Loison et al. 2015;Li et al. 2014;Hébrard et al. 2012;Krasnopolsky 2009;Lavvas et al. 2008;Hébrard et al. 2007). Past simulations provide a reasonable agreement with the Cassini atmospheric HCN measurements, and agree on the importance of H 2 CN + H −−→ HCN + H 2 as a pathway for HCN production, and HCN photolysis (HCN + hν −−→ CN + H) as a destruction process (Hébrard et al. 2012;Loison et al. 2015;Willacy et al. 2016;Vuitton et al. 2019). Rate coefficients in these networks are typically gathered from a variety of sources with differing accuracies (e.g. experiments, theoretical simulations, similar reactions, thermodynamics). Quantum chemistry methods are also occasionally used to introduce new reaction rate coefficients (e.g. Hébrard et al. (2012); Loison et al. (2015); Vuitton et al. (2019)).
Up until this point, there have been gaps in the HCN chemical data, preventing simulations from obtaining a complete picture of HCN production and destruction in Titan's atmosphere. There is a particular absence of rate coefficient data for reactions involving excited species such as doublet nitrogen atoms ( 2 N) and singlet methylene ( 1 CH 2 ): two species that are directly produced from ultraviolet (UV) and galactic cosmic ray (GCR) dissociation in Titan's upper and lower atmosphere, respectively (Vuitton et al. 2019).
We take the next step in simulating HCN chemistry in Titan's atmosphere, by finding, calculating, and including all the missing reactions relevant to the production and destruction of HCN in this environment and any direct competing reactions. Furthermore, given the complexity and computational cost of analyzing large-scale networks, we take a different approach from past Titan simulations by building and implementing a reduced network. What we have discovered, is that HCN chemistry can really be understood with a highly reduced set of chemical reactions, and that this approach is invaluable to obtaining physical insights from the results. A reduced network approach has been taken in the past for modeling warm to hot hydrogen-dominated atmospheres, and has similarly been found to conserve most of the information and insight into the dominant production and destruction pathways of observable species (Venot et al. 2019). Moreover, shifting the focus to the development of smaller chemical networks that are less computationally demanding would make them implementable in 3D atmosphere models (e.g., Lebonnois et al. (2012)), which are the only ones that are able to reproduce seasonal effects.
To summarize our method, we develop a "bottom up" theoretical approach to analyze the main formation and destruction channels of HCN in Titan's atmosphere. We start by using computational quantum chemistry simulations to scan all possible reactions between a small set of primary species. The primary species are the main reactive constituents of Titan's atmosphere (CH 4 , N 2 , H 2 , HCN), their dissociation products, and a few key intermediates (see Table 1 for list of primary species). We then use canonical variational transition state theory (Truhlar & Garrett 1984;Pearce et al. 2019) and Rice-Ramsperger-Kassel-Marcus/master equation (RRKM/ME) theory (Forst 2003) to calculate the rate coefficients for all these reactions. Many of these reactions are discovered here for the first time. We validate these calculations by comparing with experimental values in the 32% of cases that are available.
From this, we combine our calculated rate coefficients with available experimental values to obtain a consistent and complete chemical network for the reduced set of primary species. We call this network CRAHCN (Consistent Reduced Atmospheric Hybrid Chemical Network). We then couple CRAHCN with a 1D chemical kinetic model (ChemKM) (Molaverdikhani et al. 2019) in order to simulate the production and destruction of HCN in Titan's atmosphere.
This approach is beneficial in that it allows us to a) accept or reject previously reported reactions, b) discover previously unknown and potentially very important reactions, and c) develop a very fast, accurate, and consistent code to compute HCN chemistry in atmospheres.
We began this effort in Pearce et al. (2019), where we developed a feasible and accurate method to calculate a small network of 41 unique reaction rate coefficients that are directly involved with or in competition with the production of HCN in atmospheres. We focused mainly on validating the method using reactions previously studied by theory or experiment; However 15 of the calculated rate coefficients in that work had no previously known values.
In this work, we expand the network to include CN as a primary reacting species. We also explore more deeply into the unknown territory of HCN chemistry by simulating all the possible efficient interactions between the primary molecular species in the network. In this expansion, we also include three-body reactions, where an atmospheric molecule collisionally deexcites a vibrationally excited intermediate. Next, we modify some of the reactions from Pearce et al. (2019) due to new knowledge about vibrationally excited and unstable intermediates (see theoretical case studies in Appendix materials). Finally, we include two experimental spinforbidden collisionally induced intersystem crossing reactions (Douglas et al. 2018;Suzuki et al. 1993) whose rate coefficients cannot be calculated using our theoretical method. We direct the reader to the Appendix Tables for the CRAHCN rate coefficient data.
In this process, we discover 33 brand new reactions, the majority of which are based on H 2 CN, CN, CH, and electronically excited molecules ( 2 N and 1 CH 2 ). Ultimately, we finish with a consistent reduced network containing 104 reactions, which is complete for the 14 primary species in this work. In the end, we discover that only 19 reactions are at the heart of HCN production and destruction in Titan's atmosphere.
This paper proceeds as follows: in Section 2 we describe the theoretical and computational quantum chemistry methods used to calculate the reaction rate coefficients in CRAHCN, and we outline the model pa-rameters for our atmospheric numerical simulations of Titan. Then, in Section 3 we describe the results of the rate coefficient calculations, and their conformance to any experimentally measured values. We also perform a methods comparison to compare the accuracies of our chosen computational quantum chemistry method and two other widely used methods. Next, in Section 4 we analyze the results of our four numerical models of Titan's atmosphere: our fiducial model, a model with only 19 dominant reactions, a model with no GCRs, and a model with a different input for eddy diffusion. We then compare our results to the three most recent Titan models in the literature (Vuitton et al. 2019;Willacy et al. 2016;Loison et al. 2015). In this section, we also describe the two sensitivity analyses which allow us to identify the dominant pathways to HCN production and destruction on Titan. Sensitivity analyses involve running simulations where reactions from the network are excluded I) one at a time, and II) multiple at a time. Then, in Section 5 we present a step-by-step guide to the production of HCN in Titan's atmosphere, and we discuss how CRAHCN can be used for other atmospheric models. Finally, we summarize our main conclusions in Section 6.
The Appendix materials contains a) three tables comparing calculated rate coefficients at 298 K with experimental values, b) two tables containing the CRAHCN reaction rate coefficients, expressed as their Lindemann and Arrhenius parameters for temperatures from 50-400 K, c) any experimental data for the reactions calculated in this work, d) a breakdown of the calculations for some of the non-standard reactions in CRAHCN, and e) the raw computational quantum chemistry data used for rate coefficient calculations.

METHODS
Our atmospheric model can be roughly divided into two components, a) the chemical network, which is the collection of reactions and their experimental or calculated rate coefficients, and b) the chemical kinetic code, which handles radiative transfer, molecular and eddy diffusion, molecular influx and escape, photodissociation, and GCR dissociation.

Rate Coefficient Calculations
There are 104 reactions in our reduced network, representing all efficient reactions between the 14 primary species in our network (see Table 1). Our strategy for building up the network is as follows.
Firstly, if experimental data is available, we use it in this network. This accounts for 42 of the 104 reactions.
Second, we use a standard, fast, and accurate computational quantum chemistry method com- bined with standard theoretical methods to compute all possible rate coefficients and compare our results with any experimental values. This computational method is the Becke-Half-and-Half-Lee-Yang-Parr 2 (BHandHLYP) density functional and the augmented correlation-consistent polarized valence doubleζ (aug-cc-pVDZ) basis set 3 (Becke 1993;Lee et al. 1988;Dunning, Jr. 1989;Kendall & Dunning, Jr. 1992;Woon & Dunning, Jr. 1993). We call this method BH/d for short. We show that there is a good agreement between experimental and calculated rate coefficients when using this computational method, with the majority of calculations (64%) landing within a factor of 2 of experimental values and all values landing within about an order of magnitude of experimental values.
Typical uncertainties for rate coefficients-assigned in large-scale experimental data evaluations-range from a factor of 2 to an order of magnitude (Baulch et al. 1992;Tsang & Hampson 1986). As examples in our network, Baulch et al. (1992) assign factor of 2-3 uncertainties to the rate coefficients of H + H + M −−→ H 2 + M and CH 3 + H + M −−→ CH 4 + M, and order of magnitude 2 Hartree-Fock (HF) methods tend to over-estimate energy barriers, and density functional theory (DFT) methods tend to underestimate energy barriers. BHandHLYP offers a reasonable solution by using 50% Hartree-Fock and 50% density functional theory for the exchange energy calculation. 3 The basis set is the defined space for the problem, in our case it represents the atomic orbitals. The aug-cc-pVDZ basis set includes all atomic orbitals within the electron shell that is 1 above the atom's valence shell.
uncertainties to the rate coefficients of CH 4 + CH −−→ products and CN + 4 N −−→ N 2 + C. Third, we compare the accuracy of our BH/d calculations with A) a second DFT method, and B) an ab inito method. The DFT method we use for comparison is the fairly recently developed asymptotically corrected ωB97XD functional (Chai & Head-Gordon 2008). We again use the aug-cc-pVDZ basis set, and call this method ωB/d for short. The ab inito method is coupledcluster singles and doubles 4 (CCSD) with the aug-cc-pVTZ basis set (Ramabhadran & Raghavachari 2013). For convenience, we designate this as CC/t. We will show below that ωB/d and CC/t do not necessarily improve the accuracy of our calculations in terms of agreement with experimental rate coefficients. In 8 out of 12 cases, BH/d gives the best, or equal to the best agreement with experiment in comparison with the other two methods. The other two methods give the best, or equal to the best agreement in 7 out of 12 and 6 out of 12 cases, respectively. We summarize these results in detail in Section 3.1. For recent detailed reviews of DFT and coupled-cluster theory, we refer the reader to Mardirossian & Head-Gordon (2017), and Bartlett & Musia l (2007), respectively.
We have also compared the accuracy of a variety of different methods on a benchmark reaction (BHandHLYP, CCSD, HF, M06-2x, CAM-B3LYP and B3LYP) in past work, and found the BH/d computational method provided the best accuracy (Pearce et al. 2019). We note that the CAM-B3LYP/aug-cc-pVDZ method provided comparable accuracy to BH/d, and would be an interesting method to explore for future network calculations.
Finally, we compare the difference in accuracy for our benchmark method when increasing the size of the basis set from double-ζ to triple-ζ. We will show that there is no measurable improvement in accuracy for the 12 cases chosen in this work when increasing the basis set size by this level.

One-and Two-Body Reactions
One-and two-body reaction rate coefficients are calculated using CVT. This method involves varying the reaction coordinate (e.g. a bond distance) along a minimum energy path in order to find the minimum rate coefficient. This is expressed as: (Truhlar 1985) ( 1) where k GT (T, s) is the generalized transition state theory rate coefficient, T is the temperature, and s is the reaction coordinate (e.g. bond distance).
Neglecting effects due to tunneling, the generalized transition state theory (GT) reaction rate coefficient is given by: (Eyring 1935;Truhlar 1985) k where σ is the reaction path multiplicity, k B is the Boltzmann constant (1.38×10 −23 J K −1 ), T is temperature (K), h is the Planck constant (6.63×10 −34 J·s), Q ‡ is the partition function of the transition state per unit volume (cm −3 ), with its zero of energy at the saddle point, Q i is the partition function of species i per unit volume, with its zero of energy at the equilibrium position of species i, n i is the stoichiometric coefficient of species i, N is the number of reactant species, E 0 is the difference in zero-point energies between the generalized transition state and the reactants (kJ mol −1 ) (0 for barrierless reactions), and R is the gas constant (8.314×10 −3 kJ K −1 mol −1 ). To find the location along the minimum energy path where the GT rate coefficient is smallest, we use the maximum Gibbs free energy criterion, which offers a compromise of energetic and entropic effects (Truhlar 1985;Truhlar & Garrett 1980). The partition functions per unit volume have four components, where V is the volume (cm −3 ) and the t, e, v, and r subscripts stand for translational, electronic, vibrational, and rotational, respectively. The zero-point energies, Gibbs free energies, and partition functions are calculated for each reaction along its minimum energy path using the Gaussian 09 software package (Frisch et al. 2009).
This accurate yet inexpensive method was developed in Pearce et al. (2019) and typically provides rate coefficients within an order of magnitude of their published experimental values. In this work, in order to improve accuracy by a factor of 2 on average, we modify the location of the zero of energy for the vibrational partition function in Equation 2. Now, instead of being at the bottom of the internuclear potential energy well, we place it at the first vibrational level, i.e., the zero-point level. This gives the vibrational partition function the form: where N is the number of vibrational modes, Θ n is the vibrational temperature of the n th mode (Θ n = ωn k B ), and T is temperature.
Barrierless reaction rate coefficients do not typically vary by more than a factor of a few for temperatures between 50 and 400 K (Clary 1990;Li et al. 2008;Hase et al. 1987;Jasper et al. 2007;Daranlot et al. 2013), thus temperature dependences are only calculated for the reactions with barriers. This is done by calculating the rate coefficients at 50, 100, 200, 298.15, and 400 K and then fitting the results to the modified Arrhenius expression where k(T ) is the temperature-dependent second-order rate coefficient (cm 3 s −1 ), α, β, and γ are fit parameters, and T is temperature (in K).

Three-Body Reactions
When two reactants combine to form a single product, a third body is generally required to take away some of the excess vibrational energy from the reaction product, otherwise it will dissociate (Vallance 2017). The mechanism of these three-body reactions is expressed as In other words, reactants A and B combine to form a vibrationally excited product C (ν) . This product is then collisionally deexcited by species M.
The pressure-dependent rate law for this reaction is where k([M ]) is the pressure-dependent second-order rate coefficient (cm 3 s −1 ). In the high atmospheric pressure regime, the collision rate approaches 100%, and since it cannot exceed this value, the reaction rate becomes pressure independent. Conversely, in the low atmospheric pressure regime, empirical data show that the reaction rate is linear with pressure. For this reason, the pressure dependent rate coefficient is expressed as a function of the Lindemann parameters, i.e., the high-pressure limit (k ∞ [cm 3 s −1 ]) and low-pressure limit (k 0 [cm 6 s −1 ]) rate coefficients (Carstensen & Dean 2007).
It can be shown by taking the pressure limits of this equation that in the high-pressure limit, k = k ∞ , and in the low-pressure limit, We calculate the high-pressure limit rate coefficients in the same way as we calculate two-body rate coefficients, using CVT, however in this case we make use of the ktools code of the Mulitwell Program Suite (Barker et al. 2020;Barker 2001Barker , 2009). This method employed in ktools is equivalent to our method of calculating two body rate coefficients; However, we choose to use ktools as it is convenient to building the input files required for calculating low-pressure limit rate coefficients. We find deviations of < 5% between our manual CVT calculations and those performed by ktools.
To calculate the low-pressure limit rate coefficients, we make use of the Multiwell Master Equation (ME) code, which employs RRKM theory. The ME describes the interaction between collisional energy transfer (with the atmospheric "bath" gas) and chemical reaction (Pilling & Robertson 2003). In the case of our three-body reactions, the ME contains the probabilities that our vibrationally excited product will collisionally stabilize for a given atmospheric pressure and temperature. The Multiwell ME code employs Monte Carlo sampling of the ME to build a statistical average of the possible outcomes of the reaction.
The low-pressure limit rate coefficients are then calculated using the output from these stochastic trials: (Barker 2001;Akbar Ali & Barker 2015) where k ∞ is the high-pressure limit rate coefficient, f prod is the fractional yield of the collisionally deexcited product calculated by the Multiwell ME code, and [M] is the simulated concentration (cm −3 ), which is low enough for k 0 to converge.
We used N 2 as a bath gas, as it is the primary constituent of Titan's atmosphere. The energy transfer was treated with a standard exponential-down model with < ∆E > down = 0.8 T K −1 cm −1 (Gong et al. 2015;Zhao et al. 2012). The Lennard-Jones parameters 5 for N 2 and all the products were taken from the literature (Reid et al. 1977;Welty et al. 2008;Wang et al. 2010) and can be found in Table A7.
In the diffuse, upper regions of atmospheres, vibrationally excited species produced from the combination of two reactants will typically dissociate back into the original reactants. In these cases, the three-body reactions completely describe the chemistry occurring in both diffuse and dense regions of the atmosphere. However in some cases, the favourable vibrational decay products are not the original reactants (e.g. CN + 4 N −−→ CN 2(ν) · −−→ N 2 + C). In these cases, we also include the two-body reactions to these favourable decay pathways.
We would expect these two-body reactions to be less efficient in the lower, denser regions of atmospheres where the vibrationally excited intermediates can be collisionally deexcited. Regardless, we allow these reaction rate coefficients to be independent of pressure. Since radicals are typically not very abundant in the denser regions of the atmosphere, this treatment should not produce significant error.
We verify the most efficient decay pathway of vibrationally excited molecules from previous experimental studies. For more details, we refer the reader to the theoretical case studies in the Appendix materials.
The setup parameters for our fiducial Titan atmospheric model mostly match those in Hébrard et al. (2012). This includes the atmospheric temperature, pressure, and initial molar mixing ratios of N 2 , CH 4 , H 2 , and CO, the eddy diffusion profile, the influx of H 2 O from micrometeorites (5×10 6 cm −2 s −1 ), and 21 of the photochemical reactions and cross-sections. We use the solar mean for the top-of-atmosphere radiation (with solar zenith angle of 50 • ) (Thuillier et al. 2004). We also initially included the Jean's thermal escape of H and H 2 , but found it did not significantly affect our results.
Eddy diffusion describes the turbulent mixing of molecules in an atmosphere. Its form in Hébrard et al. (2012), which was originally developed by Hörst et al. (2008), is where K o is the surface eddy coefficient (400 cm 2 s −1 ), K ∞ is the top-of-atmosphere eddy coefficient (3×10 7 cm 2 s −1 ), p is the pressure (Pa), p o = 1.77×10 3 Pa and γ = 2.
In Figure 1, we plot our fiducial eddy diffusion profile along with the profiles of the three most recent Titan  Our grid is 100 uniform layers from 0 to 1300 km. To handle the galactic cosmic ray (GCR) dissociation of N 2 and CH 4 , we use the most recent GCR models by Gronoff et al. (2011Gronoff et al. ( , 2009). However, we only include the GCR reactions which produce neutral species, as our network does not contain ions. These reactions are, We do not include condensation in our models, as we are primarily interested in the gas phase chemistry leading to the production of HCN. In addition, CRAHCN does not include the heavy hydrocarbons which produce the majority of hazes, therefore we do not include haze production/destruction.

RESULTS -RATE COEFFICIENTS
In Pearce et al. (2019), we calculated 42 reaction rate coefficients involved with or in competition with the production of HCN in early Earth and Titan atmospheres. In this work, we improve on the original network by including pressure dependence on addition reactions, by modifying the vibrational partition function to improve the accuracy of our calculations (by factor of ∼2 on average), and by removing the 2 H 2 CN −−→ HCN + H 2 CNH abstraction reaction, which our recalculations show has a large barrier. We also expand the original network to 104 reactions, which is the result of exploring the entire field of possible reactions for the list of primary species in Table 1.
In Table A1, we display the 34 pressure-independent rate coefficients from Pearce et al. (2019), recalculated with a modified vibrational partition function described in Section 2. The only pressure-independent reaction not carried over from Pearce et al. (2019) is 2 H 2 CN −−→ HCN + H 2 CNH. Our recalculations show that this abstraction reaction actually has a large barrier and is therefore too inefficient to consider (k(298 K) ∼ 10 −41 cm 3 s −1 ). 20 of these 34 reactions have experimentally measured rate coefficients, and all our calculations land within an order of magnitude of these experimental values. The majority (70%) land within a factor of 2 of the experimental values. This level of accuracy is consistent with the uncertainties assigned in large-scale experimental data evaluations (Baulch et al. 1992;Tsang & Hampson 1986).
We note however that in two cases, our chosen computational method (BH/d) predicted barriers for reactions that did not have barriers. Moreover, in two other cases, this method predicted barrierless reactions for reactions with small experimental barriers. These are limitations of our chosen method, and in these few cases, we artificially remove the barriers from these calculations, or introduce experimental barriers to these calculations, respectively.
As a result of exploring the entire field of possible reactions for the primary species in this work, we calculate the rate coefficients of 36 new two-body reactions, and 32 new three-body reactions (68 total). 33 of these new reactions have no previously known rate coefficient.
In Table A2, we display the 36 new two-body reactions, along with our calculated rate coefficients at 298 K and any experimentally measured values. Seven of these reactions have experimental values, and the majority (71%) of our calculated rate coefficients are within a factor of 4 of these values. All our calculated values are within one order of magnitude of experimental values. However, in two cases our chosen computational method predicted no barrier for reactions that have small experimental barriers. As before, we artificially introduce the experimental barriers to these calculations.
In Table A3, we display the 32 calculated low pressure (k 0 ) and high pressure (k ∞ ) limit rate coefficients at 298 K for the three-body reactions in this work, as well as any experimentally measured values. 16 of the high pressure limit rate coefficients have experimental values, and the majority of cases (69%) are within a factor of 4 of these experimental values. Again, all calculations are within about an order of magnitude of experimental values.
Our calculated third-order, low pressure limit rate coefficients are within a factor of 2 of experimental values 67% of the time, and nearly within an order of magnitude in all cases. Only in the case of CN + 4 N + M −−→ CN 2 + M is our calculated rate coefficient slightly less accurate, differing from the one experimentally measured value (Provencher & McKenney 1972) by a factor of 36. This reaction is not well studied, therefore it is possible that we are not as far off from the exact value as this discrepancy implies. Calculations at the CC/t level of theory only bring this third-order rate coefficient to within a factor of 28 of the experimental value. This reaction turns out not to be important in the story of HCN in Titan's atmosphere.
In Tables A4 and A5, we display the 104 pressureand temperature-dependent rate coefficients for the fiducial chemical network used in this study. Experimentally measured rate coefficients are used when available, which is the case for 42 reactions. Sometimes experimental values are only available for one of either the high-pressure or low-pressure limit rate coefficient, in which case we use a combination of experimental and calculated values. Our calculated values are used in the majority of the network (68%).

Methods Comparison on Dominant Reactions
In a past computational methods comparison, our calculated rate coefficient for CH 4 + H −−→ CH 3 + H 2 at the CCSD/aug-cc-pVDZ level of theory was a factor of 8 smaller than the experimental values (Pearce et al. 2019). Conversely, our calculated the rate coefficient at the BH/d level of theory was within the experimental range. This, along with speed, were major motivating factors for choosing BH/d for our large-scale theoretical chemical reaction rate study.
In our sensitivity analysis in Section 4.6, we find 19 reactions dominate the production and destruction of HCN in Titan's atmosphere. 11 of these 19 reactions have experimental rate coefficients. In the case of CH 3 + H + M −−→ CH 4 + M, both the high-and low-pressure rate coefficients are experimentally measured.
In Table 2, we compare the accuracy of the BH/d, ωB/d, and CC/t for calculating the rate coefficients of these 11 dominant reactions. Based on general agreement to the experimental values, and a correct diagnosis of the reaction barrier, each method is found to have variable accuracy. In the following paragraphs, we move down through the reactions in Table 2, commenting on some of the most important results for selected reactions.
The results of this methods comparison shows a similar level of consistency in accuracy across all three methods. For the 12 coefficients, BH/d was the most accurate or tied for the most accurate method 8 times, ωB/d was the most accurate or tied for the most accurate method 7 times, and CC/t was the most accurate or tied for the most accurate method 6 times. When comparing only the DFT methods, BH/d was more accurate than ωB/d 3 times, less accurate 2 times, and similarly accurate 7 times.
In some cases, one or more methods would miss a barrier, or find one when one should not be present. In these cases, a method that correctly diagnosed the barrier was considered more accurate than one that incorrectly diagnosed the barrier, regardless of the calculated error factor. Also, all methods that incorrectly diagnosed the barrier were considered equally accurate. CC/t incorrectly diagnosed barriers five times, BH/d incorrectly diagnosed barriers four times, and ωB/d incorrectly diagnosed barriers twice.
In the case of H 2 CN + 4 N −−→ HCN + NH, the BH/d, ωB/d, and CC/t methods compute barriers of heights E 0 ∼ 15, 2, and 23 kJ mol −1 , respectively. The one experimental measurement suggests little or no barrier Table 2. A comparison of the accuracy of three methods for calculating rate coefficeints for the dominant reactions in this study. The methods are BHandHLYP/aug-cc-pVDZ (BH/d), ωB97XD/aug-cc-pVDZ (ωB/d), and CCSD/aug-cc-pVTZ (CC/t). Only 11 of the 19 dominant reactions have experimentally measured values and can be calculated using our theoretical methods. The error factor is the multiplicative or divisional factor from the nearest experimental or suggested value. For three-body reactions, the displayed rate coefficients are either the high-pressure limit (k∞) or low-pressure limit (k0). The aug-cc-pVDZ basis set is used for all calculations. Second-order rate coefficients have units cm 3 s −1 . Third-order rate coefficients have units cm 6 s −1 .

Reaction equation
We remove the barrier from this calculation as experiments predict this reaction to be barrierless or nearly barrierless (Nesbitt et al. 1990).
See Appendix materials for more details. b We introduce an experimental barrier of 8.3 kJ mol −1 (Bullock & Cooper 1972) to this calculation as no barrier is found at this level of theory.
c We remove the barrier from this calculation as experiments and theory predict this reaction to be barrierless (Miller & Bowman 1989;Marston et al. 1989). d We introduce an experimental barrier of 6.3 kJ mol −1 (Takayanagi et al. 1999) to this calculation as no barrier is found at this level of theory.
e Simulations had sporadic convergence beyond a C-N bond distance of 2.76Å. The rate coefficient is calculated with the variational transition state at this location, which has the highest ∆G. f Simulations had sporadic convergence beyond a C-N bond distance of 2.60Å. The rate coefficient is calculated with the variational transition state at this location, which has the highest ∆G. g We remove the barrier from the rate limiting step of this calculation, i.e. CH4 -CH · −−→ C2H5 · as experiments predict this reaction to be barrierless (Berman & Lin 1983;Blitz et al. 1997;Thiesemann et al. 1997). h We remove the barrier from this calculation as experiments and theory predict this reaction to be barrierless (Brouard et al. 1989).
is present (Nesbitt et al. 1990) (see Appendix materials for more details). ωB/d computes the smallest barrier; However the rate coefficient calculated using this method is a factor of 13 smaller than the experimental value, removing the barrier brings the calculated rate coefficient to within a factor of 6 of experiment. Given these discrepancies, and the lack of theoretical studies on this reaction, we recommend both a thorough theoret-ical follow-up study and additional experimental measurements.
Both the BH/d and CC/t methods find the CH 4 -CH · −−→ C 2 H 5 · step of the CH 4 + CH −−→ C 2 H 4 + H reaction to have a barrier with a height above the reactants. This differs from our ωB/d calculation, the results of experiment, and other theoretical studies, which suggest this reaction is barrierless (Ribeiro & Mebel 2015;Wang et al. 1997;Berman & Lin 1983) (see Appendix materials for further details).
All three methods find no barrier for the 2 N + CH 4 −−→ H 3 CNH · reaction step. This step is expected to have a small barrier of E 0 = 6.3 kJ mol −1 (Takayanagi et al. 1999). Similarly, Balucani et al. (2009) did not find a barrier for this reaction using the CCSD(T)/aug-cc-pVTZ level of theory.
CC/t misdiagnoses the barrier for the three-body reaction CH 3 + H + M −−→ CH 4 + M. This reaction is barrierless (Brouard et al. 1989), however CC/t calculations estimate a barrier of E 0 ∼ 62 kJ mol −1 . Similarly, CC/t misdiagnoses the barrier for 4 N + CH 3 −−→ 3 H 3 CN · , which is also barrierless (Miller & Bowman 1989;Marston et al. 1989). The barrier height for this reaction at the CC/t level of theory is 17 kJ mol −1 . BH/d and ωB/d correctly calculate no barriers for these two reactions. Conversely, BH/d and ωB/d do not calculate barriers for CN + CH 4 −−→ HCN + CH 3 , which is expected to have a barrier of E 0 = 8.3 kJ mol −1 (Bullock & Cooper 1972). Our CC/t calculations find a barrier for this reaction of 6.7 kJ mol −1 .
Lastly, the BH/d and CC/t methods compute rate coefficients for H + H + M −−→ H 2 + M that are factors of 2 from the nearest experimental value, whereas the ωB97XD method computes a rate coefficient for this reaction that is a factor of 12 smaller than the nearest experimental value.
Overall, BH/d seems to be a reasonable choice for moving forward with a large scale atmospheric study such as ours, with typical deviations from experiment of a factor of ≤ 2. ωB/d would also have been a reasonable choice moving forward, as this method correctly diagnoses barriers more frequently than BH/d and CC/t for this sample size, while maintaining accuracy nearly within an order of magnitude of experimental values. CC/t was the stand alone most accurate method in a two cases, but it was the least accurate method in six cases. Given this, and the much higher computational cost, we do not recommend CC/t for performing rate coefficient calculations for large scale atmospheric studies such as ours. For more comprehensive reaction investigations, we recommend using multiple methods, including CCSD, ωB97XD, and BHandHLYP, to verify the presence or absence of reaction barriers.

Basis Set Comparison on Dominant Reactions
In Table 3, we compare the accuracy of the augcc-pVDZ and aug-cc-pVTZ basis sets, paired with the BHandHLYP method, on the dominant reactions in this study that have experimental rate coefficients. The intent is to see, for our chosen method, if increasing the basis set size from double-ζ to triple-ζ leads to an improvement in accuracy with respect to agreement with experimental values.
What we find, is that rate coefficients calculated at the double-ζ level are generally very close to the values calculated at the triple-ζ level. Typical differences are less than 15%. In two out of twelve cases, the rate coefficients at the double-ζ level differ from the triple-ζ values by a factor of ∼2. However, in none of these twelve cases is the rate coefficient triple-ζ level more accurate than the rate coefficient at the double-ζ level with respect to experimental agreement. For this reason, and considering the added computational cost, we do not upgrade to the triple-ζ level for our large scale atmospheric study.

CRAHCN model (fiducial)
For our fiducial model, we use the CRAHCN network and the fiducial model parameters as described in Section 2.2. In Figure 2A, we see that the HCN profile from our fiducial model agrees very well with the HCN observations in Titan's lower atmosphere, landing right in the middle of the Cassini CIRS measurements. Our fiducial profile also nails the single Cassini INMS data point at 1050 km, which is the only in situ HCN measurement of Titan's atmosphere. The trade-off in agreeing so well with the INMS measurement, is that we do not agree as well with the VIMS limb measurements, or the UVIS stellar occultation measurements. This, as can be seen in Figure 2B, is standard for current state-of-the-art Titan models.

CRAHCN19 model
Our sensitivity analyses of the CRAHCN network (discussed in Section 4.6 below) revealed that 19 reactions are predominantly involved in the production and destruction of HCN in Titan's atmosphere. For the CRAHCN19 model, we use the fiducial model parameters, and a network containing only the 19 dominant reactions out of the 104 total reactions in CRAHCN. In Figure 2A, we see the HCN profile from the CRAHCN19 model almost perfectly aligns with the profile from the CRAHCN model. Maximum deviations between these Table 3. A comparison of the accuracy of BHandHLYP/aug-cc-pVDZ, and BHandHLYP/aug-cc-pVTZ for calculating rate coefficeints for the difficult dominant reactions in this study. For three-body reactions, the displayed rate coefficients are either the high-pressure limit (k∞) or low-pressure limit (k0). Second-order rate coefficients have units cm 3 s −1 . Third-order rate coefficients have units cm 6 s −1 .

Reaction equation
a We remove the barrier from this calculation as experiments predict this reaction to be barrierless or nearly barrierless. b We introduce an experimental barrier of 8.3 kJ mol −1 (Bullock & Cooper 1972) to this calculation as no barrier is found at this level of theory. c We introduce an experimental barrier of 6.3 kJ mol −1 (Takayanagi et al. 1999) to this calculation as no barrier is found at this level of theory. d Simulations had sporadic convergence beyond a C-N bond distance of 2.76Å. The rate coefficient is calculated with the variational transition state at this location, which has the highest ∆G. e Simulations had sporadic convergence beyond a C-N bond distance of 2.73Å. The rate coefficient is calculated with the variational transition state at this location, which has the highest ∆G. f We remove the barrier from the rate limiting step of this calculation, i.e. CH4 -CH · −−→ C2H5 · as experiments predict this reaction to be barrierless (Berman & Lin 1983;Blitz et al. 1997;Thiesemann et al. 1997).
model curves in the upper atmosphere are ∼10%. The total HCN produced in the CRAHCN19 model is only 6% more than the HCN produced in the fiducial model. This is mainly due to slight deviations in the lower atmosphere between models. This result suggests that the CRAHCN19 network contains nearly all that is necessary to simulate the production of HCN in Titan-like atmospheres.

CRAHCN no GCRs model
The aim of the CRAHCN no GCRs model is to examine the sensitivity of HCN production in the lower atmosphere to GCR flux. Thus, this model is similar to our fiducial model, except that all GCR reactions are removed. In Figure 2A, we see the HCN profile from the CRAHCN no GCRs model overlaps with the fiducial model in the upper atmosphere, where no GCR reactions occur. The HCN profile in the mid-lower atmosphere is only reduced by a factor of ∼2 compared to the fiducial model. Overall, the CRAHCN no GCRs model produces about one-third as much HCN in Titan's atmosphere as the fiducial model. This result suggests that GCRs are not critical for the production of HCN in Titan's atmosphere, however they boost total HCN production by about a factor of ∼3.  Loison et al. 2015). For the literature model parameters, see Table 4.

CRAHCN Li Kzz model
To investigate the effects of eddy diffusivity in the distribution of atmospheric HCN, we present a model in which we modify the eddy diffusion profile to match that of Li et al. (2014) (see Figure 1). This profile differs from the Hörst et al. (2008) profile used in all other models in that the profile inverts at 350 km, creating a low eddy diffusion zone near 550 km. This is also the profile used in the Willacy et al. (2016) Titan model. In Figure 2A, we see that using this eddy diffusion profile reduces the HCN abundance in the lower atmosphere by a factor of ∼3, and increases the HCN abundance in the mid atmosphere by a factor of ∼4 with respect to our fiducial model. Overall, the HCN profile from this model does not agree with the Cassini CIRS data as well as our fiducial model, as the former misses the range of CIRS measurements from 200-300 km by about a factor of 3.

Comparison to Other Recent Titan Models
In Figure 2B, we plot HCN profiles from the three most recent Titan models (Vuitton et al. 2019;Willacy et al. 2016;Loison et al. 2015) to compare with our fiducial HCN profile. It is important to emphasize that the three models from the literature focused on reproducing the observed profiles of many chemical species, only one of which was HCN. The differences in parameters and chemistry between these models is summarized in Table 4 and Figure 1.
Differences between all four model curves are within a factor of ∼3 in the lower atmosphere, a factor of ∼8 in the mid atmosphere, and a factor of ∼2 in the upper atmosphere. Given the differences in eddy diffusion profiles, condensation/sedimentation, photochemistry, GCR chemistry, and reaction networks, a complete explanation on the variations between these curves is not possible, however, we note a few things below.
The  sphere. This suggests that the major differences between our fiducial HCN profile and the HCN profile in Vuitton et al. (2019) are due to differences in chemical networks and photochemistry, rather than eddy diffusion. Differences between our fiducial HCN profile and the HCN profile in Loison et al. (2015) also vary by 2. However, due to lack of data, we cannot comment on a major source of the discrepancies.
Other differences between our model and those in the literature include treatments for condensation/sedimentation and haze formation. We do not include condensation/sedimentation in our fiducial model, as we are mainly interested in the gas phase chemistry leading to the production of HCN. In addition, CRAHCN does not include the heavy hydrocarbons that produce the majority of hazes. Willacy et al. (2016) find that condensation/sedimentation only affects the HCN profile below ∼100 km, which is below any Cassini measurement. Willacy et al. (2016) also included the permanent removal of HCN via haze production in one of their models (Model C), which resulted in a reduction of HCN below ∼500 km of approximately a factor of 4 compared to their models without hazes (Models A and B). Their haze model, however, does not agree with the Cassini CIRS HCN data as well as their models without haze production.
Overall, our fiducial HCN model is in general agreement with the most recent Titan models in the literature.

Sensitivity Analyses
Not every reaction in an atmospheric chemical network contributes significantly to the production and destruction of a given species. To discover which reactions in CRAHCN contributes to the fiducial HCN profile we perform two types of sensitivity analyses on our fiducial model.

Sensitivity Analysis I
The first sensitivity analysis involves running 104 additional numerical simulations of Titan's atmosphere. In each simulation, one of the 104 reactions in CRAHCN is removed, and the resultant HCN profile compared with the fiducial HCN profile. We also perform this sensitivity analysis on the CRAHCN/no GCRs model. This first sensitivity analysis revealed 17 of the 19 dominant reactions. In Figure 3A we display the changes to the HCN profiles that occur when each of these 17 reactions are removed. The removal of all other reactions did not significantly effect the fiducial HCN profile.

Sensitivity Analysis II
The second sensitivity analysis involves running a much larger number of simulations. Starting with all 104 reactions, as was done in sensitivity analysis I, each reaction is excluded in a simulation to see how it effects the HCN profile. The reaction whose removal affects the HCN profile the least is then removed, and the process is repeated with 103 reactions. The least important reaction is removed at each stage, until the exclusion of any of the remaining reactions leads to a 10% deviation from the fiducial HCN profile.
The second sensitivity analysis revealed 2 additional dominant reactions, CH 3 + CH −−→ C 2 H 2 + H + H and 3 CH 2 + CH −−→ C 2 H + H + H, bringing the total to 19 dominant reactions. In Figure 4A we display the changes to the HCN profiles that occur when each of these 19 reactions are removed from a network containing only these 19 reactions (CRAHCN19).

19 Dominant Reactions
The 19 dominant reactions are listed in Table 5 below. Five reactions dominate the production of HCN, four are critical for increasing the feedstock of precursor molecules that react to produce HCN, one reaction dominates the destruction of HCN, seven reactions reduce the key precursor molecules that produce HCN, one reaction attenuates the precursor sinks by reducing H abundance, and one reaction acts as both a precursor source and sink.
The biggest impact to the fiducial HCN profile is the removal of CN + CH 4 −−→ HCN + CH 3 . This is the key reaction that recycles CN-primarily from HCN photodissociation-back into HCN. This reaction accounts for ∼36-46% 7 of the total HCN in Titan's atmosphere, and is dominant primarily because of Titan's high atmospheric CH 4 abundance.
The next most important channel is 2 N + CH 4 −−→ H 2 CN + H 2 , followed by H 2 CN + H −−→ HCN + H 2 . Again, due primarily to the high atmospheric CH 4 concentrations, this multi-step reaction is responsible for ∼32-38% of the total HCN in Titan's atmosphere. 7 Percent contributions for the four main HCN channels are calculated by dividing the difference in the total HCN abundance (integrated over all altitudes) when removing that reaction, by the summed up total differences in HCN abundances when removing each of the four main reactions. Calculations differ when using the fiducial and CRAHCN19 models, therefore we express the values as a range.
Finally, the final dominant reaction for HCN production in Titan's atmosphere, is 2 N + CH 3 −−→ H 2 CN + H −−→ HCN + H 2 . This reaction was discovered by Pearce et al. (2019), and had no known rate coefficient prior to that work. It accounts for ∼2% of the total HCN in Titan's atmosphere.
It is worth noting that the leading reaction, CN + CH 4 −−→ HCN + CH 3 , only produces HCN if CN is present. CN primarily comes from the photodestruction of HCN, therefore this reaction is not responsible for starting HCN synthesis in Titan's atmosphere, but rather, maintaining it. Given that we ignore this maintenance reaction, the other three leading channels produce approximately 59%, 37%, and 4% of the total initial HCN in Titan's atmosphere, respectively. There are four reactions that play an important role in processing radical species to produce the precursors for HCN production reactions. The main one in the upper atmosphere is CH 4 + 1 CH 2 −−→ CH 3 + CH 3 , which provides CH 3 for reactions 69 and 74. Although this reaction removes a CH 4 molecule, which is also a reactant for HCN production via reactions 54 and 73, the rate coefficients for these CH 4 -based reactions are 2-5 orders of magnitude smaller than those for reactions 69 and 74, and thus producing more CH 3 leads to more efficient HCN production. The other upper atmospheric precursor source is CH + H −−→ C + H 2 followed by the photolysis of H 2 to form two H atoms to be used by reaction 43.
The other two processing reactions, which were only revealed by the second sensitivity analysis, produce the H atoms necessary for H 2 CN + H −−→ HCN + H 2 . These reactions dominate in the lower atmosphere, where UV light does not reach and thus H 2 photodissociation does not occur.
The dominant sink for HCN is HCN + 2 N −−→ N 2 + CH and accounts for nearly 100% of total HCN removal. Although the photodissociation reaction HCN + hν −−→ CN + H destroys HCN efficiently, the removal of this reaction does not significantly affect the total HCN abundance in Titan's atmosphere. This is because CN efficiently reacts with CH 4 to recycle back into HCN.
Several reactions reduce HCN production by acting as sinks to important HCN precursors, i.e., CH 3 , CN, CH 4 , 1 CH 2 , 3 CH 2 , 4 N, and H. One of these precursor sink reactions (no. 37) ironically produces HCN. However, since this reaction also produces NH, the dominant effect is the removal of 4 N from the atmosphere via NH + 4 N −−→ N 2 + H. One reaction, CH 4 + CH −−→ C 2 H 4 + H, seems to act as a precursor sink in the lower atmosphere, and a precursor source in the mid-upper atmosphere. However, this changes when switching from the fiducial to the CRAHCN19 models, and therefore the true role of this reaction is uncertain. Finally, one reaction is key to attenuating the effect of a precursor sink in the lower atmosphere. H + H + N 2 −−→ H 2 + N 2 reduces the H-atom abundance to attenuate the effects of CH 3 + H + N 2 −−→ CH 4 + N 2 .

The Case of No GCRs
In Figures 3B and 4B, we display the changes to the HCN profiles that occur when each the 19 dominant reactions are removed from the CRAHCN/no GCRs and CRAHCN19/no GCRs models, respectively. As a reminder, the 19 dominant reactions are the reactions which upon their removal, have the greatest effect the HCN profile. Interestingly, only 15 reactions, upon their removal, affected the no GCR HCN profiles. The removal of CH 3 + CH −−→ C 2 H 2 + H + H, 3 CH 2 + CH −−→ C 2 H + H + H, CH 3 + H + N 2 −−→ CH 4 + N 2 and H + H + N 2 −−→ H 2 + N 2 did not significantly effect the CRAHCN/no GCRs or CRAHCN19/noGCRs HCN profiles.
These four particular reactions require high abundances of CH 3 , 3 CH 2 , and H produced by the GCR destruction of CH 4 in the lower atmosphere in order to become important for HCN production and destruction. It can be seen in Figures 3A and 4A that these four reactions have no affect on the HCN profiles in the upper atmosphere, where GCR reactions do not occur.
In the upper atmospheres, the HCN difference profiles in Figures 3B and 4B look nearly identical to the likecolored profiles in Figure 3A. Differences between any two like-colored curves are more drastic in the lower atmospheres, where GCR reactions occur. Removing GCR reactions changes the feedstock of methane and nitrogen radicals, and therefore adjusts the relative importance of each of the dominant reactions that use these radicals.
Overall, these results suggest that GCR reactions do not significantly control which reactions dominate at producing and destroying HCN in Titan's atmosphere, but they do affect the relative amount that they contribute to the overall HCN abundance.

How HCN is Produced in Titan's Atmosphere
Out of the 104 chemical reactions in CRAHCN, we find only 19 reactions significantly contribute to the production and destruction of HCN in Titan's atmosphere. Most of these reactions are direct sources and sinks for HCN and sources and sinks for the precursors to HCN (e.g. CH 3 , 4 N). In the only other case, a reaction has the role of attenuating the effects of a precursor sink.
In Figure 5, we describe the step-by-step process of HCN production in Titan's atmosphere. First, UV radiation in the upper atmosphere, and GCRs in the lower atmosphere, break apart CH 4 , N 2 , and H 2 into reactive high-energy radical species. Second, these radicals get processed via chemical reactions to form HCN precursors (e.g. CH 3 and H). These processing reactions differ in the upper and lower atmosphere. For example, CH 4 + 1 CH 2 −−→ CH 3 + CH 3 is only a key processing reaction in the upper atmosphere where the reactant 1 CH 2 is produced from the UV dissociation of CH 4 . Conversely, CH 3 + CH −−→ C 2 H 2 + H + H is only important in the lower atmosphere where a key alternate H-atom source (UV dissociation of H 2 ) does not occur.
Somewhat unintuitively, too much H in the lower atmosphere can lead to less HCN, as it is a reactant in the precursor sink reaction CH 3 + H + N 2 −−→ CH 4 + N 2 . For this reason, the precusor sink attenuation reaction H + H + N 2 −−→ H 2 + N 2 is also important for increasing HCN production in the lower atmosphere.
Next, HCN production occurs through 3 main channels, with the total initial percent contributions labeled in parentheses In the upper atmosphere, where partial pressures are low, UV radiation is the main dissociating agent; it generates reactants for Channels A-C. Channel B dominates the HCN production in this region, due to its comparatively high rate coefficient. Here, UV radiation is also responsible for breaking apart HCN into CN + H. Eddy diffusion mixes species including CN from the upper atmosphere into the lower atmosphere, where Channels A and D mainly take over HCN production. In the lower atmosphere, high partial pressures screen out UV radiation and increase the probability for GCR collisions, therefore GCRs become the main dissociating agent here. In this region, concentrations of CH 4 are high enough that the relative reaction rate of Channel A surpasses that of Channel B. These high concentrations of CH 4 also drastically increase the reaction rate of Channel D, which recycles CN back into HCN. This recycling process is the overall dominant channel to HCN, accounting for 36-46% of the total HCN in Titan's atmosphere. Channel C remains the fourth most important HCN source in both areas of the atmosphere, as although the rate coefficient of 2 N reacting with CH 3 is higher than that of 2 N reacting with CH 4 , there is a much higher concentration of CH 4 compared with CH 3 in all areas of the atmosphere.

Using CRAHCN
Due to the exceptional alignment of the HCN profiles from the CRAHCN19 and CRAHCN models, we suggest the CRAHCN19 network provides a lean, accurate, fast, and intuitively clear code to calculate the HCN abundance in Titan-like atmospheres. Without data from higher pressure and temperature planetary atmospheres, we cannot be certain that these same 19 reactions would suffice for other planetary environments of Titan-like composition. For this reason, we advocate using the full (104 reaction) CRAHCN network to simulate the production of HCN in N 2 -, CH 4 -, and H 2dominated atmospheres. We emphasize that this is a reduced network to accurately model HCN chemistry, rather than an extended network to cover the chemistry of a large range of species. For this reason, CRAHCN should only be used to simulate the production HCN.

CONCLUSIONS
In this paper, we calculate the production of HCN in Titan's atmosphere using a novel quantum chemistry and atmospheric modeling strategy. This strategy has two components: 1) we use quantum chemistry simulations to scan the entire field of possible reactions for a list of primary species relevant to N 2 -, CH 4 -, and H 2dominated atmospheres. We then calculate the rate coefficients for the uncovered reactions and construct a consistent reduced atmospheric hybrid chemical network (CRAHCN). This network contains experimental rate coefficients when available (32% of cases), but is predominantly composed of our calculated values using a consistent computational and theoretical method. 2) We pair CRAHCN with a chemical kinetic code called ChemKM to model the atmosphere of Titan. HCN has been observed at a range of altitudes in Titan's atmosphere by the Cassini spacecraft, making it an excellent testbed for validating chemical networks for HCN production in atmospheres.
We list the major conclusions of this work in bullet form below.
• CRAHCN contains 104 reactions, 33 of which are newly discovered in this work.
• Our calculated rate coefficients are accurate to within about an order of magnitude of experimental values, which is consistent with the uncertainties assigned in large-scale experimental data evaluations.
• In comparison with other widely used computational quantum methods, BHandHLYP/aug-cc-pVDZ is found to provide a reasonable balance of speed and acceptable accuracy for our large scale atmospheric study. Increasing the basis set to augcc-pVTZ did not improve the accuracy of calculations with respect to agreement with experimental values.
• The HCN profile from our fiducial model of Titan's atmosphere agrees very well with the Cassini observations, and is well in line with the three most recent Titan models in the literature.
• Only 19 reactions are responsible for the production and destruction of HCN in Titan's atmosphere. These reactions are sources and sinks of HCN, sources and sinks of the precursors to HCN, and a presursor sink attenuation reaction.
• There are 4 main channels to HCN production: • The first and second reactions dominate in the lower atmosphere, whereas the second, third, and fourth reactions dominate in the upper atmosphere. In the upper atmosphere, where partial pressures are low, relatively high rate coefficients tend to dictate the dominant reactions. In the lower atmosphere, the high CH 4 partial pressure increases the reaction rates for the first two reactions, which is the reason for their dominance here.
• The fourth dominant source of HCN in Titan's atmosphere is a reaction first discovered in our recent work (Pearce et al. 2019).
• HCN + 2 N −−→ N 2 + CH is the main sink for HCN. Conversely, HCN + hν −−→ CN + H is not an effective HCN sink, because it produces a CN molecule that reacts with CH 4 to form back into HCN.
• GCRs triple the total production of HCN in Titan's atmosphere, however they do not affect which reactions dominate HCN production and destruction.
Our work suggests that chemical networks of hundreds or thousands of reactions are not necessary to accurately simulate the production of HCN in N 2 -, CH 4 -, and H 2dominated atmospheres. Instead, using our novel strategy of exploring the entire field of possible reactions for a short list of primary atmospheric species has proven to be valuable at uncovering the dominant chemical pathways to producing HCN in Titan's atmosphere.
In upcoming work, we will use this strategy to expand CRAHCN to explore the production of HCN in the early Earth atmosphere, which, along with N 2 , CH 4 , and H 2 , is expected to have contained oxygen-based primary species such as CO 2 and H 2 O, as well as their dissociation fragments.

RATE COEFFICIENT CALCULATIONS
In Table A1, we display the recalculated two-body reaction rate coefficients from Pearce et al. (2019), with the modified vibrational partition function described in Section 2. Table A1. Re-calculated two-body reaction rate coefficients at 298 K from Pearce et al. (2019) with new vibrational partition function model. For these calculations, the vibrational partition functions are adjusted so that the zero of energy is at the first vibrational level, i.e., the zero-point level. Calculations are performed at the BHandHLYP/aug-cc-pVDZ level of theory. The presence or absence of an energy barrier in the rate-limiting step of the reaction is specified. The error factor is the multiplicative or divisional factor from the nearest experimental or suggested value.

Reaction equation
Forw./Rev. Barrier? k(298) calculated k(298) experimental Error factor a We remove the barrier from this calculation as experiments predict this reaction to be barrierless or nearly barrierless (Nesbitt et al. 1990).
See case study 1 for more details. b The existence of a barrier for this reaction is dependent on computational method. We introduce an experimental barrier of 6.3 kJ mol −1 (Takayanagi et al. 1999) to this calculation as at the BHandHLYP/aug-cc-pVDZ level of theory, we find this reaction to be barrierless. X Simulations had sporadic convergence beyond a C-N bond distance of 2.76Å. The rate coefficient is calculated with the variational transition state at this location, which has the highest ∆G. c The existence of a barrier for this reaction is dependent on computational method. We introduce an experimental barrier of 7.3 kJ mol −1 (Suzuki et al. 1993) to this calculation as at the BHandHLYP/aug-cc-pVDZ level of theory, we find this reaction to be barrierless. d We remove the barrier from the rate limiting step of this calculation, i.e. CH4 -CH · −−→ C2H5 · as experiments predict this reaction to be barrierless (Berman & Lin 1983;Blitz et al. 1997;Thiesemann et al. 1997).
In Table A2, we list the new two-body reaction rate coefficients calculated in this work at 298 K, along with any experimental values. Table A2. New two-body reaction rate coefficients calculated in this work at 298 K. Calculations are performed at the BHandHLYP/aug-cc-pVDZ level of theory. Reactions with rate coefficients slower than k = 10 −21 cm 3 s −1 are not included in this network. The presence or absence of an energy barrier in the rate-limiting step of the reaction is specified. The error factor is the multiplicative or divisional factor from the nearest experimental or suggested value.

Reaction equation
Forw./Rev. Barrier? k(298) calculated k(298) experimental Error factor a The existence of a barrier for this reaction is dependent on computational method. We introduce an experimental barrier of 8.3 kJ mol −1 (Bullock & Cooper 1972) to this calculation as at the BHandHLYP/aug-cc-pVDZ level of theory, we find this reaction to be barrierless. b Simulations had sporadic convergence beyond a H-N bond distance of 1.45Å. The rate coefficient is calculated with the variational transition state at this location, which has the highest ∆G. c The existence of a barrier for this reaction is dependent on computational method. We introduce an experimental barrier of 13.4 kJ mol −1 (ter Horst et al. 1996) to this calculation as at the BHandHLYP/aug-cc-pVDZ level of theory, we find this reaction to be barrierless. d Simulations did not converge beyond a N-C bond distance of 3.33Å. The rate coefficient is calculated with the variational transition state at this location, which has the highest ∆G.
In Table A3, we display the reduced Lindemann parameters for the three-body reactions calculated in this work, along with any experimental values. Table A3. Lindemann coefficients for the three-body reactions calculated in this work at 298 K. Calculations are performed at the BHandHLYP/aug-cc-pVDZ level of theory, and are valid within the 50-400 K temperature range. k∞ is the second-order rate coefficient in the high pressure limit with units cm 3 s −1 . k0 is the third-order rate coefficient in the low pressure limit with units cm 6 s −1 . These values fit into the pressure-dependent rate coefficient equation k = k 0 [M ]/k∞ 1+k 0 [M ]/k∞ k∞. Reactions with rate coefficients slower than k∞ = 10 −13 cm 3 s −1 are not included in this network. The error factor is the multiplicative or divisional factor from the nearest experimental or suggested value.

Reaction equation
k∞ (298) Table A4, we display the Lindemann coefficients for the three-body reactions in CRAHCN, along with their sources. When available, rate coefficients in CRAHCN are experimental values; However, the majority of rate coefficients have not been experimentally measured, and thus we use the consistently calculated values from this work. These three-body reactions, along with the one-and two-body reactions in Table A5, make up CRAHCN. Table A4. Lindemann coefficients for the three-body reactions in the consistent reduced atmospheric hybrid chemical network (CRAHCN), valid within the 50-400 K temperature range. Experimental values are used when available, and calculated rate coefficients from this work are used otherwise. k∞ is the second-order rate coefficient in the high pressure limit with units cm 3 s −1 . k0 is the third-order rate coefficient in the low pressure limit with units cm 6 s −1 . These values fit into the pressure-dependent rate coefficient equation In Table A5, we display the temperature-dependent Arrhenius parameters calculated for the one-and two-body reactions in CRAHCN. Table A5. Arrhenius coefficients for the one-and two-body reactions in the consistent reduced atmospheric hybrid chemical network (CRAHCN). Experimental values are used when available, and calculated rate coefficients from this work are used otherwise. For the reactions with barriers from this work, rate coefficients are calculated at 50, 100, 200, 298, and 400 K, and are fit to the modified Arrhenius expression k(T ) = α T 300 β e −γ/T . Barrierless reaction rate coefficients do not typically vary by more than a factor of 1-3 for temperatures between 50 and 400 K (Clary 1990;Li et al. 2008;Hase et al. 1987;Jasper et al. 2007;Daranlot et al. 2013). Intermediate molecules are labelled with a bullet, and are included to describe the precise reaction pathway for multi-step reactions. First-and second-order reactions with rate coefficients slower than k = 10 −21 cm 3 s −1 are not included in this network. First-order rate coefficients have units s −1 . Second-order rate coefficients have units cm 3 s −1 . CH3  Vibrationally excited molecules are labeled with a (ν) subscript.

EXPERIMENTAL DATA
In Table A6, we list all experimental rate coefficients for the reactions calculated in this work.

THEORETICAL CASE STUDIES
The following theoretical case studies provide additional details for some of the non-standard reactions in CRAHCN. For example, these reactions might have an excited intermediate, or may have a barrier that isn't detected by our chosen computational method.

CASE STUDY 1: H2CN + 4 N −−→ HCN + NH
One experiment has measured the rate coefficient for this reaction to have a value of k(298 K) = 4.4×10 −11 cm 3 s −1 (Nesbitt et al. 1990). No isotope effect was observed, which is consistent with a barrierless reaction. Using three temperature data points (200K, 298K, 363K), the authors suggest Arrhenius parameters indicative of a very small barrier. However, with the small number of data points, and the data uncertainties, complete temperature independence of this reaction would also fit these data points (Nesbitt et al. 1990). Nesbitt et al. (1990) suggest this reaction either proceeds through the N -CH 2 N complex, or via direct abstraction. We find the addition reaction forming the N -CH 2 N complex to have a large barrier at the BH/d level of theory (82 kJ mol −1 ). Furthermore, we find the addition reaction forming 3 H 2 CNN to efficiently decay into the N 2 and 3 CH 2 products. 3 H 2 CNN isomerization barriers proceeding to HCN + NH decay are too large to consider this pathway. At the BH/d and CC/t levels of theory, we find the direct abstraction reaction to have barriers of 15.4 and 23.0 kJ mol −1 , respectively. On the other hand, at the B3LYP/aug-cc-pVDZ and ωB/d levels of theory, we find the abstraction reaction to be barrierless and nearly barrierless (E 0 = 1.9 kJ mol −1 ), respectively.
Including the barriers, the rate coefficients for H 2 CN + 4 N −−→ HCN + NH at the BH/d and CC/t levels of theory are 9.5×10 −15 and 7.1×10 −14 cm 3 s −1 , respectively. These are several orders of magnitude lower than the experimentally measured value.
We find it likely that BH/d and CC/t calculate barriers when there should not be any. When removing the barriers from the calculation, the rate coefficients at the BH/d and CC/t levels of theory are 4.7×10 −12 and 5.1×10 −12 cm 3 s −1 , respectively. This reduces the discrepancy between experiment and calculation to factors of 9.
The rate coefficients at the ωB/d and B3LYP/aug-cc-pVDZ levels of theory are 3.4 and 9.4×10 −12 cm 3 s −1 , which are factors of 13 and 5 smaller than the experimental value. Given all these discrepancies, we recommend further experimental and theoretical analyses be carried out for this reaction.  Atakan et al. 1992). Baulch et al. (1992) reviewed both high-and low-temperature measurements and recommended a value of k = 3.0×10 −10 cm 3 s −1 .
The mechanism of the reaction is not previously well understood, however authors have suggested this reaction may pass through the CN 2 · intermediate (Atakan et al. 1992;Baulch et al. 1992). Our theoretical simulations suggest this is indeed the correct mechanism.
We calculate the rate coefficient of CN + N −−→ CN 2 at 298 K using the BH/d level of theory to be 3.4×10 −11 cm 3 s −1 . We find the decay of CN 2(ν=0) into N 2 + C to be slow (∼10 −41 s −1 ). This suggests that this reaction likely proceeds through vibrationally excited CN 2 , as is to be expected when two reactants combine to form a single product (Vallance 2017).
The experimental barrier for this reaction is about 8.3 kJ mol −1 (Bullock & Cooper 1972). This matches well with the MP4 theoretical barriers calculated by Yang et al. (1993) which range from 7.5-8.8 kJ mol −1 . We calculate a barrier for this reaction at the CC/t level of theory to be 6.7 kJ mol −1 .
We do not calculate a barrier for this reaction at the BH/d or ωB/d levels of theory. We calculate a barrierless rate coefficients at 298 K at to be 2.2 and 3.7 ×10 −11 cm 3 s −1 , respectively, which are over an order of magnitude greater than the nearest experimental value. This disagreement with experiment is due to the lack of barriers calculated at these levels of theory. For this reason, we introduce the experimental barrier of 8.3 kJ mol −1 (Bullock & Cooper 1972) to our BH/d calculation. This produces a rate coefficient at 298 K of 7.7×10 −13 cm 3 s −1 , which is within the experimental range.
Experiments generally agree on a barrier for this reaction of ∼16.7 kJ mol −1 (Che & Liu 1995). However, an ab initio theoretical study suggests a barrier of 13.4 kJ mol −1 provides a much better agreement between theory and experiment (ter Horst et al. 1996).
We calculate no barrier for this reaction at the BH/d level of theory. The barrierless rate coefficient we calculate at 298 K is 8.7×10 −13 cm 3 s −1 , which is a factor of 18 larger than the nearest experimental value. This discrepancy is due to the lack of barrier in our calculation when one should be present.
Similarly to ter Horst et al. (1996), when including a barrier of ∼16.7 kJ mol −1 , the calculated rate coefficient at 298 K is too low with respect to experimental values. Therefore, we include an experimental barrier of 13.4 kJ mol −1 in our calculation and obtain a rate coefficient of k(298 K) = 3.9×10 −15 cm 3 s −1 . This value is only a factor of 3 smaller than the nearest experimental value. CASE STUDY 5: CH4 + 1 CH2 −−→ C2H6 (ν) · −−→ CH3 + CH3 1 CH 2 insterts into the C-H bond of CH 4 to produce the C 2 H 6 intermediate, and subsequently dissociates into CH 3 + CH 3 (Langford et al. 1983;Ashfold et al. 1981;Tsang & Hampson 1986).
We calculate the decay of C 2 H 6 in the ground vibrational state into CH 3 + CH 3 to be slow (< 10 −37 s −1 ), which suggests this reaction proceeds through an excited vibrational state.
Previous theoretical studies find this reaction proceeds through three steps, first forming a CH 4 -CH complex, then C 2 H 5 , and finally, decay into C 2 H 4 + H (Ribeiro & Mebel 2015;Wang et al. 1997). At the CCSD(T)-F12/CBS//B3LYP/6-311G(d,p) + ZPE(B3LYP/6-311G(d,p)) level of theory, Ribeiro & Mebel (2015) found the first step of this reaction to be barrierless, and the second step to have a barrier 2.3 kJ mol −1 lower than the reactants. At the UMP2/6-31G(d,p) and UMP4/6-311++(2d,p)//UMP2/6-31G(d,p) levels of theory, Wang et al. (1997) found similar results, with a barrierless first step, and second step barriers 0.3 and 0.2 kJ mol −1 below the reactants, respectively. At the BH/d level of theory, we also find the first step to be barrierless; However, we find the CH 4 -CH · −−→ C 2 H 5 · step to have a barrier 11.5 kJ mol −1 above the reactants. This is similar to our result at the CC/t level of theory, where we find this step to have a barrier 7.7 kJ mol −1 above the reactants.
We also adjust our previously calculated reaction CH + H −−→ 1 CH 2 to be a three-body reaction with a low pressure rate coefficient of 4.4×10 −31 cm 6 s −1 . We assume that in low pressures the 1 CH 2(ν) molecule vibrationally decays back into CH + H. Table A8. Quantum Chemistry Data at the BHandHLYP/aug-cc-pVDZ level of theory. Cartesian coordinates are in angstroms. Energies are in kJ mol −1 . Ee is the electronic energy, ZPE is the zero point energy, and qx are the partition functions.