Ab initio study of nitrogen and position-specific oxygen kinetic isotope effects in the NO + O₃ reaction

Nitrogen oxides (NO_x = NO + NO₂) are important trace gases that influence the concentrations of atmospheric oxidants that drive tropospheric and stratospheric chemistry.^1–4 During the daytime, NO_x exists in a closed photochemical cycle between NO–O₂–O₃–NO₂ in the atmosphere, known as the Leighton cycle.^1,5 This cycle is initiated when NO₂ is photolyzed by UV-visible light in the blue region of the spectrum (< 400 nm) yielding O(³P) the ground state of the oxygen atom. This liberated oxygen atom can combine with O₂ to form O₃, which then oxidizes NO back to NO₂,^1,5

The analysis of the oxygen and nitrogen stable isotopes of NO_x and its oxidation product, atmospheric nitrate, may help in our understanding of this photochemical cycling^6–10 and sources of NO_x.^11–13 Variations in oxygen and nitrogen isotope compositions are reported using $δ(‰$$)$ notation where $δxO(‰$$)$ and $δ15N(‰$$)$ = (R_sample/R_ref $−$ 1) × 1000, where R_sample and R_ref denote the ^xO/¹⁶O (x = 17 or 18) or ¹⁵N/¹⁴N in the sample or reference, respectively. The oxygen isotopic reference is Vienna Standard Mean Ocean Water (VSMOW) and the nitrogen isotopic reference is atmospheric air. However, $δ18O$ and $δ15N$ of NO_x may also be influenced by isotopic fractionation processes associated with the Leighton cycle ((R1)–(R4)),^14,15 which may be propagated into atmospheric nitrate,¹⁶ yet few of these fractionation factors have been determined.

The photochemical cycling of NO_x is rapid, and prior experimental investigations of this cycling have suggested that isotopic equilibrium is achieved between O₃ and NO_x, erasing any original O isotopic NO_x signatures.¹⁰ Several studies have shown that atmospheric O₃ has an elevated $δ18O$ and a strong mass-independent component that is quantified by $Δ17O$ notation,^17–22 

In Eq. (1) the $λ$ is mass-dependent coefficient, which may be approximated as 0.52.²³ During the photochemical cycling of NO_x, these elevated $δ18O$ and $Δ17O$ signatures of O₃ are transferred to NO_x as a result of (R4).^8,10,24 While $Δ17O$ of the transferred O atom from O₃ to NO₂ should be minimally impacted, $δ18O$ may be significantly altered as a result of the mass-dependent fractionation associated with the kinetic isotope effect (KIE) of (R4). Thus, the kinetic isotope effect associated with (R4) may play an important role in the $δ18O$ of NO_x, which may be propagated into atmospheric nitrate, yet this fractionation factor is relatively unknown. Post-deposition isotope effects of nitrate such as photolysis of nitrate in a particle may alter the isotopic composition of nitrate,^25,26 but the focus on this study is to determine the isotope effect of (R4) which has implications for the isotopic composition of NO₂, which serves as precursor to atmospheric nitrate.

The kinetic isotope effect associated with (R4) may also play an important role in the $δ15N$ of NO₂. Previously, $δ15N$ of NO_x and atmospheric nitrate has been suggested to provide information about NO_x sources.^25–28 Numerous studies have quantified $δ15N$ from various NO_x sources, and these results indicate that soil emission (denitrification), transportation related sectors, and coal-fired power plants have relatively distinctive $δ15N$ values.^{11–13,25,29–35} These works have motivated several $δ15N$ studies of atmospheric nitrate as a way to partition NO_x sources to evaluate local/regional changes in NO_x source budgets.^27,28,36 However, the isotopic fractionation processes associated with the photochemical cycling of NO_x such as (R4), NO₂ photolysis, and NO_x isotopic exchange may alter the N isotopic composition of NO and NO₂ relative to total NO_x;¹⁵ however, except for NO_x isotope exchange, these fractionation processes’ impact on $δ15N$ is relatively unknown.¹⁵ If these fractionation processes are significant, daytime $δ15N$–$NO2$ may not equal $δ15N$–$NOx$⁠,^15,16 which has important implications for atmospheric nitrate formed during the daytime, because it is primarily formed through the reaction between NO₂ and photochemically produced $∙$OH.³⁷ Therefore, it is important to understand the kinetic isotope effect associated with (R4) and its implications for $δ15N$ of daytime produced atmospheric nitrate.

Bigeleisen demonstrated kinetic isotope effects (KIEs) could be approximated for reactions such as (R4) using transition state theory if the vibrational frequencies of the reactants and transition state are known.³⁸ Unfortunately, transition state frequencies for many isotopologues in NO_x related reactions, such as (R4), are unknown. Previously, ab initio methods have been used to calculate the reaction mechanism, thermochemistry, and vibrational frequencies of (R4),³⁹ but only for the main isotopologues. The present study builds on Ref. 39 by employing ab initio methods to calculate nitrogen and position-specific oxygen KIEs associated with (R4). This will allow for an understanding of the impact (R4) has on the $δ18O$ value of the transferred O atom from O₃ to NO and on $δ15N$ as NO is oxidized to NO₂.

The assumed reaction mechanism of (R4) is shown in Fig. 1, which is based on the results from Ref. 39. In the previous study, geometries and vibrational frequencies for all stationary points along the potential energy surface (PES) have been calculated by UHF, UMP2, and UMP4 methods with the 6-31G(d), 6-311G(d), and 6-311G(d,f) basis sets.³⁹ Unfortunately the UMP2 and UMP4 levels predicted anomalous vibrational frequencies for the radicals NO and NO₂ due to spin contamination and for O₃ due to its multirefernce character.³⁹ These inaccurate vibrational frequencies will have a significant impact on the accuracy of the calculated KIEs. Therefore, we have recalculated the geometries for a portion of the potential energy surface (PES) that included the reactants, products, and rate-determining transition state, which has been previously determined to be the NO radical approach to O₃ (TS1, Fig. 1),³⁹ using the high-level quantum mechanical method CCSD(T) with the 6-31G(d) and 6-311G(d) basis sets. The reactants and products geometry optimizations were carried out using default convergence criteria, while TS1 was optimized from calculated force constants at the previously reported UMP2(full)/6-31G(d) geometry.³⁹ Vibrational frequency analysis was performed to confirm that the obtained structures are stationary points along the PES that correspond to either a local minimum (3n-6 or 3n-5 real normal modes of vibration) or a transition state (exactly one imaginary frequency). Systematic model errors in the calculated harmonic frequencies have been accounted for by applying a constant scale factor of 0.9899 and 0.9542 for CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d), respectively (Fig. S1 of the supplementary material). All calculations were performed with the Gaussian09 program package revision D.01⁴⁰ on the Purdue Radon cluster.⁴¹ 

Nitrogen and position-specific oxygen KIEs were calculated in accordance with the Born-Oppenheimer, rigid-rotor, and harmonic approximations, using the complete Bigeleisen equations³⁸ as implemented in the ISOEFF program⁴² at 220, 250, 273, 298, and 320 K. Both the Bigeleisen equations and the ISOEFF program define the calculated KIE as the ratio of reaction coefficients of the light to heavy isotope. However, since we are interested in the relative rate of the heavy to light isotope, we report all calculated KIEs as the ratio of the heavy to the light isotope,

where $α$ is the calculated KIE, which is also referred to as the isotopic fractionation factor, and H and L refer to the heavy and light isotopes of a particular isotopologue pair, respectively. Corrections for tunneling and anharmonicity were neglected as these isotope effects tend to be small⁴³ and tend to cancel each other out for heavy atoms as anharmonicity tends to lower the magnitude of the calculated KIE while tunneling tends to increase the magnitude of the calculated KIE.⁴⁴ Additionally, neglect of tunnel effects may be validated due to the involved heavy atoms⁴⁵ and the observed low frequency modes along the reaction coordinate. The activation entropy $(ΔS‡)$ for this reaction is negative (−28.7 and −30.6 Cal*Mol⁻¹*K⁻¹ for CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d)), so that the reaction is characterized by a “tight transitions state,”⁴⁶ which should satisfy the rigid-rotor approximation used within the Bigeleisen equations. Overall, we expect the assumed approximations in to have a minimal impact on the calculated KIEs while pointing out that the calculated KIEs are an approximation.

Table I compares the geometries and harmonic frequencies for the most abundant isotopologues of NO, O₃, NO₂, and O₂ calculated by CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d) with experimental data.^47–54 Overall, the calculated geometries for the reactants and products are in excellent agreement with experimental data, as calculated bond lengths and bond angles are within 0.023 Å and 0.4° of experimental data^49,52,54 (Table I). Additionally, the calculated harmonic frequencies (scaled) are within a maximum difference of 92.2 cm⁻¹ and an average error of 29.9 cm⁻¹ and 43.9 cm⁻¹ relative to experimental data^47,48,50,51 at CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d), respectively. Cartesian coordinates of the optimized geometries and the scaled harmonic frequencies calculated for the major isotopologue of NO, O₃, NO₂, and O₂ from CCSD(T)-631G(d) and CCSD(T)/6-311G(d) computed force constants are available in the supplementary material (Tables S1 and S2).

An important point in accurately calculating KIEs is accounting for relative changes in vibrational energies due to substitution of a heavier isotope. Table II compares the differences in calculated vibrational zero point energies (ZPE = $12$h$∑νi$⁠) for various reactant and product ¹⁵N and ¹⁸O isotopologues relative to the major isotopologues for CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d) with available experimental data.^47,48,53,55 The difference in vibrational energies for each minor isotopologues relative to the most abundant is shown in parentheses in Table II. Our calculated relative difference in ZPE due to substitution of a heavier isotope is within 1.2 cm^$−$1 of experimental data and an average relative difference of 0.34 and 0.61 cm^$−$1 for CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d), respectively (Table II). This comparison shows the ability of our chosen level of theories to accurately reproduce changes in ZPE due to isotopic substitution. Therefore, CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d) should be sufficient methods to calculate the relative changes in thermodynamic properties of the N and O isotopologues necessary to determine the KIEs of (R4).

The calculated geometries of TS1 and atom labeling scheme is displayed in Fig. 2. As expected for an exothermic reaction, TS1 shows very small perturbations of geometrical parameters relative to the reactants.⁵⁶ Except for O₄, the atoms of TS1 are nearly planar as evident from the dihedral angle <(O₁–N₁–O₂–O₃) of −169.6° and −171.3° at CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d), respectively (Fig. 2), which is a similar finding to geometry of TS1 calculated at the UMP2/6-31(d) level of theory.³⁹ The imaginary frequency corresponding to the reaction coordinate is calculated to be 337.5i and 335.1i for CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d), respectively, which is in close agreement with the previous calculated value of 352.9i at UMP2/6-31G(d).³⁹ Cartesian coordinates of the optimized geometries and the scaled harmonic frequencies calculated for the major isotopologues of TS1 from CCSD(T)-631G(d) and CCSD(T)/6-311G(d) computed force constants are available in the supplementary material (Tables S3 and S4).

Calculated nitrogen and position-specific oxygen $α$ values at 298 K are reported in Table III as isotopic enrichment factors in units of per mil (%o),

Both the CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d) methods calculate similar $ε(‰)$ values that differed by no more than 1.6‰(Table III). The largest KIE is observed to occur for the substitution of ¹⁸O along the reaction coordinate (the ON + ¹⁸OOO $→$ ON¹⁸O + OO) that is calculated to be −43.9 and −44.7‰using CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d) calculated frequencies, respectively (Table III). The magnitude of the KIE is observed to decrease as the isotopic substitution position is further away for the reaction coordinate (i.e., secondary KIEs) as expected (Table III). Isotopic enrichment factors have also been calculated at 220, 250, 273, and 320 K and are provided in the supplementary material (Tables S5 and S6). It is important to note that we cannot claim perfection in these calculated KIEs due to numerous uncertainties associated with transition state theory (e.g., unknown tunneling effects, rigid-rotor and harmonic oscillator approximations) and the computational limitations (e.g., incomplete basis sets, uncertainties in the predicted transition state geometry). Despite these uncertainties, the excellent agreement between predicted geometries, vibrational frequencies, and shifts in ZPEs due to isotopic substitution with experimental data (Tables I and II) lend confidence in the calculated KIEs. Thus, these values should serve as a decent approximation that will be useful for future isotopic modeling of NO_x chemistry.

Oxygen mass-dependent relationships (ln(¹⁷α)/ln(¹⁸α))²³ for position-specific oxygen isotopic substitution are displayed in Table IV. Generally, (ln(¹⁷α)/ln(¹⁸α)) is found to range between 0.524 and 0.527, which is close to the generally accepted value of 0.520²³. However, (ln(¹⁷α)/ln(¹⁸α)) is found to have an anonymously large deviation from 0.520 for the ^xON + OOO $→$ ^xONO + O₂ reaction that is calculated to be 0.7072 and 0.4822 for CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d), respectively (Table V), which occurs due to both ¹⁸α and ¹⁷α being close to unity.⁵⁷ Despite this deviation, this reaction will have a minimal impact on $Δ17O$ (<0.05‰) as shown in Table IV since both ¹⁸$α$ and $α17$ are close to unity, where $Δ17O$ is calculated as

In Eq. (4), $λ$ is assumed to be 0.52.

Previously, Ref. 58 experimentally determined the kinetic isotopic fractionation associated with the NO + O₃ $→NO2+$O₂ reaction by reacting NO with excess O₃ and measuring the O isotopic composition of the unreacted O₃ and reported an overall $ε(‰$$)$ of −30.5‰ for ¹⁸O. This experimentally determined $ε(‰$$)$ value does not correspond to a singular KIE, rather it is roughly a statistical average for all the KIEs associated with the various ¹⁸O isotopomers of O₃ and their reactions with NO. In order to quantitatively evaluate our calculated KIEs with respect to the results obtained by Ref. 58, we modeled the kinetics for reactions of the various O₃ isotopomers with NO using a subset of a previously published NO_x cycle chemical kinetics model,¹⁰ utilizing Kintecus, a chemical kinetics compiler.⁵⁹ 

The NO_x cycling kinetics model previously published by Ref. 10 contains numerous facets of the NO_x cycle including NO₂ dissociation, oxygen isotope exchange, ozone formation, ozone dissociation, NO oxidation by O₃, NO oxidation by O-atom, NO₂ reaction with O-atom, NO₂ exchange with O-atom, NO exchange with O-atom, NO₂ exchange with NO, and NO oxidation by O₂. Here, we are only interested in the KIE associated with the reaction of the various isotopomers of O₃ with NO, so only the reactions pertaining to this reaction were used in the present study, which are displayed in Table V. In our model, we have only considered ¹⁸O isotopic substitution due to this reaction being a mass-dependent fractionation process. Since we are interested in the relative change of the ¹⁸O isotopic composition of O₃, the starting composition is unimportant, so we specified a starting $δ18O(O3)$ of 0‰ that was assumed to be uniformly distributed within O₃. While there is experiment evidence that $δ18O$ is not uniformly distributed within O₃,⁶⁰ this should have a minor impact on the relative change of bulk isotopic composition of $δ18O$ in our model. The reaction of NO with the main isotopologue of O₃ (¹⁶O¹⁶O¹⁶O) was set to 1.73 × 10⁻¹⁴ cm³*molecules⁻¹*s⁻¹,⁶¹ and the reactions of NO with the ¹⁸O substituted O₃ isotopologues/isotopomers were scaled by their position-specific oxygen $α$ value. Additionally, the reactions involving the asymmetric isotopic substituted O₃ isotopologue (¹⁸OOO) were scaled by 0.5 to account for the reaction channel symmetry (Table V). Our model was initiated with a NO:O₃ ratio of 0.95:1 so that O₃ was slightly in excess as in the experiment conducted by Ref. 58.

From the model output, a Rayleigh-type distillation curve of the $δ18O$ of unreacted O₃ was constructed in the form of

where $δf$ and $δ0$ are the initial and final $δ18O$ of O₃, respectively, and f is the fraction of unreacted O₃. Fig. 3 displays our calculated Rayleigh type distillation curve for O₃, where the slope of the line indicates an $ε(‰$$)$ value of −18.6‰ and −19.6‰ for CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d), respectively, which differs from the experimentally determined value of −30.5‰.⁵⁸ However, re-analysis of the experimental data from Ref. 59 indicates that there might be a high-leverage, influential x data point at ln(f) = −3.047. Fig. 4 compares two linear regression models of the experimental data from Ref. 59 which includes all of the data points (model A) and the omission of the data point at ln(f) = −3.047 (model B). While the R² values from the two models do not vary by much (R² = 0.966 and 0.969), the slopes drastically change from −0.0305 ± 0.003 to −0.0209 ± 0.003 for model A and model B, respectively (Fig. 4). The calculated $ε(‰)$ from model B is −20.9 ± 0.003‰, which is in excellent agreement with our calculated and modeled value at both CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d). This might suggest that the data point at ln(f) = −3.047 might be influenced by other isotopic fractionation processes such as the formation of higher order nitrogen oxide species such as NO₃ and N₂O₅ that might have a different $α$⁠. Also, it is also important to point out the limitations in our calculated $α$ values from the Bigeleisen equations that are obtained within the conventional transition state theory with harmonic normal modes and rigid rotor approximations.³⁸ Additionally, the experimental α may be influenced by the NO + O₃ reaction pathway proceeding through the O extraction from the apex O atom position of O₃ that recent experimental data shows may have a branching ratio of 8% ± 5%, which was neglected in the theoretial α.⁸ However, our calculated $α$ values tend to agree quite well with a significant portion of the experimental data.

Another important KIE to consider is the O transfer during the oxidation of NO by O₃ to NO₂, which will be referred to as $ε(O−trans)18$⁠,

where $δ18O(NO2(O−trans))$ and $δ18O(aO3)$ are the $δ18O$ values of the transferred O atom in NO₂ and asymmetric O₃, respectively. Previously, in Ref. 8, NO reacted with O₃ at a 1:1 ratio at room temperature and the $δ18O$ of the transferred O atom in NO₂ was measured. Using their experimental data, statistical models of previous studies of the intermolecular isotope distributions of O₃,⁶² and the enrichment of asymmetric and symmetric O₃ isotopologues in the stratosphere as a function of altitude,⁶³ $ε(O−trans)18$ was estimated to be −23.9 and −20.8‰.⁸ Under the experimental conditions (i.e., 1:1 ratio of NO:O₃), the products, NO₂ and O₂, reflect the partitioning of ¹⁸O based on the position-specific $α$ values for the ¹⁸O isotopomers of O₃, as there is a minimal isotopic fractionation impacting the residual O₃ since it nearly completely reacts.

Using our Kintecus model described in Sec. III C 1, we have estimated the $ε(O−trans)18$ using our O₃ isotopologue dependent KIEs and a NO:O₃ ratio of 1:1. The ¹⁸O isotopic composition of O₃ was assumed to be uniformly distributed with a starting $δ18O(O3)$ = $δ18O(aO3)$ = 0‰. From the output of our model, we estimate $δ18O(NO2(O−trans))$ and thus $ε(O−trans)18$ (Eq. (6)) to be −22.8‰ for both CCSD(T)/6-31G(d) and CCSD(T)/6-311G(d) calculated $α$ values, respectively. These values are in excellent agreement with the experimental determined values that are estimated to range between −23.9‰ and −20.8‰.⁸ 

Ab initio calculations have been carried out to investigate the nitrogen and oxygen KIEs associated with NO reaction with O₃. The calculated KIEs were generally close to unity except for primary KIEs, in which relatively large enrichment factors were calculated to be −44.7 and −6.7‰ for the ¹⁶O¹⁴N + ¹⁸O¹⁶O¹⁶O $→$ ¹⁶O¹⁴N¹⁸O + ¹⁶O¹⁶O and ¹⁶O¹⁵N + ¹⁶O¹⁶O¹⁶O $→$ ¹⁶O¹⁶N¹⁶O + ¹⁶O¹⁶O reactions, respectively, at 298 K (CCSD(T)/6-311G(d)). The reported KIEs may be limited due to the approximations used within the Bigeleisen equations and the computational limitations; however, excellent agreement between predicted geometries, vibrational frequencies, and shifts in ZPE due to isotopic substitution with experimental data suggests that our applied methods should provide a decent KIE approximation. Additionally, kinetic modeling of our calculated oxygen position-specific KIEs indicates the excellent agreement between our values and prior experimental measurements, which indicates that our applied methods are generally capturing the underlying physics correctly. Our calculations indicate that the O mass-dependent relationship $(ln(17α)/ln(18α))$ is generally near the expected value of 0.52. Cases in which $(ln(17α)/ln(18$α)) deviates from the expected value occur when the calculated KIEs are close to unity and thus have a minor impact on $Δ17O$⁠, as expected for a mass-dependent fractionation process. This indicates that NO reaction with O₃ may play a significant role in the $δ15N$ and $δ18O$ values of NO and NO₂, without altering $Δ17O$⁠. This has important implications for utilizing $δ15N$ and $δ18O$ as tools for NO_x source partitioning and for understanding NO_x photochemical cycling. The calculated KIEs will be useful for future work aimed at modeling NO_x isotope chemistry and will help guide future ambient NO_x isotopic measurements. Additionally, future work should aim to estimate the branching ratio of the NO reaction with O₃ at the apex O-atom position using ab initio methods.

See supplementary material for harmonic frequency scale factors, Cartesian coordinates for all calculated stationary points, calculated isotopologue harmonic frequencies, and calculated temperature dependent KIEs.

One of the authors (W.W.W.) was a National Science Foundation Graduate Research Fellow during the course of the study. We would like to thank the Purdue Climate Change Research Center (PCCRC) graduate fellowship program and the Purdue Bilsland Dissertation Fellowship for supporting this work. The authors would like to thank two anonymous reviewers for their helpful comments that significantly improved this manuscript.