© ESO, 2012

1. Introduction

Based on the fact that HCNH⁺ was considered as an important ionospheric species in Titan’s atmosphere (Banaszkiewicz et al. 2000), Petrie (2001) made the hypothesis that hydrogen isocyanide (HNC) could be formed by the dissociative recombination of HCNH⁺ in Titan’s upper atmosphere where it might be detectable and might also play a part in the formation of more complex nitriles found on Titan. Recently, observations of Titan performed with the HIFI heterodyne submillimeter instrument aboard the Herschel Space Observatory (as part of the guaranteed time key program “Water and related chemistry in the solar system” (HssO), Hartogh et al. 2009), allowed the first detection of HNC in Titan’s atmosphere through the measurement of its emission line from the J = 6  →  5 rotational transition at 543.897 GHz (Moreno et al. 2010). Their preliminary analysis suggests that the bulk of this emission must originate at altitudes above 300 km. However, the observations cannot strictly establish a HNC vertical profile (Moreno et al. 2011). Petrie (2001) argued that HNC was likely to be formed almost entirely by an ion-molecule mechanism and that neutral formation pathways were not viable. As a consequence, Petrie (2001) suggested that the concentration profile for HNC as a function of altitude would follow the typical profile for a polyatomic ion rather than a neutral molecule. So, according to Petrie (2001), HNC might be located mainly in the ionosphere with a peak of abundance around 1200 km.

Since the neutral composition of the atmosphere is a prerequisite for ionospheric models, it is of prime importance to study carefully the production of neutral species production in addition to loss processes for HNC and HCN, since HCN is strongly related to HNC. Since HNC was not considered in previous photochemical models, it is necessary to build a new chemical scheme devoted to HNC chemistry and its putative interactions with other species. The chemistry of nitrogen compounds at low temperatures (between 100 and 200 K), either in cold planetary atmospheres or in protoplanetary disks, is not well known (Hébrard et al. 2006; Vasyunin et al. 2008; Hébrard et al. 2009). Therefore, it is also critical to evaluate the uncertainties attached to the reviewed reaction rates.

In the present paper, we investigate the production of HNC via neutral reactions. A careful investigation of neutral production and loss processes for both HNC and HCN has been carried out. A 1D photochemical model is used to infer the abundance of HNC as a function of altitude. Our study includes an uncertainty propagation study and a sensitivity analysis to determine the key reactions of our chemical scheme. We briefly present our photochemical model in Sect. 2. The methodology we have adopted to build the chemical scheme for HNC and HCN is presented in detail in Sect. 3. In Sect. 4, we compare the computed abundance profiles of HNC and HCN with the available observations. We present the main reactions for the neutral production and loss of these compounds in the current model. A local sensitivity study is carried out in Sect. 5 to highlight the importance of the HNC isomerization reaction. An uncertainty propagation study is presented in Sect. 6 and a global sensitivity analysis is performed subsequently in Sect. 7 to determine the key reactions which are important to study in priority to improve the predictivity of our model. The main conclusions of our work are summarized in Sect. 8. Some selected reactions (reactions that are important for the production of HNC or HCN and reactions that contribute significantly to the uncertainties on their abundances) are commented in Sect. A.

2. Photochemical model

The photochemical model is derived from the Hébrard et al. (2007) model with a modification of the numerical solver presented in Dobrijevic et al. (2010a). Instead of using a classical Crank-Nicholson method, we now use the ODEPACK library, which implements Hindmarsh’s solvers for ordinary differential equations (Hindmarsh 1983). Our photochemical 1D model uses a constant background atmosphere with constant boundary conditions. Atmospheric parameter inputs (T, P, n) were taken from Yelle et al. (1997) recommended engineering model. We use a non-uniform grid of altitude with 125 levels from the ground to 1300 km. Two consecutive levels are separated by a distance smaller than H(z)/5, where H(z) is the atmospheric scale height at altitude z. A zero flux was assumed as an upper boundary condition for most of the species, except for atomic hydrogen H and molecular hydrogen H₂, which were allowed to escape with velocities following Jean’s thermal escape mechanism, and for water H₂O, which exhibits an external influx equal to 5  ×  10⁶ cm^-2 s^-1 to account for the water influx arising from micrometeorites initiating oxygen chemistry in Titan’s atmosphere (Feuchtgruber et al. 1997). At the lower boundary, methane CH₄ abundances was set to its tropopause abundance 1.41  ×  10^-2. The abundance of carbon monoxide CO at the surface was assumed to be 5. 2  ×  10^-5 on the basis of the (Gurwell & Muhleman 2000) high resolution ground-based interferometric observations. Molecular hydrogen H₂ abundance was assumed to be equal to 1.1  ×  10^-3 following (Samuelson et al. 1997). Calculations are performed with a solar zenith angle of 50° to account for diurnally averaged conditions at the equator. The eddy diffusion coefficient K(z) is a free parameter of 1D photochemical models which is not well constrained. We use the value derived by Hörst et al. (2008) as a mean profile among all the different profiles that have been published so far (see Fig. 1).

Fig. 1 Eddy diffusion profiles from various photochemical models (Yung et al. 1984; Strobel et al. 1992; Toublanc et al. 1995; Lara et al. 1996; Wilson & Atreya 2004; Vinatier et al. 2007; Lavvas et al. 2008; Hörst et al. 2008; Krasnopolsky 2009).

3. Chemical scheme

The basis of the chemical scheme is presented in Hébrard et al. (2006) and Hébrard et al. (2009). In the present work, many rate constants have been updated and numerous reactions have been added. Our new chemical scheme includes 137 compounds and 788 reactions (91 photodissocation processes, 2 dissociation processes of N₂ by cosmic rays, 694 bimolecular reactions and 94 termolecular reactions). The complete list of the reactions is available upon request or can be downloaded from the KInetic Database for Astrochemsity (KIDA, http://kida.obs.u-bordeaux1.fr). In the present paper, we only present reactions which are essential for the study of HCN and HNC, either because they are important for their neutral production and/or loss or because they contribute significantly to the uncertainties on their abundances.

The methodology we have adopted to build a chemical scheme for HNC and HNC is the following. As a first step, all the reactions related to HCN found in the chemical schemes published by Hébrard et al. (2006) and other photochemical models were checked and updated, by including in addition new reactions found in the literature. For all these reactions, we investigated their efficiency to produce HNC. In a second step, we completed the HNC scheme by introducing all neutral reactions producing of destroying HNC that we could find in the literature. Then, based on the relative abundance of the various chemical species computed by our photochemical model, we estimated what reactions currently missing in our chemical scheme might be important by systematically evaluating the cross-reactions involving HCN or HNC between the main species of Titan’s atmosphere: C, H, N(⁴S), N(²D), CH, CH₃, CH₄, C₂H₂, C₂H₄, NH, CN, CH₂NH, H₂CN, etc. Finally, we also estimated what reactions might be relevant to improve the chemical scheme, in particular to avoid the artificial formation of important sink species (species which are efficiently produced but not destroyed in the model). At the end, the chemical scheme consists of 56 reactions with HCN and 22 reactions with HNC (including photolysis processes). Among these reactions, we have estimated in the present work the rate constants of 12 reactions for HCN and of 14 reactions for HNC. This illustrates the lack of information regarding HNC reactivity in the literature.

3.1. Estimation of reaction rates for new processes

When introducing reactions with unknown reaction rates (without experimental measurements and/or theoretical calculations to rely upon) in the temperature range of interest (T ∈ [100−200]  K), one major concern is to make a reasonble estimation of their reaction rates. In the following, we explain how we have estimated the rate constants and branching ratios for various reactions. For the evaluation of chemical rate constants between 150 K and 200 K the presence of an energetic barrier is critical. When no information was available, the presence and the values of any energetic barriers to the entrance valley for the important reactions (in terms of production and/or loss rate) were calculated at the M06-2X/cc-pVTZ level using the Gaussian09 software package (Frisch et al. 2009) except for the H + H₂CN, N(⁴S) + H₂CN and N(²D) + HCN reactions for which calculations were performed also at MRCI+Q/vqz level using the Molpro software package (Werner et al. 2010). For the other reactions, identified as being less important, the presence or the absence of a barrier was deduced from general considerations. When a reaction was thought to proceed through direct abstraction, its rate constant has been estimated by comparison with known similar reactions. When a reaction was thought to proceed through addition (as is the case with most of the radical-radical reactions), in the absence of theoretical calculations, we have considered that there was no barrier when the ground state of the adduct arises from pairing up electrons on the two radicals reactants, whereas the surface was likely be repulsive if all the electrons remain unpaired. Therefore, doublet + doublet reactions were considered to have no barrier for the singlet surface but a barrier for the triplet surfaces (which is in good agreement with experimental and theoretical results for H + alkyl or alkyl + alkyl reactions for example, Harding et al. 2005; Klippenstein et al. 2006).

When no energetic barrier was found to be present, the value of the rate constant was estimated using long-range forces, mainly through dispersion interactions (Stoecklin & Clary 1992; Georgievskii & Klippenstein 2005), and by taking into account the electronic degeneracy γ_el. The electronic degeneracy factors were calculated by applying the spin and orbital correlation rules to the potential energy surfaces that correlate the separated reactants with the separated products. Using a capture rate, k_capture, generally overestimates the rate constant by a uncertainty factor F_capture = 3 at most (Georgievskii & Klippenstein 2005) except for specific mechanisms for which the potential coexistence of a van der Waals complex and a submerged barrier could play an important role and lead to low rate constant values at 300 K (like for OH + alkenes or CN + alkenes reactions). Accordingly, the nominal rate constant used in the model for such reactions is k = γ_el  ×  ${\mathit{k}}_{\mathrm{capture}}\mathit{/}\sqrt{{\mathit{F}}_{\mathrm{capture}}}$ to which is associated an uncertainty factor $\mathit{F}\mathrm{=}\sqrt{{\mathit{F}}_{\mathrm{capture}}}$. As a result, the minimum and maximum values expected for this effective rate constant are k_min = γ_el  ×  k_capture/F_capture and k_max = γ_el  ×  k_capture, respectively. Values for the effective uncertainty factor F were estimated mainly by comparison with similar reactions and also by taking into account the uncertainties on electronic degeneracy when no ab-initio calculations were available. Branching ratios were estimated from ab-initio calculations except for the key H + H₂CN reaction for which we performed statistical calculations of the microcanonical rate constants of the various steps of the mechanism (Bergeat et al. 2009).

3.2. Photolysis processes

In our model we consider that the photolysis of HNC is similar to the photolysis of HCN (using the identical absorption cross sections, dissociative thresholds and quantum yields). Compared to the Hébrard et al. (2009) model, we add several photolysis processes which are presented in Table 1.

4. Photochemical model results

4.1. Comparison with observations

4.1.1. HNC observations

Very recently, Moreno et al. (2011) reported the first identification of HNC in Titan’s atmosphere from observations using the HIFI instrument on the Herschel¹ Space Observatory. The column density of HNC inferred from these observations is in the range (0.6−1.5)  ×  10¹³ cm^-2 for altitudes between 400 and 1000 km, but the authors noticed that it was not possible to constrain the vertical profile of HNC from these data. Several constant profiles of HNC give a satisfactory agreement with the observations depending on the mixing ratio and the altitude cut-off. These profiles are presented in Fig. 2 and are compared with our model. The column density we obtain is 3.4  ×  10¹³ cm^-2 above 500 km whereas the value derived by Moreno et al. (2011) at this altitude is 1.2  ×  10¹³ cm^-2 (about 3 times lower). We will see in the following that this discrepancy can be simply explained by our poor knowledge of some key reactions.

Fig. 2 Mixing ratio of HNC in the upper atmosphere of Titan. Nominal model (solid line) and the different acceptable profiles derived from recent Herschel observations (Moreno et al. 2011) (grey dotted lines). Local sensitivity analysis for the isomerization reaction H + HNC  →  HCN + H is also illustrated by changing its rate by a factor of 10 or 100.

4.1.2. HCN observations

The neutral composition of Titan’s upper atmosphere between 1000 and 1100 km has been inferred from the interpretation of Cassini Ion and Neutral Mass Spectrometer (INMS) measurements by Magee et al. (2009). In particular, the global average mixing ratio of HCN at 1050 km is (2.44  ±  0.10)  ×  10^-4. Previous interpretation of INMS data from Vuitton et al. (2007) gave a similar abundance of HCN of 2.0  ×  10^-4 at 1100 km, with an uncertainty factor of 2−3. Geballe et al. (2003) detected HCN emission features in a high resolution spectrum near 3 μm acquired at the Keck II telescope. These data were re-analyzed by Yelle & Griffith (2003) with a model for fluorescence of sunlight in the ν3 band of HCN and by Kim et al. (2005) with an updated model. Their results were in agreement in the upper atmosphere. The HCN mixing ratio is from (1−3)  ×  10^-3 around 1000 km. The density of HCN between 600 and 1000 km has been also inferred by Shemansky et al. (2005) from observations of stellar occultations by the atmosphere of Titan using the Cassini Ultraviolet Imaging Spectrometer (UVIS). The HCN mixing ratio at 1000 km is about 4  ×  10^-3. Recently, Adriani et al. (2011) used the limb observations of the Visual and Infrared Mapping Spectrometer (VIMS) onboard the Cassini spacecraft to retrieve vertical profiles of HCN from its 3 μm non-LTE emission in the region from 600 to 1100 km altitude at daytime. The mixing ratio of HCN is about (5.5  ±  1.5)  ×  10^-3 at 1050 km. This result is in agreement with the upper value of the HCN profile derived by Yelle & Griffith (2003) and Kim et al. (2005). HCN has been also detected in Titan’s upper atmosphere in the ultraviolet by the UltraViolet Spectrometer (UVS) instrument aboard Voyager 1 (Vervack et al. 2004). The mixing ratio of HCN is around 10^-4 at 500 km (about 20 times lower than other observations) and then increases with altitude leading to a value in agreement with all the other observations (but with large uncertainties).

It is difficult to compare all these observations taken at different times, using different techniques and corresponding to different spatial resolutions. In addition, these observations are more or less model-dependent and a direct comparison with our nominal profile is not straightforward. Due to the inconsistency of the published observational data (error bars do not overlap), we can conclude that our nominal neutral model, which does not include ions, is roughly consistent with these observations.

Fig. 3 Mixing ratio of HCN in the upper atmosphere of Titan. Model (solid line) and main observations in the upper atmosphere (see text for references).

4.1.3. Sensitivity to eddy diffusion

If restricted above 1000 km, the mixing ratio of HNC derived from observation of Moreno et al. (2011) is about 6  ×  10^-5 and the mixing ratio of HCN is in the range (0.2−6)  ×  10^-3 according the various observations. This leads to a HNC/HCN ratio of about 0.01−0.3. In our current model, the HNC/HCN ratio is about 0.1. At 500 km, the HNC/HCN ratio from the various observations is around 5  ×  10^-3 (assuming a constant mixing ratio of HNC above 500 km). In our model, the ratio is 2  ×  10^-5, more than 100 times lower.

Since the eddy diffusion coefficient is not well constrained in the atmosphere of Titan, we have tested other K(z) profiles in order to test the sensitivity of this ratio to transport. For instance, in the case of the eddy diffusion profile K(z) derived by Hörst et al. (2008), the methane homopause is located around z_h = 800 km and K_h = 3  ×  10⁷ cm² s^-1. In the case of the eddy diffusion profile K(z) obtained by Strobel et al. (1992), the methane homopause is around z_h = 1100 km and K_h = 10⁹ cm² s^-1. Figure 4 shows that the HNC/HCN ratio is not sensitive to K(z) down to 700 km but is different by a factor of 100 at 500 km. At 500 km, the discrepancy between our model using the K(z) of Strobel et al. (1992) and observations is more pronounced since the ratio is only 3  ×  10^-7.

HNC and HCN mixing ratios for the two eddy profiles are also presented in Fig. 4. Their abundances differ by a factor of 3 above 1000 km. As explained in the following section, both HNC and HCN are mainly produced from the reaction of H₂CN, which is produced by the reaction between N(⁴S) and CH₃. So, these differences are directly linked to the mixing ratio of CH₃ which depends on the photolysis and the molecular diffusion of CH₄ above the homopause.

Consequently, the sensitivity of HNC and HCN to K(z) is lower than the uncertainties on the model (see also Figs. 8 and 9) and can not be used to constrain the eddy diffusion coefficient profile in the upper atmosphere of Titan.

Fig. 4 Mixing ratios of HNC and HCN and HNC/HCN ratio as a function of altitude for two eddy diffusion coefficients.

4.2. Main production and loss processes for HNC

The main production and loss processes for HNC (for the current chemical scheme) are presented in Fig. 5 and are listed in Table 2. Here we summarize the main processes leading to the production and loss of HNC from the primary radicals.

4.2.1. Production

At 1300 km, HNC is produced from:

At 1000 km, HNC comes mainly from the reaction:

H₂CN is mainly produced by two reactions:

while CH₂NH comes from:

and NH comes from:

At 600 km, the scheme for the production of HNC is a little bit different; HNC comes fully from the reaction:

but H₂CN is produced from:

and CH₂NH comes from:

As a conclusion, H₂CN and CH₂NH are key species in the production of HNC in the upper atmosphere of Titan.

Fig. 5 Reaction rate profiles for most important reactions leading to the production and loss of HNC (see Table 2 for reactions).

4.2.2. Loss

At 1300 km, HNC is mainly destroyed by its reaction with N(²D):

At 1000 km, the isomerization of HNC into HCN begins to contribute to the HNC loss as well:

At 600 km, HNC is fully destroyed through its isomerization into HCN:

4.3. Main production and loss processes for HCN

The main production and loss processes for HCN are presented in Fig. 6 and are listed in Table 3. Here we summarize the main processes leading to the production and loss of HCN from the primary radicals.

4.3.1. Production

At 1300 km, several reactions contribute to the production of HCN:

At 1000 km, the situation is simpler, HCN comes from:

In a similar manner to HNC, the production of HCN is strongly related to the production of H₂CN. Again, the reaction between N(⁴S) and CH₃ is important since it produces both H₂CN and HCN.

At 600 km, several reactions contribute equally to the production of HCN. Many of them involve CN radicals.

Where C₂H₃CN comes from:

And CN comes from:

As a conclusion, H₂CN (z  >  1000 km) and CN (z  <  1000 km) are important intermediate compounds for the production of HCN in the upper atmosphere of Titan.

4.3.2. Loss

At 1300 km, HCN is mainly destroyed by its reaction with N(²D) and by its photolysis:

At 1000 km, the situation is quite similar:

At 600 km, HCN is partly destroyed through its reaction with C₃N radical:

Fig. 6 Reaction rate profiles for most important reactions leading to the production and loss of HCN (see Table 3 for reactions).

Remark.

While studying the relative importance of the individual contribution of each reaction to the production and loss rates of a given compound, it is important to keep in mind the following fact: due to the large current uncertainties in reaction rates, the relative importance of each reaction might change drastically as we go along in the improvement of reaction rates accuracy at low temperature. As a consequence, the set of reactions given here for our current nominal model might not be the set of reactions that would dominate had some of the initial rate constants been different in the production or loss of HNC and HCN (see Sect. 6 for an illustration of the impact of uncertainties on reaction rates in model outputs).

5. Local sensitivity analysis: a critical reaction study

Figure 5 shows that the most important loss reaction for HNC around 500 km is the isomerization reaction H + HNC  →  HCN + H (see also Table 3). The reaction rate of this process is not very well known at low temperature (≈150 K) with an uncertainty factor that we estimate to be about 10 (see Appendix A). So, this reaction is clearly critical for the abundance of HNC in our model. In order to pinpoint the importance of this reaction, we performed a local sensitivity analysis, which consists of changing the rate of this reaction by a factor of 10 whilst keeping all other rates unchanged. It is important to note that this kind of study gives only incomplete information about the importance of this reaction. A global sensitivity analysis is required to study how the uncertainty of this reaction rate propagates in the model through the set of strongly coupled and non-linear differential equations (Dobrijevic et al. 2010b). Results are presented in Fig. 2. Whilst there is only a factor of 2 between the extreme HNC mole fractions at 1300 km, the difference reaches a factor of 10⁵ at 500 km. As a consequence, the poor knowledge of the rate of this reaction limits strongly the accuracy of photochemical models.

6. Uncertainty propagation study

6.1. Method

The methodology used to study the propagation of uncertainties in the model is described in Hébrard et al. (2007) and Hébrard et al. (2009). There are many sources of uncertainty in a 1D photochemical model. In the present study, we focus exclusively on its chemical sources through the uncertainties of the photodissociation and reaction rates. These uncertainties originate in their experimental or theoretical determination, and are generally quantified by a standard deviation or a relative uncertainty. Because of the profoundly non-linear nature of the photochemical model and the potentially large uncertainties displayed by many parameters, a linear uncertainty propagation is not expected to produce valid results. Propagation of distributions by Monte Carlo sampling is better adapted to such problems (BIPM et al. 2008, 2006). Due to the positivity constraint on these properties, their distributions are modeled by lognormal probability density functions: $\mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}\mathit{x\sigma }}\exp \left(\mathrm{-}\frac{\mathrm{(}\ln \mathit{x}\mathrm{-}\mathit{\mu }{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}\right)$(1)with μ = lnk(T), the logarithm of the nominal value of the reaction rate at temperature T, and σ = lnF(T), the logarithm of the geometric standard uncertainty F(T) of the lognormal distribution. With these notations, the 67% confidence interval for a reaction rate at a given temperature is given as [k(T)/F(T),  k(T) × F(T)].

Estimation of the uncertainty factor F(T) of a reaction rate k(T) at any given temperature follows an expression adapted from KIDA (Wakelam et al. 2012): $\mathit{F}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\mathit{F}\mathrm{(}\mathrm{300} \mathrm{K}\mathrm{)}\exp \left|\mathit{g}\left(\frac{\mathrm{1}}{\mathit{T}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{300}}\right)\right|$(2)where F(300   K) is the uncertainty in the rate constant k(T) at T = 300 K and g is the “uncertainty-extrapolating” coefficient defined for use with F(300   K) in the above expression to obtain the rate constant uncertainty F(T) at different temperatures. These assigned uncertainty factors F(300   K) and g are evaluated to construct the appropriate uncertainty factor, F(T), following an approach based on the fact that rate constants are almost always known with a minimum uncertainty at room temperature. The knowledge of both F(300   K) and g parameters allows to quantify the temperature-dependent uncertainties carried by each reaction rate present in the standard sets of reaction rates in a temperature range adequate for Titan’s atmosphere.

Most of the reaction rate coefficients and their associated uncertainty factors used in the present study are extracted from our previous reviews (Hébrard et al. 2006; Hébrard et al. 2009). Reaction rates and uncertainties for new processes are estimated according to the methodology introduced in Sect. 3. The uncertainty factor is set to 1.2 for all the photodissociation rates for simplicity (see Peng et al. 2012, for a valuable discussion about uncertainties on photodissociation rates). We perform 1000 runs to have statistically significant results. The integration time for each run is set to 10¹¹ s for simplicity and to limit the computation time. This time is sufficient to reach a steady state in the upper atmosphere. For instance, the most important relative variation of the HNC mole fractions is lower than 1.0  ×  10^-2 for all the runs with very few exceptions.

Figure 7 presents an example of a rate constant histogram generated by our Monte-Carlo procedure for a reaction rate constant estimated by capture theory. As explained in Sect. 3, this kind of rate constant is in practice $\begin{array}{ccc}\mathit{k}\mathrm{=}{\mathit{\gamma }}_{\mathrm{el}}\mathrm{\times }{\mathit{k}}_{\mathrm{capture}}\mathit{/}\sqrt{{\mathit{F}}_{\mathrm{capture}}}& & \end{array}$and its uncertainty factor $\mathit{F}\mathrm{=}\sqrt{{\mathit{F}}_{\mathrm{capture}}}$. Our simulation is in good agreement with what was expected and we see that only a limited fraction of runs ( < 15%) give rates greater than the k_capture value.

Fig. 7 Histograms of reaction rates obtained after 1000 runs for reaction N(4S) + C2H  →  C2N + H. The reaction rate estimated for this reaction is kcapture = 6.0  ×  10-11 and its uncertainty factor is Fcapture = 16.

6.2. Results: HCN and HNC profiles

The type of abundance distribution depends on the compound and can vary with altitude: distributions are not always normal or log-normal. In this case, quantiles are useful measures to represent the distributions. Figures 8 and 9 show the 5th and 15th of the 20-quantiles and the 1st and 19th of the 20-quantiles which give the intervals containing respectively 50% and 90% of the profiles.

The HCN and HNC profiles obtained from the uncertainty propagation study are presented in Fig. 8. We obtain large uncertainties for HNC, especially below 900 km, whereas uncertainties are quite reasonable for HCN throughout the upper atmosphere. At 1300 km, the mean value of the HNC mole fraction is 10^-4 and 50% of the profiles lie between 5.3  ×  10^-5 and 2.0  ×  10^-4. For HCN, the mean mole fraction is 1.2  ×  10^-3 and 50% of the profiles lie between 8.0  ×  10^-4 and 1.6  ×  10^-3.

Fig. 8 Abundance profiles of HNC and HCN obtained after 1000 runs. Black solid line: initial profile. Black dotted line: median profile obtained from the uncertainty propagation study. Black dashed-dotted lines: 5th and 15th 20-quantiles of the distribution. Black long-dashed lines: 1st and 19th 20-quantiles of the distribution.

The HNC/HCN ratio as a function of altitude is given in Fig. 9. The determination of this ratio is important with respect to the detection of HCN by the Cassini Ion Neutral Mass Spectrometer (INMS) (Magee et al. 2009). Since HNC and HCN have the same mass, it is likely that both species contributed significantly to the same peak in this instrument. We see in Fig. 9 that the ratio HNC/HCN can statistically reach a value as high as 0.5. More precisely, 90% of the HNC/HCN profiles are within the range [0.015, 0.5] and 50% in the range [0.04, 0.17]. By comparison with observations, Moreno et al. (2011) estimated that this ratio is about 0.3 (based on the abundance of HCN retrieved by INMS).

Fig. 9 Abundance profiles of HNC/HCN ratios obtained after 1000 runs. Black solid line: initial profile. Black dotted line: median profile obtained from the uncertainty propagation study. Black dashed-dotted lines: 5th and 15th 20-quantiles of the distribution. Black long-dashed lines: 1st and 19th 20-quantiles of the distribution.

The main conclusion of these results is that the current accuracy of photochemical models for HNC is very poor, especially below 800 km. Theoretical and experimental studies are required to improve the situation. We see in the next section that a global sensitivity analysis can pinpoint the key reactions that are responsible for these huge uncertainties.

7. Global sensitivity analysis: determination of key reactions

The technique we use is based on the computation of the Rank Correlation Coefficients (RCCs) between rate constants and abundances at different altitudes. It has been previously used and described in Carrasco et al. (2007); Dobrijevic et al. (2008); Hébrard et al. (2009). The greater the absolute value of a RCC is, the more important the contribution a given reaction rate has on the abundance uncertainty for a given species. As a consequence, all reactions with a high absolute value of RCC should be studied in priority to improve the accuracy of the model concerning HNC and HCN in the upper atmosphere of Titan. The power of this technique, based on the uncertainty propagation study and the knowledge of RCCs, to improve photochemical models has been demonstrated for Titan by Hébrard et al. (2009) and for Neptune by Dobrijevic et al. (2010a).

Among the few key reactions we give here, H₂CN + H, HNC + H and N(²D) + HNC/HCN are clearly the most important ones for HNC and HCN in the upper atmosphere of Titan. In Appendix A, we discuss in more details our current knowledge about these reactions.

8. Conclusions

Following the prediction of Petrie (2001), Moreno et al. (2011) concluded that a purely ionospheric source may be quantitatively viable for HNC, provided that the protonation rates are not too high. Here we show that a purely neutral source is efficient enough to produce HNC and HCN in the upper atmosphere of Titan in agreement with current observations. Our study does not mean that ion-molecule mechanisms are not relevant to produce HNC but it shows that neutral reactions are competitive processes in the upper atmosphere and in the middle atmosphere as well.

The aim of the present work is to investigate as exhaustively as possible the neutral chemistry of HNC and HCN in order to construct a chemical scheme as complete as possible, to evaluate the uncertainties on the results that originate from the ones attached to the neutral reactions and to determine the key reactions that should be studied in priority to improve the model’s precision regarding neutral chemistry. We show that the precision on the HCN abundance predicted by photochemical models is currently strongly limited by the poor knowledge of some reaction rates, especially the isomerization reaction H + HNC  →  HCN + H. Further studies of the reactivity of H₂CN with H and the reactivity of N(²D) with HNC and HCN are also very important.

Our results are a basis for studying the production and loss of HNC and HCN in the ionosphere of Titan. For instance, Krasnopolsky (2009) stated that 24% of the production of HCN comes from ion reactions, while ion reactions account only for 17% of the loss of HCN. More recently, Plessis et al. (2012) found that the production of HCN and HNC are about 5.9  ±  2.5 molecule cm^-3 s^-1 and that 65% of this production comes from the dissociative recombination reaction HCNH⁺ + e⁻. So, we can expect that the production of HCN and HNC should be greater than the production we find in the present study taking into account only neutral chemistry. However, we can also predict that uncertainties in ion-reaction rates might increase the uncertainty on HNC and HCN mole fractions obtained by the models. In conclusion, a coupled neutral and ion model is required to investigate the impact of ion chemistry on these compounds, to determine the uncertainties on the computed abundances and to pinpoint the key reactions involving ions.

Acknowledgments

The authors acknowledge partial financial support from the Observatoire Aquitain des Sciences de l’Univers (OASU), from the CNRS “Programme National de Planétologie” and from the CNRS interdisciplinary program “Environnements Planétaires et Origines de la Vie (EPOV)”.

References

Appendix A: Comments on selected reactions

A.1. N(⁴S) + CH₃

We have shown that the reaction N(⁴S) + CH₃  →  products is of prime importance in the production of H₂CN, which appears to be a key compound for initiating both HCN and HNC chemistries in Titan’s upper atmosphere. For this reason, we give in the following details our recommendation for the rate constant and the branching ratios.

The N(⁴S) + CH₃ reaction has been studied experimentally down to 200 K (Marston et al. 1989a) showing a high value for the rate constant at 300 K (8.5  ±  2.0)  ×  10^-11 cm³ molecule^-1 s^-1 with a complex negative temperature dependence. The authors propose two expressions for the global rate constant, k(200 − 423   K)  =  6.2  ×  10^-11 + 2.2  ×  10^-9exp( − 1250/T) cm³ molecule^-1 s^-1 or k(T) = 4.3  ×  10^-10exp( − 420/T) cm³ molecule^-1 s^-1. The second one deviates from the experimental results at low temperature and may underestimate the rate constant at temperatures representative of Titan’s atmosphere. The possible exothermic exit channels are: $\begin{array}{ccc}\begin{array}{c}\\ & & \\ & & \\ & & \\ & & \end{array}& & \end{array}$Product branching ratios have been obtained for the N(⁴S) + CH₃ and N(⁴S) + CD₃ reactions (Marston et al. 1989b) leading mainly to H₂CN + H formation (85 − 100%) with some HCN formation (0 − 15%). These authors suggest that HCN formation is associated with H₂, however HCN + H₂ production is spin-forbidden and needs intersystem crossing to occur. Additionally recent ab-initio calculations (Cimas & Largo 2006) found almost 100% of H₂CN production in good agreement with previous calculations (Nguyen et al. 1996). As H₂CN may have enough internal energy (153 kJ mol^-1) to overcome the dissociation barrier for C-H dissociation (130 kJ mol^-1), some HCN + H + H may be produced (Nguyen et al. 1996). The HCN obtained by Marston et al. (1989a) may also come from reaction of H₂CN with atomic nitrogen (used in excess in their experiment) through H₂CN + N(⁴S)  →  HCN + NH or more likely through H₂CN + N(⁴S)  →  CH₂ + N₂ followed by N(⁴S) + CH₂  →  HCN + H. We choose to recommend a rate constant value of 6.2  ×  10^-11 cm³ molecule^-1 s^-1 between 150 K and 200 K with a branching ratio equal to 90% for H₂CN + H formation and 10% for HCN + H + H. Thus:

A.2. HCN  →  HNC isomerization

A specific problem in HCN/HNC formation is the possibility of isomerization. The HCN  →  HNC isomerization barrier is calculated equal to 186 kJ mol^-1, and the HNC  →  HCN isomerization barrier is calculated equal to 124 kJ mol^-1 at the RCCST(T)/cc-pVTZ level (DePrince III & Mazziotti 2008). Some reactions producing HCN and all reactions producing HNC are highly exothermic. In Titan’s atmosphere, relaxation occurs through collisional stabilization. The collisional stabilization time may be estimated (Forst 2003) as equal to $\frac{\mathrm{1}}{\mathrm{3}\mathrm{\times }{\mathrm{10}}^{-10}\mathrm{\times }\mathrm{\left[}\mathit{M}\mathrm{\right]}}\mathrm{\approx }{\mathrm{10}}^{-7}$ s at 1 Torr, much greater than the interconversion time-scale estimated as  < 10^-13 s (Herbst et al. 2000). Thus, as relaxation slowly occurs, isomerization leads to equilibrated isomeric abundances at each internal energy. The final balance is determined at or near the effective barrier for isomerization. As the available energy of the exothermic reactions producing HCN and HNC can be distributed among the translational and vibrational-rotational modes of both products, a fraction of HCN/HNC molecules will possess enough internal energy to overcome the isomerization barrier. To a first approximation we assumed roughly statistical distribution of energy in the fragments and consider that all of the HCN or HNC produced above the isomerization barrier will lead to equal amounts of HCN and HNC (Herbst et al. 2000).

The HNC + H  →  HCN + H reaction (ΔH_r(298 K) =  − 53 kJ mol^-1) has been calculated to possess a low barrier. The first estimation of the barrier was performed by Talbi & Ellinger (1996) at a relatively high level of theory (MP4/6-311++(3df,2p)//MP3/6-311++G(d,p)) with a best estimate value equal to 17.6  ±  4 kJ mol^-1. Sumathi & Nguyen (1998) found a barrier equal to 13.8 kJ mol^-1 at the CCSD(T)/6-311++(3df,3pd)//CCSD(T)/6-311++G(d,p) level and more recently Petrie (2002) found a value equal to 8.0 kJ mol^-1 at the CBS/RAD or B3LYP/6-311G** level and 12.9 kJ mol^-1 at the QCISD/6-311G* level. With estimated values between 8.0 kJ mol^-1 (960 K) and 18.0 kJ mol^-1 (2160 K), the rate will be low, but not negligible, for relaxed HNC. Moreover, as HNC is produced mainly with high internal energy there is a possibility of rate enhancement. For relaxed HNC we recommend the rate constant calculated by Sumathi & Nguyen (1998) with an average value of the various calculations for the barrier including Petrie (2002):

A.3. N(²D) + CH₄

The N(²D) + CH₄ reaction has been studied experimentally (Takayanagi et al. 1999; Umemoto et al. 1998) and theoretically (Ouk et al. 2011; Takayanagi & Kurosaki 1999; Takayanagi et al. 1999) and a review has been performed by Herron (1999). Theoretical calculations suggest two pathways for this reaction, direct H atom abstraction and N(²D) insertion in one C-H bond, both mechanism presenting a barrier in the entrance valley (Ouk et al. 2011; Takayanagi & Kurosaki 1999; Takayanagi et al. 1999). It should be noted that ab-initio calculations with DFT and CCSD methods lead to the absence of a barrier for the insertion (Balucani et al. 2009). The N(²D) + CH₄ reaction has been studied experimentally in detail by Takayanagi et al. (1999) leading to k(223−292 K) = 7.13  ×  10^-11  ×  exp( − 755/T) cm³ molecule^-1 s^-1. $\begin{array}{ccc}\begin{array}{c}\\ & & \\ & & \end{array}& & \end{array}$Umemoto et al. (1998) found a ratio between H₂CNH + H/NH + CH₃ equal to 0.8/0.3. Balucani et al. (2009) studied this reaction in a crossed beam experiment at high collision energy (above 22 kJ mol^-1 which corresponds to T  >  2600 K) suggesting an increasing NH branching ratio with the temperature and also H₂CNH and CH₃N production. As their results are ambiguous and correspond to high collisional energies, we preferred to use the Herron (1999) average value for the global rate constant (there is a typographical error in the Herron paper: it should be A = 4.8  ×  10^-11 instead of 4.8  ×  10^-12 in Table 3) associated with Umemoto et al. (1998) branching ratios:

A.4. H + H₂CN

The H + H₂CN reaction is a key reaction for both HCN and HNC production. The rate constant for this reaction has been measured to be greater than 7  ×  10^-11 cm³ molecule^-1 s^-1 (Nesbitt et al. 1990) and the authors also determined HD/(HCN+HNC) branching ratios for the H + D₂CN reaction, leading to HD/(HCN+HNC) = 5  ±  3. This reaction has two exothermic bimolecular exit channels: $\begin{array}{ccc}\begin{array}{c}\\ & & \\ & & \end{array}& & \end{array}$Previous theoretical calculations (Nguyen et al. 1996; Larson et al. 2006) show no barrier in the entrance valley for H₂CNH formation. We performed new ab-initio calculations at the DFT/M06-2X/vtz and MRCI+Q/vqz levels (Loison et al., in prep.) showing unambiguously that H₂CNH formation occurs in addition to the direct H atom abstraction without barriers in the entrance valley. The only possible evolution of the H₂CNH molecule is back-dissociation or H₂ + HNC formation, the H₂ + HCN exit channel involving a high energy exit barrier located well above the reactant energy and therefore forbidden at 150−200 K. H₂ + HNC production also involves a high exit barrier calculated close to the energy level of the reactants. The position of the transition state TS varies with the calculation level and has been found to range between +15 kJ mol^-1 and  − 9 kJ mol^-1 by comparison with the reactants energy (Zhou & Schlegel 2009; Nguyen et al. 1996). As this energy is critical to estimate HNC production, we calculate it at the MRCI+Q/vqz + ZPE(DFT/M06-2X/vtz) level leading to a TS energy equal to  −9 kJ mol^-1 below the reactants energy. To estimate the amount of HNC formed, we performed RRKM calculations (Loison et al., in prep.) with a TS energy given by the MRCI+Q/vqz calculations and the geometry and vibrational frequencies obtained at the DFT/M06-2X/vtz level. The main exit channel of the H₂CNH adduct is found to be back-dissociation with a yield greater than 90% except at high pressure when the three body stabilisation plays a role. This corresponds to pressures above those relevant for Titan’s upper atmosphere where H₂CN is supposed to be present at high altitude. For the H + D₂CN results of Nesbitt et al. (1990) we attribute the m/e = 3 signal, due to the HD molecule, to direct D atom abstraction, and the m/e = 27 signal, due to HNC (and also HCN coming from HNC isomerization), to the H + H₂CN  →  H₂CNH  →  H₂ + HNC pathway. Our RRKM calculations are then in reasonable agreement with experimental results considering the various uncertainties, particularly regarding the TS energy value. However additional uncertainties arise from the difference between the H₂CNH  →  H₂ + HNC and D₂CNH  →  D₂ + HNC dissociation rate constants, and also from the fact that HCN and HNC will be produced partly above the HCN  ⇌  HNC isomerization barrier. Considering a roughly statistical energy distribution in HNC and HCN, we estimate 20% of total HNC production. For the global rate constant, we chose the experimental one even if it is only a minimum value.

A.5. N(⁴S) + H₂CN

There is one experimental determination of the rate constant between 200 K and 363 K (Nesbitt et al. 1990) and one indirect branching ratio determination (Marston et al. 1989b). The possible exit channels are: $\begin{array}{ccc}\begin{array}{c}\\ & & \\ & & \\ & & \end{array}& & \end{array}$Cimas & Largo (2006) performed ab-initio calculations and found a small barrier for direct H atom abstraction, located at 11 kJ mol^-1 at the CCSD/p-vtz level but at  − 13 kJ mol^-1 at the G2 level. This negative barrier is due to the ZPE variation. We performed calculations at the UHF-M06-2X/VTZ level and also found a barrier for direct abstraction (+6.4 kJ mol^-1 without ZPE and  − 3.6 kJ mol^-1 including ZPE). No barrier was found for H₂CNN adduct formation quickly leading to CH₂ + N₂ formation. The good agreement between the calculated reaction enthalpy and the value derived from thermochemical data at 298 K (Baulch et al. 2005) gives us confidence in the calculations and in our conclusion that CH₂ + N₂ is an open exit channel. There is a large uncertainty about NH + HCN production however which is likely to be a minor but non negligible exit channel between 150 and 200 K.

Marston et al. (1989b) found that HCN (and/or HNC) formation is the main product in the N(⁴S) + CH₃ system, HCN being attributed to the result of the reaction sequence N(⁴S) + CH₃  →  H + H₂CN followed by N(⁴S) + H₂CN  →  NH + HCN. However in their experiment, N(⁴S) was in excess and the HCN could have been the result of the reaction sequence N(⁴S) + CH₃  →  H + H₂CN followed by N(⁴S) + H₂CN  →  CH₂ + N₂ and N(⁴S) + CH₂  →  H + HCN/HNC. The experimental rate constant (k(T) = 1.0  ×  10^-10  ×  exp( − 200/T) cm³ molecule^-1 s^-1 in the 200 − 363 K range) has large uncertainties (Nesbitt et al. 1990) and cannot be extrapolated to low temperature as it may be the result of the two reaction channels including the direct abstraction with a barrier. We point out that a simple capture rate model, with the dispersion term only, leads to a rate constant equal to k(N(⁴S)+H₂CN) = 4.5  ×  10^-10  ×  (T/300)^0.17  ×  (3/8) = 1.7  ×  10^-10  ×  (T/300)^0.17 cm³ molecule^-1 s^-1, 4 times higher than the experimental rate constant equal to 4.4  ×  10^-11 cm³ molecule^-1 s^-1 at 300 K. Considering the uncertainties in the experimental results we propose:

A.6. HCN + CH

The total rate constant for this reaction has been measured between 296 and 674 K equal to (5.0 ± 0.4)  ×  10^-11  ×  exp((500 ± 30)/T) cm³ molecule^-1 s^-1 (Zabarnick et al. 1991). This expression cannot be extrapolated to low temperature as it yields unrealistically large values. The measured k(296 K) is equal to 2.7  ×  10^-10 cm³ molecule^-1 s^-1. A simple capture treatment leads to a high rate constant (close to 6.0  ×  10^-10 cm³ molecule^-1 s^-1 at 298 K) due to the strong dipole-dipole interaction. We propose to scale the capture rate constant to the value at room temperature and to conserve the temperature dependence appropriate for the dipole-dipole interaction, leading to k(T) = 2.7  ×  10^-10  ×  (T/300)^-0.17 cm³ molecule^-1 s^-1. Theoretical calculations (Du & Zhang 2006) show no barrier for H₂CCN formation through CH + HCN  →  HCC(H)N  →  H₂CCN. H₂CCN can lead subsequently to $\begin{array}{ccc}\begin{array}{c}\\ & & \\ & & \end{array}& & \end{array}$Osamura & Petrie (2004) have determined that H + HCCN production has no barrier in the exit channel and that H₂ + CCN production exhibits a tight TS located only 6 kJ mol^-1 above this exit channel. These two exit channels should therefore have similar branching ratios to a first approximation, even if RRKM calculations are necessary to estimate more precisely their value.

A.7. N(²D) + HCN

The first electronically excited state of atomic nitrogen, N(²D) is known to react quickly with most radicals and molecules (Herron 1999). We performed ab-initio calculations at various levels for N(²D) attack on HCN in C_s and C_2v symmetry. The most attractive pathways have been found at MRCI+Q/vtz and R-CCSD(T)/vtz levels for C_2v geometry for which the ²B₁ and ²B₂ states show no barrier for HCNN adduct formation. The similar results for both methods coupled to the very good agreement between experiments (Ralchenko et al. 2011) and our calculations at the MRCI+Q/vtz level for N(⁴S)/N(²D)/N(²P) relative energies lead us to predict no barrier for this reaction. Further evolution of the adduct can lead to: $\begin{array}{ccc}\begin{array}{c}\\ & & \\ & & \\ & & \end{array}& & \end{array}$To estimate the branching ratio we can use the various calculations performed for the CH + N₂ system (Moskaleva et al. 2000; Berman et al. 2007) which strongly suggest that CH + N₂ is the most favored exit channel. As there is no barrier for only two surfaces, there is an electronic factor equal to 2/5. The simple capture rate constant model based on dispersion only leads to a rate constant close to (2/5) × 4.0 × 10^-10 × (T/300)^-0.17 = 1.6 × 10^-10 × (T/300)^-0.17 cm³ molecule^-1 s^-1. This rate constant is higher than those for similar reactions, N(²D) + singlet state molecules such as N(²D) + H₂O and N(²D) + NH₃ reactions which are also likely to occur without a barrier considering their high rate constant values at 300 K, 4.0  ×  10^-11 and 5.0  ×  10^-11 cm³ molecule^-1 s^-1 respectively (Takayanagi et al. 1998; Herron 1999). To evaluate more carefully the N(²D) + HCN rate constant, we performed MRCI+Q/vtz and R-CCSD(T)/vtz calculations for the N(²D) + NH₃ reaction also showing no barrier for ²B₁ and ²B₂ surfaces. The difference between the experimental value and the higher capture rate value may be due to several reasons such as a bottleneck in the entrance valley (the interaction potential is attractive mainly around linear attack) and/or eventually to a possible relaxation (N(²D)  →  N(⁴S)). The capture rate constant is a maximum value and therefore, the actual rate constant is likely to be smaller, around (4−6)  ×  10^-11 cm³ molecule^-1 s^-1 by comparison with N(²D) + H₂O and N(²D) + NH₃ reactions.

All Tables

All Figures

	Fig. 1 Eddy diffusion profiles from various photochemical models (Yung et al. 1984; Strobel et al. 1992; Toublanc et al. 1995; Lara et al. 1996; Wilson & Atreya 2004; Vinatier et al. 2007; Lavvas et al. 2008; Hörst et al. 2008; Krasnopolsky 2009).
In the text

	Fig. 2 Mixing ratio of HNC in the upper atmosphere of Titan. Nominal model (solid line) and the different acceptable profiles derived from recent Herschel observations (Moreno et al. 2011) (grey dotted lines). Local sensitivity analysis for the isomerization reaction H + HNC  →  HCN + H is also illustrated by changing its rate by a factor of 10 or 100.
In the text

	Fig. 3 Mixing ratio of HCN in the upper atmosphere of Titan. Model (solid line) and main observations in the upper atmosphere (see text for references).
In the text

	Fig. 4 Mixing ratios of HNC and HCN and HNC/HCN ratio as a function of altitude for two eddy diffusion coefficients.
In the text

	Fig. 5 Reaction rate profiles for most important reactions leading to the production and loss of HNC (see Table 2 for reactions).
In the text

	Fig. 6 Reaction rate profiles for most important reactions leading to the production and loss of HCN (see Table 3 for reactions).
In the text

	Fig. 7 Histograms of reaction rates obtained after 1000 runs for reaction N(4S) + C2H  →  C2N + H. The reaction rate estimated for this reaction is kcapture = 6.0  ×  10-11 and its uncertainty factor is Fcapture = 16.
In the text

	Fig. 8 Abundance profiles of HNC and HCN obtained after 1000 runs. Black solid line: initial profile. Black dotted line: median profile obtained from the uncertainty propagation study. Black dashed-dotted lines: 5th and 15th 20-quantiles of the distribution. Black long-dashed lines: 1st and 19th 20-quantiles of the distribution.
In the text

	Fig. 9 Abundance profiles of HNC/HCN ratios obtained after 1000 runs. Black solid line: initial profile. Black dotted line: median profile obtained from the uncertainty propagation study. Black dashed-dotted lines: 5th and 15th 20-quantiles of the distribution. Black long-dashed lines: 1st and 19th 20-quantiles of the distribution.
In the text