# Mass-dependent and independent isotopic fractionation in ozone - theoretical modeling of experimental data

Peter Reinhardt, François Robert

► To cite this version:

Peter Reinhardt, François Robert. Mass-dependent and independent isotopic fractionation in ozone - theoretical modeling of experimental data. Chemical Physics Impact, 2026, 12, pp.101022. 10.1016/j.chphi.2026.101022. hal-05524544

HAL Id: hal-05524544

https://hal.sorbonne-universite.fr/hal-05524544v1

Submitted on 24 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Full length article

# Mass-dependent and independent isotopic fractionation in ozone — theoretical modeling of experimental data

Peter Reinhardt<sup>a, b,*</sup>, François Robert<sup>b</sup>

## ARTICLE INFO

### Keywords:

Ozone

Mass-independent isotope fractionation

Quantum mechanics

Model chemistry

Molecular dynamics

## ABSTRACT

An ad-hoc correction of classical-trajectory scattering simulations fully taking into account the quantum-mechanical principle of indiscernibility, produces a mass-independent isotope fractionation as observed in ozone.

Variations with respect to temperature are reproduced for two different potential surfaces. Coupling the mass-independent fractionation to mass-dependent corrections gives us the possibility to estimate the kinetic constant for ozone-stabilizing collisions with a third body, which in turn allows to access the pressure dependence of the isotope fractionation.

Obtained results are in qualitative agreement to published experimental data.

## Introduction

In chemical 2-body reactions, the partitioning of isotopes always obeys the Mass Dependent Fractionation theory (referred here to as MDF) ultimately related to the difference in mass of the different isotopes of the same element [1–3]. A prominent violation of this theoretical expectation was discovered by Thiemens and Heidenreich [4,5] for ozone, where the three oxygen isotopes were found to exhibit a Mass Independent Fractionation (MIF). The effect, which we will designate in the present paper as Chemical MIF (Chem-MIF), has been reproduced in laboratory experiments for oxygen isotopes [6–8] as well as for other chemical elements [9–11]. It has to be distinguished from other mass independent effects such as a Magnetic Isotope Effect (MIE) through coupling of electronic and nuclear spins [12,13] or the Nuclear Volume effect (NVE), [14] or other mass-independent effects [15–17].

In an isotopic exchange between two chemical species in equilibrium or not (noted Sample 1 and 2), the isotopic fractionation factor $\alpha_{A-B}$ is defined as the ratio of the isotopes A and B of Sample 1 with respect to Sample 2:

$$ \alpha_{A-B}^{1-2} = \frac{([A]/[B])_{\text{Sample 1}}}{([A]/[B])_{\text{Sample 2}}} \quad (1) $$

In the case of oxygen isotopes and for practical purposes, the Standard Mean Ocean Water (SMOW) is the reference ratio to define the $\delta$ units used in geochemistry:

$$ \delta(\textperthousand) = \left( \frac{R_{\text{Sample}}}{R_{\text{SMOW}}} - 1 \right) \times 1000 \quad (2) $$

with R being the isotopic ratios of the sample and of the SMOW, respectively. By definition, for the reference sample the molecular fractions $x_{16}$, $x_{17}$, and $x_{18}$ (or relative isotopic abundances of $^{16}$O, $^{17}$O and $^{16}$O) have to sum to 1.

Using the three isotopes the MDF effect can be introduced by

$$ \alpha_{17-16} = \alpha_{18-16}^{\beta} \quad (3) $$

For equilibrium exchange or for kinetic fractionation, small but measurable differences in the $\beta$ values are observed ($0.510 \le \beta \le 0.526$) which can be related [3] to the difference in isotope masses with respect to mass 16 as<fnref n="1" />

$$ \delta^{17}\text{O} = 0.52 \delta^{18}\text{O} \quad (4) $$

In contrast to this theoretical expectation, the Chem-MIF effect shows simply $\delta^{17}\text{O} \approx \delta^{18}\text{O}$.

<footnote><sup>a</sup> Laboratoire de Chimie Théorique, Sorbonne Université et CNRS, UMR7616, 4 place Jussieu, 75252 Paris CEDEX 05, France</footnote>

<footnote><sup>b</sup> Institut Origine et Evolution (O & E), Muséum National d'Histoire Naturelle et Sorbonne Université, IMPMC-UMR 7590 CNRS, 57 rue Cuvier, 75005 Paris, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: Peter.Reinhardt@sorbonne-universite.fr (P. Reinhardt).</footnote>

<footnote n="1">¹ $\delta^{17}\text{O} = (\alpha_{18-16}^{\beta} - 1) \times 1000 = 1000 \left( \left( 1 + \frac{\delta^{18}\text{O}}{1000} \right)^{\beta} - 1 \right) = 1000 (1 + \epsilon)^{\beta} - 1 \approx 1000((1 + \beta \epsilon) - 1) = \beta \delta^{18}\text{O}$</footnote>

<footnote n="2">² Multiplying $\alpha$ with a factor of 1.18 results immediately in huge $\delta$ values, shifting the origin in the 3-isotope plot to 180 ‰ on both axes. The data points are close to the $\delta^{17}\text{O} = \delta^{18}\text{O}$ MIF straight line.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. The Lindemann mechanism of ozone formation, with the two exit channels for an atom-molecule collision, ozone stabilized by a third body or re-dissociation of the intermediate complex.</figcaption></figure>

Since its discovery, the origin of the Chem-MIF has been the subject of extensive numerical computations. Hathorn and Marcus [18,19] calculated reaction rates for different isotope combinations, based on symmetry arguments [20]. In order to match experimental results, they suggest to multiply [21–24] the mass-dependent fractionation factor² with a constant $\eta = 1.18$:

$$ \alpha = \alpha_{\text{MD}} \times \eta \quad (5) $$

Detailed calculations by Schinke et al. [25] show that the data are better accounted for $\eta = 1.14$ rather than $\eta = 1.18$ because a mass dependent isotopic fractionation takes place during the dissociation of the complex.

Mauersberger et al. [26,27], and Janssen et al. [23] showed experimentally that the fractionation of asymmetric ozone isotopomers is correlated to the difference in zero-point energies of the two possible fragments A–B and B–C of a molecule A–B–C containing atoms of two isotopes (like 16-16-17 or 667).

Kinetic considerations, as advanced in the literature already, may bring us to a possible explanation of this enigma. Ozone is formed by a 3-body (Lindemann) reaction as shown in Fig. 1.

In steady-state conditions, the rate of formation of $O_3$ is

$$ \begin{gathered} \frac{d[O_3^*]}{dt} = 0 \\ k_f [O][O_2] = [O_3^*] k_D + [O_3^*] k_M [M] \\ \frac{d[O_3]}{dt} = [O_2][O] \frac{k_f k_M [M]}{k_D + k_M [M]} \end{gathered} \quad (6) $$

where $1/k_D$ is the lifetime of the complex. A lifetime sufficiently long to encounter a third body for stabilization of ozone constitutes channel 1 in the figure, otherwise channel 2 is followed, resulting in a re-decomposition of the complex $O_3^*$ (with the subscripts f, D and M for the formation, dissociation and stabilization rate constants). The low pressure approximation ($k_D \gg k_M [M]$) gives

$$ \frac{d[O_3]}{dt} = [O_2][O] \frac{k_f k_M [M]}{k_D} \quad (7) $$

In this paper, this approximation is valid for all the experimental data to which we compare our numerical simulations.

However, $\eta$ was introduced as an arbitrary parameter fitting the data and standing for the lifetime ratio of the various asymmetric isotopomers relative to $^{16}\text{O}^{16}\text{O}^{16}\text{O}^*$. Up to now, no convincing quantum mechanical calculations were able to reproduce the experimentally observed $\eta$ value [28]. As well, adding quantum ingredients of the potential surface (e.g. ZPE-weighted exit channels [29]) or selection rules [30,31] together with classical molecular-dynamics simulations [32] were not able to result in the mass-independent fractionations.

In contrast to these interpretations, we argued [33,34] that not only the quantum mechanics of the potential surface preparation influence the present system, but the fundamental principle of indistinguishability of identical quantum particles should apply to the scattering process $O + O_2 \rightarrow O_3^* \rightarrow O + O_2$. For classical particles we know exactly whether the outgoing atom is the incoming atom (i.e. scattering without isotope exchange; hereafter designated by NR for “non-reactive”) or whether the outgoing atom comes from the molecule (scattering

with isotope exchange; designated by R for “reactive”). However, in quantum mechanics, if the same isotopes are present as the incoming atom and as part of the molecule, the two processes R and NR have to be superposed for the scattering $A + AB \rightarrow A + AB$ (note that $A + AB \rightarrow B + AA$ is a classical R scattering) and for $A + AA \rightarrow A + AA$.

The ratio of the average lifetime of a “classical” (or distinguishable) scattering A + B-C, and the average lifetime of a “quantum” (or indistinguishable) scattering leads directly to the mass-independent $\eta$ factor, as we will show in the next section.

In the following we recall our suggested modification of “classical” lifetimes to “quantum” lifetimes, in order to arrive at the experimental data within classical molecular dynamics in a given potential surface. In the present study we will make use of two different potential surfaces, one of Bitterova and Schinke [35], and another of Tyuterev's et al. [36,37], including variations of temperature as measured by Morton et al. [38].

A numerical quantum counterpart of this interpretation is still pending.

## 1. Reaction model

In order to obtain a mass-independent reaction model, based on the quantum principle of indistinguishability only, we set in following all masses to 16 atomic mass units, and we label isotopes instead. Let us consider collisions between distinguishable isotopes:

$$ \begin{align} {}^{16}\text{O} + {}^{17}\text{O}^{18}\text{O} & \rightarrow ({}^{16}\text{O}^{17}\text{O}^{18}\text{O})^* \rightarrow {}^{16}\text{O} + {}^{17}\text{O}^{18}\text{O} \quad (\text{NR}) \nonumber \\ {}^{16}\text{O} + {}^{17}\text{O}^{18}\text{O} & \rightarrow ({}^{16}\text{O}^{17}\text{O}^{18}\text{O})^* \rightarrow {}^{17}\text{O} + {}^{16}\text{O}^{18}\text{O} \quad (\text{R}) \tag{8} \end{align} $$

where we do not consider the relative positions of the different isotopes in the complex (OOO)*. The reaction noted NR stands for a “non-reactive” collision i.e. a reaction through which isotopes are not exchanged while R stands for “reactive” i.e. an isotope exchange reaction between O and $O_2$. For now, we disregard the reaction giving atomic $^{18}\text{O}$. Although we cannot tell for a stabilized complex forming ozone whether the collision belongs to the R or the NR branch (no redissociation having taken place), we may yet calculate average lifetimes separately, and deduce, if a difference occurs, that the R or the NR full collision (if left to do so) will be preferentially stabilized, supposing that a formal relation may exist between initial conditions of a trajectory and the lifetime of the trajectory without the presence of an interfering third body for stabilization.

Initial conditions are selected in our simulations via a temperature-dependent distribution of relative velocities (atom – molecule) and rotational states of the $O_2$ molecule. The impact parameter obeys a geometrical consideration ($n(b) = 2\pi b$), with a completely random initial orientation of the $O_2$ molecule.

From the distribution $f(\tau)$ of obtained lifetimes $\tau$, we have in fact two separate distributions, $f_{\text{NR}}(\tau)$ and $f_{\text{R}}(\tau)$, respectively.<fnref n="3" /> Claiming

<footnote n="3">³ The function $f$ may be globally normalized ($\int_0^\infty f(\tau)d\tau = 1$) or normalized for each $\tau_{\text{min}}$ separately as $\int_{\tau_{\text{min}}}^\infty f(\tau)d\tau = 1$. Even if the definition of P uses integrals of different limits, i.e. the functions in numerator and denominator are differently normalized, quantities within same limits cancel out.</footnote>
that trajectories may result in stabilized ozone, if the lifetime is longer than a minimum lifetime $\tau_{min}$, we can form two separate averages

$$ 
\begin{aligned}
\bar{\tau}_R(\tau_{\min}) &= \int_{\tau_{\min}}^{\infty} \tau f_R(\tau) d\tau \bigg/ \int_{\tau_{\min}}^{\infty} f_R(\tau) d\tau \quad \text{and} \\
\bar{\tau}_{NR}(\tau_{\min}) &= \int_{\tau_{\min}}^{\infty} \tau f_{NR}(\tau) d\tau \bigg/ \int_{\tau_{\min}}^{\infty} f_{NR}(\tau) d\tau \quad .
\end{aligned}
\quad (9) $$

Combining these separate averages with appropriate weight functions, we arrive easily at the classical overall average $\bar{\tau}$

$$ 
\begin{aligned}
\bar{\tau}(\tau_{\min}) &= x(\tau_{\min}) \bar{\tau}_{NR}(\tau_{\min}) + (1 - x(\tau_{\min})) \bar{\tau}_R(\tau_{\min}) \\
&= \bar{\tau}_D(\tau_{\min})
\end{aligned}
\quad (10) $$

with

$$ x(\tau_{\min}) = \frac{X_{NR} P_{NR}(\tau_{\min})}{X_R P_R(\tau_{\min}) + X_{NR} P_{NR}(\tau_{\min})} \quad (11) $$

where we introduced two auxiliary functions, X and P defined as

$$ P_{NR}(\tau_{\min}) = \frac{\int_{\tau_{\min}}^{\infty} f_{NR}(\tau) d\tau}{\int_{0}^{\infty} f_{NR}(\tau) d\tau} \quad \text{and} \quad (12) $$

$$ X_{NR} = \frac{\int_{0}^{\infty} f_{NR}(\tau) d\tau}{\int_{0}^{\infty} (f_{R}(\tau) + f_{NR}(\tau)) d\tau} \quad (13) $$

and corresponding quantities $X_R$ and $P_R(\tau_{\min})$. Trivially, $X_R + X_{NR} = 1$.

The relation Eq. (10) is exact for calculating an average lifetime in any ensemble of trajectories, as we can write

$$ 
\begin{aligned}
\bar{\tau}(\tau_{\min}) &= \frac{1}{N} \sum_{i, \tau_i > \tau_{\min}} \tau_i = \frac{\int_{\tau_{\min}}^{\infty} \tau f(\tau) d\tau}{\int_{\tau_{\min}}^{\infty} f(\tau) d\tau} \\
&= \frac{\int_{\tau_{\min}}^{\infty} \tau (f_R(\tau) + f_{NR}(\tau)) d\tau}{\int_{\tau_{\min}}^{\infty} (f_R(\tau) + f_{NR}(\tau)) d\tau} \quad .
\end{aligned}
\quad (14) $$

Let us consider now the reaction between identical isotopes such as

$$ ^{16}\text{O} + ^{16}\text{O}^{16}\text{O} \rightarrow (^{16}\text{O}^{16}\text{O}^{16}\text{O}) * \rightarrow ^{16}\text{O} + ^{16}\text{O}^{16}\text{O} \quad (15) $$

For this kind of collision it is impossible to identify whether a R or a NR branch is followed. Quantum mechanics demands a superposition of both processes. We may introduce this by suggesting to replace $P_R(\tau_{\min})$ and $P_{NR}(\tau_{\min})$ by one common (or averaged) quantity $P_I(\tau_{\min})$.

Inserting this in our expressions for the average lifetime, Eq. (10), $P_I(\tau_{\min})$ drops out, and we arrive at

$$ \bar{\tau}_I(\tau_{\min}) = X_{NR} \bar{\tau}_{NR}(\tau_{\min}) + (1 - X_{NR}) \bar{\tau}_R(\tau_{\min}) \quad (16) $$

For an isotope exchange reaction passing through an ozone molecule as

$$ ^{17}\text{O} + ^{16}\text{O}^{16}\text{O} \leftrightarrow (^{16}\text{O}^{16}\text{O}^{17}\text{O}) \leftrightarrow ^{16}\text{O} + ^{17}\text{O}^{16}\text{O} $$

one may collect $O_3$ samples and analyze different masses, leading to another $\alpha$ notation than Eq. (1) by looking at a particular scrambled combination of isotopes XYZ:

$$ \alpha(\text{XYZ}) = \frac{\left( \frac{\left( [X^O Y^O Z^O O] / [^{16}\text{O}^{16}\text{O}^{16}\text{O}] \right)_{\text{sample}}}{\left( [X^O Y^O Z^O O] / [^{16}\text{O}^{16}\text{O}^{16}\text{O}] \right)_{\text{stat}}} \right)}{\left( \frac{\left( [X^O Y^O Z^O O] / [^{16}\text{O}^{16}\text{O}^{16}\text{O}] \right)_{\text{sample}}}{\left( [X^O Y^O Z^O O] / [^{16}\text{O}^{16}\text{O}^{16}\text{O}] \right)_{\text{stat}}} \right)} \quad (17) $$

Scrambled means to consider all possible isotopomers XYZ, YZX and ZXY together.

Assuming that oxygen isotopes in free atoms and $O_2$ molecules are statistically distributed, as $[^{16}\text{O}] = x_{16}[\text{O}]$, $[^{16}\text{O}^{17}\text{O}] = 2x_{17}x_{16}[\text{O}_2]$, $[^{16}\text{O}^{16}\text{O}] = x_{16}x_{16}[\text{O}_2]$ etc., with $x_{16} + x_{17} + x_{18} = 1$, we expect that the statistical ratio $^{17}\text{O}/^{16}\text{O}$ in ozone is the same as in the reservoir, i.e. $x_{17}/x_{16}$.

From Eq. (7) we assimilate rate constant ratios to abundance ratios as

$$ \left( \frac{[^{16}\text{O}^{17}\text{O}^{18}\text{O}]}{[^{16}\text{O}^{16}\text{O}^{16}\text{O}}} \right)_{\text{sample}} = \frac{d \left( \frac{[^{16}\text{O}^{17}\text{O}^{18}\text{O}]}{d\tau} \right)}{d \left( \frac{[^{16}\text{O}^{16}\text{O}^{16}\text{O}]}{d\tau} \right)} \quad (18) $$

<caption>Table 1</caption>

Possible situations of an isotope A scattering with an $O_2$ molecule. The type of scattering (I for indistinguishable and D for distinguishable), and the number of involved A atoms and B atoms are indicated in the last three columns.

<table><thead><tr><th>atom</th><th>molecule</th><th>[O]</th><th>[O<sub>2</sub>]</th><th>type</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>A</td><td>A-A</td><td>x<sub>A</sub></td><td>x<sub>A</sub><sup>2</sup></td><td>I</td><td>3</td><td>0</td></tr><tr><td>A</td><td>A-B</td><td>x<sub>A</sub></td><td>2x<sub>A</sub>x<sub>B</sub></td><td>D</td><td>2</td><td>1</td></tr><tr><td>A</td><td>A-C</td><td>x<sub>A</sub></td><td>2x<sub>A</sub>x<sub>C</sub></td><td>I</td><td>2</td><td>0</td></tr><tr><td>A</td><td>B-B</td><td>x<sub>A</sub></td><td>x<sub>B</sub><sup>2</sup></td><td>D</td><td>2</td><td>0</td></tr><tr><td>A</td><td>B-C</td><td>x<sub>A</sub></td><td>2x<sub>B</sub>x<sub>C</sub></td><td>D</td><td>1</td><td>1</td></tr><tr><td>A</td><td>C-C</td><td>x<sub>A</sub></td><td>x<sub>C</sub><sup>2</sup></td><td>D</td><td>1</td><td>0</td></tr></tbody></table>

and divide by the statistically expected abundances. The mass-independent part of the fractionation, $\eta$, becomes then just the ratio of the dissociation rates $k_D$ in Eq. (7) for which we insert the inverse lifetimes of the complexes. For the scrambled molecules 16-17-18, 17-18-16 and 18-16-17 from six possible collisions A + B-C we have then, together with the statistically expected abundances

$$ \alpha(\{16, 17, 18\}) = \eta(\tau_{\min}) = \frac{\bar{\tau}_{D}(\tau_{\min})}{\bar{\tau}_{I}(\tau_{\min})} \quad (19) $$

This result assumes the validity of the low pressure approximation as a limit to the more general expression (see Eq. (6)). The term $k_m[M]$ to be added in numerator and denominator in Eq. (19) attenuates the fractionation, leading for high pressures to a complete disappearance of an observable MIF [38,39].

For a scrambled mixture of 16-16-17 and 16-17-16 isotopomers we have to take into account a distinguishable and an indistinguishable component, leading to

$$ \alpha(\{16-16-17, 16-17-16\}) = \frac{1}{3} \frac{2 \bar{\tau}_D(\tau_{\min}) + \bar{\tau}_I(\tau_{\min})}{\bar{\tau}_I(\tau_{\min})} = \frac{2}{3} \eta(\tau_{\min}) + \frac{1}{3} \quad (20) $$

Exploiting Table 1 for accounting correctly atoms present and reaction types, we find (see B) for the fractionation of two isotopes in a completely scrambled ozone sample with all possible isotopomers

$$ \alpha_{17-16} = \frac{\bar{\tau}_I k_1 + 2\bar{\tau}_D k_2}{\bar{\tau}_I k_3 + 2\bar{\tau}_D k_4} \quad \text{where} \quad (21) $$

$$ k_1 = x_{16}^2 + 3x_{17}^2 + x_{18}^2 + 2x_{17} (x_{16} + x_{18}) \quad (22) $$

$$ k_2 = x_{16}^2 + 3x_{16}x_{18} + x_{18}^2 + 2x_{17}(x_{16} + x_{18}) \quad (23) $$

$$ k_3 = 3x_{16}^2 + x_{17}^2 + x_{18}^2 + 2x_{16}(x_{17} + x_{18}) \quad (24) $$

$$ k_4 = x_{17}^2 + 3x_{17}x_{18} + x_{18}^2 + 2x_{16}(x_{17} + x_{18}) \quad (25) $$

For natural abundances ($x_{16} \approx 1$, $x_{17} \approx 0$, $x_{18} \approx 0$) the expression reduces to

$$ \alpha_{17-16} = \frac{2 \bar{\tau}_D}{3 \bar{\tau}_I} + \frac{1}{3} = \frac{2}{3} \eta + \frac{1}{3} \quad (26) $$

Exactly the same holds for $\alpha_{18-16}$ as isotope masses are set identically to 16 amu, and thus $\delta^{18}\text{O} = \delta^{17}\text{O}$.

However, the value of $\eta$ is sensitive to several parameters: (i) the shape of the interaction potential describing the reaction O + $O_2$, (ii) the minimum lifetime $\tau_{\min}$ and (iii) the geometrical criterion to define the lifetime of the complex formed of the three oxygen atoms. These we will study in the next section.

## 2. Results

Classical trajectory simulations were carried out in two different potential surfaces, one given by Bitterova and Schinke [35] and one by Tyuterev et al. [36,37]. In order to be able to measure lifetimes, a criterion has to be chosen when to speak of a complex. We opted
<caption>**Table 2**</caption>

Example or data for a complex-defining radius of 4.6 bohr, the Bitterova-Schinke surface and a temperature of 300 K $X_{NR}$ and $X_{R}$ amount to 0.9045 and 0.095, respectively. In total we have 141321 trajectories, out of which 13493 are of R type and 127828 without particle exchange.

<table><thead><tr><th>τ<sub>min</sub></th><th>P<sub>NR</sub></th><th>P<sub>R</sub></th><th>τ&#773;<sub>NR</sub></th><th>τ&#773;<sub>R</sub></th><th>τ&#773;<sub>D</sub></th><th>τ&#773;<sub>I</sub></th><th>η</th></tr></thead><tbody><tr><td>0.00</td><td>1.0000</td><td>1.0000</td><td>2.9131</td><td>27.6849</td><td>5.2783</td><td>5.2783</td><td>1.0000</td></tr><tr><td>0.01</td><td>0.9947</td><td></td><td>2.9287</td><td></td><td>5.3038</td><td>5.2924</td><td>1.0022</td></tr><tr><td>0.02</td><td>0.9780</td><td></td><td>2.9783</td><td></td><td>5.3852</td><td>5.3373</td><td>1.0090</td></tr><tr><td>0.03</td><td>0.9500</td><td></td><td>3.0655</td><td></td><td>5.5275</td><td>5.4161</td><td>1.0206</td></tr><tr><td>0.04</td><td>0.9089</td><td></td><td>3.2024</td><td></td><td>5.7498</td><td>5.5399</td><td>1.0379</td></tr><tr><td>0.05</td><td>0.8539</td><td></td><td>3.4058</td><td></td><td>6.0769</td><td>5.7239</td><td>1.0617</td></tr><tr><td>0.06</td><td>0.7866</td><td></td><td>3.6924</td><td></td><td>6.5311</td><td>5.9832</td><td>1.0916</td></tr><tr><td>0.07</td><td>0.7115</td><td></td><td>4.0751</td><td></td><td>7.1252</td><td>6.3293</td><td>1.1257</td></tr><tr><td>0.08</td><td>0.6348</td><td></td><td>4.5590</td><td></td><td>7.8564</td><td>6.7670</td><td>1.1610</td></tr><tr><td>0.09</td><td>0.5604</td><td></td><td>5.1528</td><td></td><td>8.7244</td><td>7.3041</td><td>1.1944</td></tr><tr><td>0.10</td><td>0.4949</td><td></td><td>5.8216</td><td></td><td>9.6648</td><td>7.9090</td><td>1.2220</td></tr><tr><td>0.11</td><td>0.4381</td><td></td><td>6.5635</td><td></td><td>10.6647</td><td>8.5802</td><td>1.2429</td></tr><tr><td>0.12</td><td>0.3873</td><td></td><td>7.4096</td><td></td><td>11.7524</td><td>9.3455</td><td>1.2575</td></tr><tr><td>0.13</td><td>0.3444</td><td></td><td>8.3163</td><td></td><td>12.8600</td><td>10.1655</td><td>1.2651</td></tr><tr><td>0.14</td><td>0.3065</td><td></td><td>9.3285</td><td></td><td>14.0310</td><td>11.0811</td><td>1.2662</td></tr><tr><td>0.15</td><td>0.2753</td><td></td><td>10.3686</td><td></td><td>15.1679</td><td>12.0220</td><td>1.2617</td></tr><tr><td>0.16</td><td>0.2484</td><td></td><td>11.4738</td><td></td><td>16.3080</td><td>13.0216</td><td>1.2524</td></tr><tr><td>0.17</td><td>0.2245</td><td></td><td>12.6767</td><td></td><td>17.4760</td><td>14.1097</td><td>1.2386</td></tr><tr><td>0.18</td><td>0.2039</td><td></td><td>13.9442</td><td></td><td>18.6316</td><td>15.2561</td><td>1.2213</td></tr><tr><td>0.19</td><td>0.1867</td><td></td><td>15.2103</td><td></td><td>19.7160</td><td>16.4014</td><td>1.2021</td></tr><tr><td>0.20</td><td>0.1716</td><td></td><td>16.5330</td><td></td><td>20.7806</td><td>17.5978</td><td>1.1809</td></tr><tr><td>0.22</td><td>0.1516</td><td></td><td>18.6826</td><td></td><td>22.3777</td><td>19.5421</td><td>1.1451</td></tr><tr><td>0.24</td><td>0.1404</td><td></td><td>20.1627</td><td></td><td>23.3916</td><td>20.8809</td><td>1.1202</td></tr><tr><td>0.26</td><td>0.1330</td><td>1.0000</td><td>21.2637</td><td>27.6849</td><td>24.1050</td><td>21.8768</td><td>1.1018</td></tr><tr><td>0.28</td><td>0.1278</td><td>0.9999</td><td>22.1183</td><td>27.6889</td><td>24.6380</td><td>22.6502</td><td>1.0878</td></tr><tr><td>0.30</td><td>0.1237</td><td>0.9997</td><td>22.8446</td><td>27.6930</td><td>25.0767</td><td>23.3075</td><td>1.0759</td></tr><tr><td>0.32</td><td>0.1201</td><td>0.9991</td><td>23.5108</td><td>27.7093</td><td>25.4735</td><td>23.9117</td><td>1.0653</td></tr><tr><td>0.34</td><td>0.1174</td><td>0.9986</td><td>24.0404</td><td>27.7235</td><td>25.7825</td><td>24.3921</td><td>1.0570</td></tr><tr><td>0.36</td><td>0.1155</td><td>0.9972</td><td>24.4302</td><td>27.7621</td><td>26.0185</td><td>24.7483</td><td>1.0513</td></tr><tr><td>0.38</td><td>0.1136</td><td>0.9953</td><td>24.8411</td><td>27.8131</td><td>26.2690</td><td>25.1248</td><td>1.0455</td></tr><tr><td>0.40</td><td>0.1119</td><td>0.9931</td><td>25.2120</td><td>27.8745</td><td>26.4998</td><td>25.4662</td><td>1.0406</td></tr></tbody></table>

for the second-largest interatomic distance, $r_2$, creating thus a two-dimensional parameter space with $τ_{min}$. For a given temperature we run half a million of trajectories and collect data for establishing the lifetime distribution functions $f_R(τ)$ and $f_{NR}(τ)$. Only trajectories with a non-zero lifetime participate in the averaging process, excluding thus trajectories with too large impact parameters $b$ for generating an ozone molecule. For a temperature of 300 K Table 2 gives some details of the ingredients of $η$.

In Fig. 2 we trace our obtained fractionation $η$ against $τ_{min}$, for different limiting radii $r_{limit}$.

Tables 2 and 3 present the different ingredients for the calculation of $η$, for the two potential surfaces considered, and one limiting radius of $r_{limit} = 4.6$ bohr.

Several observations can be formulated:

* Limiting radii are quite large (between 4 and 5 bohr) for observing a notable fractionation, and significantly larger than the O–O bond in ozone, which is about half of this size.

* Minimum lifetimes of around 0.1 to 0.2 ps are shorter than minimal lifetimes for a collision with particle exchange (R collisions). This means that all R collisions are considered for the stabilization of ozone, and a significant part of rapid NR collisions which do not form a real three-atom molecule with two comparable (covalent) O–O bond lengths

* As all R collisions are taken into account, the average lifetime $\bar{τ}_R$ is nearly invariant in the region of minimum lifetimes of interest, and much longer than the “classical” ($\bar{τ}_D$) and the “quantum” ($\bar{τ}_I$) average lifetime (see Tables 2 and 3).

* As a consequence the average NR lifetime is much shorter than the obtained “classical” average lifetime and the “quantum” average lifetime.

<caption>**Table 3**</caption>

Example or data for a complex-defining radius of 4.6 bohr, the Tyuterev surface and a temperature of 300 K $X_{NR}$ and $X_{R}$ amount to 0.7964 and 0.2036, respectively. In total we have 51259 trajectories, out of which 10439 are of R type and 40820 without particle exchange.

<table><thead><tr><th>τ<sub>min</sub></th><th>P<sub>NR</sub></th><th>P<sub>R</sub></th><th>τ&#773;<sub>NR</sub></th><th>τ&#773;<sub>R</sub></th><th>τ&#773;<sub>D</sub></th><th>τ&#773;<sub>I</sub></th><th>η</th></tr></thead><tbody><tr><td>0.00</td><td>1.0000</td><td>1.0000</td><td>5.5838</td><td>23.5757</td><td>9.2479</td><td>9.2479</td><td>1.0000</td></tr><tr><td>0.01</td><td>0.9944</td><td></td><td>5.6153</td><td></td><td>9.2894</td><td>9.2730</td><td>1.0018</td></tr><tr><td>0.02</td><td>0.9783</td><td></td><td>5.7076</td><td></td><td>9.4105</td><td>9.3464</td><td>1.0069</td></tr><tr><td>0.03</td><td>0.9526</td><td></td><td>5.8608</td><td></td><td>9.6101</td><td>9.4685</td><td>1.0150</td></tr><tr><td>0.04</td><td>0.9190</td><td></td><td>6.0739</td><td></td><td>9.8841</td><td>9.6382</td><td>1.0255</td></tr><tr><td>0.05</td><td>0.8779</td><td></td><td>6.3559</td><td></td><td>10.2404</td><td>9.8627</td><td>1.0383</td></tr><tr><td>0.06</td><td>0.8333</td><td></td><td>6.6930</td><td></td><td>10.6574</td><td>10.1312</td><td>1.0519</td></tr><tr><td>0.07</td><td>0.7815</td><td></td><td>7.1322</td><td></td><td>11.1863</td><td>10.4809</td><td>1.0673</td></tr><tr><td>0.08</td><td>0.7282</td><td></td><td>7.6493</td><td></td><td>11.7888</td><td>10.8927</td><td>1.0823</td></tr><tr><td>0.09</td><td>0.6796</td><td></td><td>8.1894</td><td></td><td>12.3961</td><td>11.3229</td><td>1.0948</td></tr><tr><td>0.10</td><td>0.6365</td><td></td><td>8.7381</td><td></td><td>12.9908</td><td>11.7598</td><td>1.1047</td></tr><tr><td>0.11</td><td>0.5968</td><td></td><td>9.3129</td><td></td><td>13.5914</td><td>12.2175</td><td>1.1125</td></tr><tr><td>0.12</td><td>0.5602</td><td></td><td>9.9130</td><td></td><td>14.1951</td><td>12.6954</td><td>1.1181</td></tr><tr><td>0.13</td><td>0.5287</td><td></td><td>10.4967</td><td></td><td>14.7607</td><td>13.1603</td><td>1.1216</td></tr><tr><td>0.14</td><td>0.5022</td><td></td><td>11.0442</td><td></td><td>15.2727</td><td>13.5962</td><td>1.1233</td></tr><tr><td>0.15</td><td>0.4803</td><td></td><td>11.5406</td><td></td><td>15.7223</td><td>13.9916</td><td>1.1237</td></tr><tr><td>0.16</td><td>0.4610</td><td></td><td>12.0174</td><td></td><td>16.1416</td><td>14.3713</td><td>1.1232</td></tr><tr><td>0.17</td><td>0.4430</td><td></td><td>12.4978</td><td></td><td>16.5522</td><td>14.7539</td><td>1.1219</td></tr><tr><td>0.18</td><td>0.4274</td><td></td><td>12.9485</td><td></td><td>16.9269</td><td>15.1127</td><td>1.1200</td></tr><tr><td>0.19</td><td>0.4129</td><td></td><td>13.3960</td><td></td><td>17.2894</td><td>15.4691</td><td>1.1177</td></tr><tr><td>0.20</td><td>0.4006</td><td></td><td>13.8021</td><td></td><td>17.6103</td><td>15.7925</td><td>1.1151</td></tr><tr><td>0.22</td><td>0.3768</td><td>1.0000</td><td>14.6611</td><td>23.757</td><td>18.2654</td><td>16.4766</td><td>1.1086</td></tr><tr><td>0.24</td><td>0.3600</td><td>0.9999</td><td>15.3338</td><td>23.780</td><td>18.7577</td><td>17.0128</td><td>1.1026</td></tr><tr><td>0.26</td><td>0.3486</td><td>0.9995</td><td>15.8267</td><td>23.5869</td><td>19.1096</td><td>17.4071</td><td>1.0978</td></tr><tr><td>0.28</td><td>0.3397</td><td>0.9990</td><td>16.2340</td><td>23.5981</td><td>19.3950</td><td>17.7337</td><td>1.0937</td></tr><tr><td>0.30</td><td>0.3322</td><td>0.9984</td><td>16.5962</td><td>23.6137</td><td>19.6460</td><td>18.0253</td><td>1.0899</td></tr><tr><td>0.32</td><td>0.3256</td><td>0.9968</td><td>16.9258</td><td>23.6496</td><td>19.8785</td><td>18.2951</td><td>1.0865</td></tr><tr><td>0.34</td><td>0.3206</td><td>0.9938</td><td>17.1833</td><td>23.7215</td><td>20.0743</td><td>18.5148</td><td>1.0842</td></tr><tr><td>0.36</td><td>0.3159</td><td>0.9900</td><td>17.4352</td><td>23.8097</td><td>20.2713</td><td>18.7334</td><td>1.0821</td></tr><tr><td>0.38</td><td>0.3120</td><td>0.9856</td><td>17.6496</td><td>23.9145</td><td>20.4494</td><td>18.9255</td><td>1.0805</td></tr><tr><td>0.40</td><td>0.3078</td><td>0.9808</td><td>17.8804</td><td>24.0294</td><td>20.6412</td><td>19.1326</td><td>1.0788</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Comparison of results for the potential surfaces of Bitterova et al. (black lines), and Tyuterev et al. (red lines), for limiting radii of 4.2, 4.4, 4.6, and 4.8 bohr, respectively. $η$ is increasing for increasing $r_{limit}$. Data are for 300 K.</figcaption></figure>

Looking at the two potential surfaces, see Fig. 3, the complex-limiting radius of 4.6 bohr seems reasonable as at that distance the potentials present the entry region for the stabilizing valley, with (Bitterova–Schinke) or without (Tyuterev) a barrier.

Temperature dependence is added to our model via initial conditions. Fixing the limiting radius at 4.6 bohr and a minimum lifetime for stabilization at 0.14 ps (the maximum in the two tables), we see a temperature dependence of our model as displayed in Fig. 4. It should be noticed that, contrary to all known isotopic exchange reactions, the isotopic fractionation of ozone increases with the temperature.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Plot of the two potential surfaces for a given O–O distance of 2.409 bohr.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Temperature dependence of the calculated mass-independent fractionation for r<sub>limit</sub> = 4.6 bohr and a minimum lifetime of τ<sub>min</sub> = 0.14 ps.</figcaption></figure>

We may as well look into the pressure dependence, at a given temperature, by adding k<sub>M</sub> [M]. For a given pressure and temperature we can easily calculate [M] assuming for instance an ideal gas. However, we do not dispose of the constant k<sub>M</sub>, which we may extract from experimental data. To this end we have to look first at corrections for mass-dependent effects of our determined MIF.

### 3. Introducing the mass-dependent fractionation

In the previous section we derived the principal fractionation effect, which is mass-independent. However, experimental results of Morton et al. [38] show as well differences between <sup>17</sup>O and <sup>18</sup>O, in a smaller range. These experiments were conducted at temperatures between 150 and 350 K, for an O<sub>2</sub> pressure of 50 Torr. In this domain the low-pressure approximation applies [38]. To model these variations with temperature we call upon experimental rate coefficients relative to <sup>16</sup>O + <sup>16</sup>O<sup>16</sup>O measured at 300 K and 200 Torr [23] for the reaction

$$ \begin{aligned} & <sup>16</sup>O + <sup>16</sup>O<sup>18</sup>O \rightarrow <sup>16</sup>O<sup>16</sup>O<sup>18</sup>O \quad (i) \\ & <sup>18</sup>O + <sup>16</sup>O<sup>16</sup>O \rightarrow <sup>16</sup>O<sup>16</sup>O<sup>18</sup>O \quad (ii) \end{aligned} \quad (27) $$

Corresponding rate coefficient ratios are α<sub>6–68</sub> = 1.45 and α<sub>8–66</sub> = 0.92, respectively.

### 3.1. Eliminating the chem-MIF effect

These rates were corrected by Janssen et al. [23] for the expected contribution of the mass dependent isotopic exchange between O<sub>2</sub> and O<sub>3</sub>. A measured η<sub>ABB</sub> of 1.14 corresponds thus to an overall η of 1.21 (see Eq. (26)).

Following the procedure of Gao and Marcus, we calculate the mass dependent rate coefficients ratio with respect to 16 + 16–16, α<sub>MD</sub>, for the reactions (i) and (ii) as

$$ \begin{aligned} \alpha_{6-68,MD} &= 1.45/1.14 = 1.272 \\ \alpha_{8-66,MD} &= 0.92/1.14 = 0.807 \end{aligned} \quad (28) $$

The subscript MD indicates that the α values are the mass-dependent part of the fractionation. In the natural abundances of oxygen isotopes (<sup>16</sup>O ≈ 100%), the reactions involving <sup>18</sup>O are limited to the two reactions (i) and (ii). Similarly, for <sup>17</sup>O the two dominant reactions are

$$ <sup>16</sup>O + <sup>16</sup>O<sup>17</sup>O \rightarrow <sup>16</sup>O<sup>16</sup>O<sup>17</sup>O \quad (29) $$

$$ <sup>17</sup>O + <sup>16</sup>O<sup>16</sup>O \rightarrow <sup>16</sup>O<sup>16</sup>O<sup>17</sup>O \quad (30) $$

Applying the mass-dependent slope of 0.52 in the δ plot, we calculate δ<sup>17</sup>O from δ<sup>18</sup>O, and obtain the corresponding α<sub>MD</sub> as

$$ \begin{aligned} \alpha_{6-67,MD} &= 1 + (1.272 - 1) \times 0.52 = 1.141 \\ \alpha_{7-66,MD} &= 1 + (0.807 - 1) \times 0.52 = 0.900 \end{aligned} \quad (31) $$

In a scrambled mixture, the two reactions (i) and (ii) contribute equally in the formation of <sup>18</sup>O<sup>16</sup>O<sup>16</sup>O*. Hence the overall coefficients α<sub>m–6,MD</sub> are

$$ \begin{aligned} \alpha_{7-6,MD} &= \frac{1}{2} \cdot 1.119 + \frac{1}{2} \cdot 0.885 = 1.021 \\ \alpha_{8-6,MD} &= \frac{1}{2} \cdot 1.272 + \frac{1}{2} \cdot 0.807 = 1.039 \end{aligned} \quad (32) $$

which is expressed in the δ notation as δ<sup>17</sup>O = +21‰ and δ<sup>18</sup>O = +39‰, respectively. These purely mass-dependent fractionations are indeed much smaller than the Chem-MIF effect.

Following the detailed calculation by Schinke et al. [29], α<sub>m–6,MD</sub> is supposed to be approximately constant with temperature and may be caused by the rapid rotations of the complex on itself triggering its dissociation.

But we see as well that the δ<sub>m–6,MD</sub> values (i.e. +21‰ and +39‰) remain markedly different from those expected for the isotopic exchange between O<sub>3</sub> and O<sub>2</sub>, for which Kaye [40] obtained δ<sup>17</sup>O = –12‰ and δ<sup>18</sup>O = –24‰, at 300 K. We interpret these differences as a contribution of an additional mass-dependent effect linked to the mass-independent fractionation. We illustrate this interpretation by the following calculations.
<figure><img src="image_7.png" /><figcaption>Fig. 5. The mass-corrected fractionation according to Eq. (35) for ¹⁷O (blue) and ¹⁸O (orange), together with the uncorrected data (green). Same color code for the data points from Morton et al. [38]. The upper lines correspond to data generated from the surface of Bitterova and Schinke [35], the lower curves to Tyuterev's surface [36].</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. The pressure dependence of the isotope fractionation, for ¹⁷O (left panel) and ¹⁸O (right panel). Black points from [38], blue points from [39], green points from [41], and red points from [42]. Pressure is given in a logarithmic scale.</figcaption></figure>

## 3.2. Coupling of MIF and MDF

To compare our numerical data to the experimental results of Morton for which no mass correction was applied, we introduce the mass-dependent isotopic fractionation factor $\alpha^*_{m-8,MD}(T)$ standing for the isotopic exchange between $O_3$ and $O_2$ as calculated by Kaye [40]

$$ \alpha^*_{m-6,MD}(T) = \left( \frac{[^{m}\text{O}]}{[^{16}\text{O}]} \right)_{O_3} / \left( \frac{[^{m}\text{O}]}{[^{16}\text{O}]} \right)_{O_2} \\ = \frac{2}{3} \frac{1 + k_{m-6}(T)}{k_{m-6}(T)} \qquad (33) $$

with $m = 7$ or 8 standing for 17 or 18, respectively, with the following relation with temperature:

$$ k_{7-6}(T) = 1.96 \exp\left(\frac{16.8}{T}\right) \\ k_{8-6}(T) = 1.94 \exp\left(\frac{32}{T}\right) \qquad (34) $$

The overall isotopic fractionation factor for ozone — which we may call $\bar{\alpha}_{m-6}(T)$ — becomes then

$$ \bar{\alpha}_{m-6}(T) = \left( \frac{1}{3} + \frac{2}{3} \eta(T) \alpha_{m-6,MD} \right) \times \alpha^*_{m-6,MD}(T) \qquad (35) $$

Reinhardt et al. [34] have shown that the term 1/3 in Eq. (35) has a negligible mass effect on the lifetime ratio $\tau_{l,6-6m}/\tau_{l,6-66}$ ($\tau_I/\tau_I = 1$); hence this term is not corrected for mass effects. The other term

(2/3) corresponds to the two reactions contributing equally for the formation of $^{m}\text{O}^{16}\text{O}^{16}\text{O}$; hence the mass correction factors $\alpha^*_{m-6,MD}$ and the temperature-dependent Chem-MIF factor $\eta(T)$.

The isotopic fractionation factor $\alpha^*_{m-6,MD}(T)$ applies to all the possible ozone isotopomers $O_3^*$ exchanging with $O_2$; hence, the mass-correction factor in Eq. (35), which applies to the overall complex $O_3^*$.

We take now the resulting $\eta(T)$ of our simulations (Fig. 4), with either potential surface, and apply the mass-dependent corrections of Eq. (35). Due to the mass-dependent functions $\alpha^*$ which multiply the nearly degenerate mass-corrected Chem-MIF curves, we generate in this way different temperature-dependent fractionations. Fig. 5 shows the results of the exercise, together with the experimental data points.

The experimental data seem to show a crossover between fractionations of ¹⁷O and ¹⁸O, which is reproduced by both surfaces. For low temperatures the fractionation of ¹⁷O is larger than that of ¹⁸O, which is explained by the fact that the mass-dependent fractionation factors between $O_2$ and $O_3$ ($\alpha^*$ in Eq. (35)) become commensurate with the mass-independent isotopic fractionation factor $\eta(T)$.

## 3.3. Pressure dependence

In order to estimate $k_M$ we look at available experimental data as published in Refs. [38,39,41,42]. These data may be compared to Eq.
(35), by correcting now $\eta(T)$ as

$$ \eta_M(p, T) = \frac{1/\bar{\tau}_I(T) + k_M[M]}{1/\bar{\tau}_D(T) + k_M[M]} \quad (36) $$

where we take the corresponding data for $\bar{\tau}_I$ and $\bar{\tau}_D$ from Tables 2 and 3. From the ideal gas we have $[M](p, T) = n/V = p/(RT)$. If we trace

$$ \bar{\alpha}_{m-6}(p, T) = \left( \frac{1}{3} + \frac{2}{3} \times \frac{(1/\bar{\tau}_I(T)) + k_M[M](p, T)}{(1/\bar{\tau}_D(T)) + k_M[M](p, T)} \alpha_{m-6,MD} \right) \alpha_{m-6,MD}^* (T) \quad (37) $$

for $T = 300$ K and varying pressure $p$, we still have $k_M$ as free parameter. Taking $k_M$ as $10^{-8.9} = 1.26 \times 10^{-9} \text{cm}^3\text{s}^{-1}$, we arrive at a figure following best the trend of the experimental data, see Fig. 6.

## 4. Conclusions

In this paper we present a model with the correct asymptotics, capable of producing a mass-independent isotope fractionation, solely based on the fundamental quantum-mechanical principle of the indistinguishability of identical particles. Particles are here complete atoms, not single elementary particles. That complex quantum particles obey this quantum-mechanical principle has been shown for instance for nuclear collisions [43].

Chem-MIF effects are in general much larger than mass-dependent fractionations – not only we find this by our simulations, but as well from the experimental data with a little statistics and by subtracting the MIF effect. This leads immediately to still puzzling, but experimentally observed orderings of data sets.

Despite the simplicity of our simulations, we were able to correctly reproduce temperature and pressure dependence of the isotope fractionation in ozone. We expect that a fully quantum-mechanical calculation, including right from the beginning the principle of indistinguishability of identical particles, corroborates our lifetime-based approach.

### CRediT authorship contribution statement

**Peter Reinhardt:** Writing – review & editing, Visualization, Software, Resources, Data curation. **François Robert:** Writing – original draft, Validation, Project administration, Methodology, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: We would like to avoid Chr. Janssen or Reinhardt Schinke as referees for this article.

### Acknowledgments

The authors thank R. Schinke and V. Tyuterev for having made available their potential surfaces. Computer time was provided by the Laboratoire de Chimie Théorique at Sorbonne University in Paris, France.

## Appendix A. Computational details

Simulations are carried out with a simple, home-developed code for classical trajectories, using the Verlet algorithm [44] with a fixed timestep of $10^{-4}$ ps. As total energy should be a constant of motion we scale at each timestep the kinetic energy of all particles with a common factor in order to maintain the initial total energy. Temperature is introduced via the distribution of initial conditions following Boltzmann statistics for discrete rotational energies of the molecule and continuous translational energies of the atom.

The initial distance of the atom from the molecule is set to 100 bohr, and a trajectory is run until the intermediate of the three interatomic distances is again at that value. Calculations are carried out on parallel LINUX machines (32–96 processes) using OpenMPI. The simulation program is written entirely in Fortran90, using the IRPF90 preprocessor [45].

## Appendix B. Derivation of $\alpha_{17-16}$

The statistically expected abundance ratio of $^{17}$O and $^{16}$O in an $O_3$ sample is just $x_{17}/x_{16}$.

The production rate constant of an isotopomer A-B-C is $k_f k_M [M] / k_D$ (see Eq. (6)), where we distinguish two different dissociation constants, $k_D(D)$ and $k_D(I)$, standing for “distinguishable” (D) ozone molecules A-X-B, and “indistinguishable” (I) ones, A-X-A. The dissociation constants are inverse lifetimes, $K_D(D) = 1/\tau_D$ and $K_D(I)$.

The overall production rate of $^{17}$O in an ozone sample is then the sum of the individual rate constants of each isotopomer, weighted with the statistical collision probability and the number of $^{17}$O isotopes involved. In the final division the common factor $k_f k_M [M]$ drops out. We assume that isotopes are statistically distributed among single atoms and $O_2$ molecules.

We have then

$$ \begin{aligned} d(^{17}\text{O})/dt &= k_f k_M [M] \left( \underbrace{2 \tau_D x_{16}^2 x_{17}}_{16-17-16} + \underbrace{\tau_I x_{16}^2 x_{17}}_{16-17-16} + 4 \tau_D x_{17}^2 x_{16} + \underbrace{2 \tau_I x_{17}^2 x_{16} + 3 \tau_I x_{17}^3}_{17-16-17} + \underbrace{2 \tau_D x_{16} x_{17} x_{18}}_{16-17-18} + \underbrace{2 \tau_D x_{16} x_{17} x_{18} + \underbrace{2 \tau_D x_{16} x_{17} x_{18} + \underbrace{2 \tau_D x_{16} x_{17} x_{18} + \underbrace{\tau_I x_{17}^2 x_{18}}_{17-16-18} + \underbrace{4 \tau_D x_{17}^2 x_{18}}_{17-17-18} + \underbrace{2 \tau_I x_{17}^2 x_{18}}_{17-18-17} + \underbrace{2 \tau_D x_{18}^2 x_{17}}_{17-18-18} + \underbrace{\tau_I x_{18}^2 x_{17}}_{18-17-18}}_{16-18-17} \right) \\ &= k_f k_M [M] x_{17} \{ \tau_I (x_{16}^2 + 2x_{16}x_{17} + 3x_{17}^2 + 2x_{17}x_{18} + x_{18}^2) + 2\tau_D (x_{16}^2 + 2x_{16}x_{17} + 3x_{16}x_{18} + 2x_{17}x_{18} + x_{18}^2) \} \end{aligned} \quad (B.1) $$

Dividing by the same expression with $x_{17}$ and $x_{16}$ exchanged and dividing the ratio by the statistically expected abundances leads to Eq. (19) in the text.

## Appendix C. Error estimation

We calculate $\eta$ via a purely statistical process — all trajectories are completely independent. This allows on one hand for an embarrassingly parallel code, but asks on the other hand for the calculation of an error bar. The problem here lies in the fact that we cannot calculate $\eta$ from one single trajectory, like a local energy in a Quantum Monte-Carlo simulation. $\eta$ is composed of several averages, and not every trajectory contributes to every average.

For the average lifetime we may apply the standard formula

$$ \begin{aligned} \bar{X}_n &= \frac{1}{n} \sum_{i=1}^{n} X_n \\ \sigma(X)_n &= \frac{1}{\sqrt{n-1}} \sqrt{\bar{X}^2 - \bar{X}^2} \end{aligned} \quad (C.1) $$

We may thus calculate $\bar{X}_{(n+1)}$ from $\bar{X}_n$ as

$$ \bar{X}_{(n+1)} = \frac{n}{n+1} \bar{X}_n + \frac{1}{n+1} X_{(n+1)} \quad (C.2) $$

This applies to $\bar{\tau} = \bar{\tau}_D$, $\bar{\tau}_R$ and $\bar{\tau}_{NR}$. If we assume that $X_R$ and $X_{NR}$ have no error, we may calculate the errors of the lifetime averages and continue as

$$ \begin{aligned} \Delta\eta &= \eta \left( \frac{\Delta\bar{\tau}_D}{\bar{\tau}_D} + \frac{\Delta\bar{\tau}_I}{\bar{\tau}_I} \right) \\ &= \left( \frac{\Delta\bar{\tau}_D}{\bar{\tau}_D} + \frac{X_R \Delta\bar{\tau}_R + X_{NR} \Delta\bar{\tau}_{NR}}{X_R \bar{\tau}_R + X_{NR} \bar{\tau}_{NR}} \right) \end{aligned} \quad (C.3) $$
## Appendix D. The Bitterova-Schinke and Tyuterev potential surfaces

Bitterova and Schinke published a surface as a number of points in a 3-dimensional grid, together with an asymptotic extrapolation procedure based on analytical functions outside the grid. Stored are all combinations of first-order derivatives in each direction [46], i.e. $f(x, y, z)$, $\partial f(x, y, z)/\partial x$, ..., $\partial^2 f(x, y, z)/(\partial x \partial y)$, ..., up to $\partial^3 f(x, y, z)/(\partial x \partial y \partial z)$. The 3-dimensional interpolation is then a cubic spline in all three directions, given by the 4 characteristic cubic polynomials $u_1(x)$, ..., $u_4(x)$.

$f(x, y, z)$ is then the linear combinations of the polynomials times the corresponding derivatives at the corner points of the cube around the point $(x, y, z)$.

The data are furnished in one file, together with a routine for initialization, and a second one for the evaluation of the potential.

Tyuterev's surface is a long analytical expansion in exponentials, sines and cosines in hyperspherical coordinates, assuring the correct permutational symmetry between the three coordinates $\rho, \theta$ and $\phi$. All data are described in the two cited publications [36,37], and the potential has been given to the authors as one single Fortran90 subroutine, without additional data.

## Appendix E. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.chphi.2026.101022.

### Data availability

Data will be made available on request.

### References

[1] J. Bigeleisen, The relative reaction velocities of isotopic molecules, J. Chem. Phys. 17 (1949) 675-678.

[2] H.C. Urey, The thermodynamic properties of isotopic substances, J. Chem. Soc. (Res.) (1947) 562-581.

[3] E.D. Young, A. Galy, H. Nagahara, Kinetic and equilibrium mass-dependent isotope fractionation laws in nature and their geochemical and cosmochemical significance, Geochim. Cosmochim. Acta 66 (2002) 1095-1104.

[4] M.H. Thiemens, J.E. Heidenreich, The mass-independent fractionation of oxygen – a novel isotope effect and its possible cosmochemical implications, Science 219 (1983) 1073-1075.

[5] J.E. Heidenreich, M.H. Thiemens, A non-mass-dependent oxygen isotope effect in the production of ozone from molecular oxygen – the role of molecular symmetry in isotope chemistry, J. Chem. Phys. 84 (1986) 2129-2136.

[6] N. Asset, M. Chaussidon, G. Lombardi, J. Villeneuve, R. Tartèse, S. Mostefaoui, F. Robert, Mass-independent fractionation of oxygen isotopes during high-temperature condensation in cosmochemical plasmas, Proc. Nat. Acad. Sci. 122 (2025) e2426711122.

[7] M.H. Thiemens, S. Chakraborty, G. Dominguez, The physical chemistry of mass-independent isotope effects and their observation in nature, Ann. Rev. Phys. Chem. 63 (2012) 155-177.

[8] S. Chakraborty, M. Ahmed, T.L. Jackson, M.H. Thiemens, Experimental test of self-shielding in vacuum ultraviolet photodissociation of CO, Science 321 (2008) 1328-1331.

[9] F. Robert, M. Chaussidon, A. Gonzalez-Cano, S. Mostefaoui, Oxygen and magnesium mass independent fractionation induced by chemical reactions in plasma, Roc. Nat. Acad. Sci. 118 (2021) 52.

[10] F. Robert, R. Tartese, G. Lombardi, P. Reinhardt, M. Roskosz, B. Doisneau, Z. Deng, M. Chaussidon, Mass-independent fractionation of titanium isotopes and its cosmochemical implications, Nat. Astron. 4 (2020) 762-768.

[11] HBao J. Farquhar, M.H. Thiemens, Atmospheric influence of earth's earliest sulfur cycle, Science 289 (2000) 756-758.

[12] A.L. Buchachenko, Magnetic isotope effect: nuclear spin control of chemical reactions, J.Phys. Chem.A 105 (2001) 996-1011.

[13] A.L. Buchachenko, Mass-independent isotope effects, J.Phys. Chem.B 117 (2013) 2231-2238.

[14] J. Bigeleisen, Second-order correction to the bigeleisen-mayer equation due to the nuclear field shift, Proc. Nat. Acad. Sci. 98/99 (1998) 4808-4809.

[15] J.R. Lyons, E.D. Young, CO self-shielding as the origin of oxygen isotope anomalies in the early solar nebula, Nature 435 (2005) 317-320.

[16] T.R. Rao, G. Guillon, S. Mahapatra, P. Honvault, Huge quantum symmetry effect in the o + $\bar{o}_2$ exchange reaction, J.Phys. Chem.Lett. 6 (2015) 633-636.

[17] H. Yurimoto, K. Kuramoto, Molecular cloud origin for the oxygen isotope heterogeneity in the solar system, Science 305 (2004) 1763-1766.

[18] B.C. Hathorn, R.A. Marcus, An intramolecular theory of the mass-independent isotope effect for ozone, J. Chem. Phys. 111 (1999) 4087-4100.

[19] B.C. Hathorn, R.A. Marcus, An intramolecular theory of the mass-independent isotope effect for ozone. II. Numerical implementation at low pressures using a loose transition state, J. Chem. Phys. 113 (2000) 9497-9509.

[20] G.I. Gellene, An explanation for symmetry-induced isotopic fractionation in ozone, Science 274 (1996) 1344-1346.

[21] Y.Q. Gao, W.C. Chen, R.A. Marcus, A theoretical study of ozone isotopic effects using a modified ab initio potential energy surface, J. Chem. Phys. 117 (2002) 1536-1543.

[22] Y.Q. Gao, R.A. Marcus, Strange and unconventional isotope effects in ozone formation, Science 293 (2001) 259-263.

[23] Ch. Janssen, J. Guenther, K. Mauersberger, D. Krankowsky, Kinetic origin of the ozone isotope effect: a critical analysis of enrichments and rate coefficients, Phys. Chem.Chem.Phys. 3 (2001) 4718-4721.

[24] Ch. Janssen, J. Guenther, K. Mauersberger, D. Krankowsky, Relative formation rates of $^{50}\text{O}_3$ and $^{52}\text{O}_3$ in $^{16}\text{O}-^{18}\text{O}$ mixtures, J. Chem. Phys. 111 (1999) 7179.

[25] R. Schinke, P. Fleurat-Lessard, The effect of zero-point energy differences on the isotope dependence of the formation of ozone: A classical trajectory study, J. Chem. Phys. 122 (2005) 94317.

[26] K. Mauersberger, B. Erbacher, D. Krankowsky, J. Guenther, R. Nickel, Ozone isotope enrichment: Isotopomer-specific rate coefficients, Science 283 (1999) 370-372.

[27] K. Mauersberger, D. Krankowsky, C. Janssen, R. Schinke, Assessment of the ozone isotope effect, in: H. Bederson (Ed.), Advances in Atomic, Molecular, and Optical Physics, vol. 50, 2005, pp. 1-54.

[28] J.M. Carlstad, K. Boering, Isotope effects and the atmosphere, Ann.Rev. Phys. Chem. 74 (2023) 439-465.

[29] R. Schinke, P. Fleurat-Lessard, S.Y. Grebenshchikov, Isotope dependence of the lifetime of ozone complexes formed in o + $\bar{o}_2$ collisions, Phys. Chem. Chem. Phys. 5 (2003) 1966.

[30] M.V. Ivanov, D. Babikov, On molecular origin of mass-independent fractionation of oxygen isotopes in the ozone forming recombination reaction, Proc. Nat. Acad. Sci. 110 (2013) 17708-17713.

[31] D. Babikov, B.K. Kendrick, R.B. Walker, R.T. Pack, P. Fleurat-Lessard, R. Schinke, Formation of ozone: Metastable states and anomalous isotope effect, J. Chem. Phys. 119 (2003) 2577-2589.

[32] T.A. Baker, G.I. Gellene, Classical and quasi-classical trajectory calculations of isotope exchange and ozone formation proceeding through O+O-2 collision complexes, J. Chem. Phys. 117 (2002) 7603-7613.

[33] P. Reinhardt, F. Robert, Mass-independent isotope fractionation in ozone, Earth Planet. Sci. Lett. 368 (2013) 195-203.

[34] P. Reinhardt, F. Robert, On the mass independent isotopic fractionation in ozone, Chem. Phys. 513 (2018) 287-294.

[35] R. Siebert, R. Schinke, M. Bitterova, Spectroscopy of ozone at the dissociation threshold: quantum calculations of bound and resonance states on a new global potential energy surface, Phys. Chem. Chem. Phys. 3 (2001) 1795.

[36] V. Tyuterev, R. Kochanov, S. Tashkun, F. Holka, P. Szalay, New analytical model for the ozone electronic ground state potential surface and accurate ab initio vibrational predictions at high energy range, J. Chem. Phys. 139 (2013) 134307, http://dx.doi.org/10.1063/1.4821638.

[37] D. Lapierre, A. Alijah, R. Kochanov, V. Kokouline, V. Tyuterev, Lifetimes and wave functions of ozone metastable vibrational states near the dissociation limit in a full-symmetry approach, Phys. Rev. A 94 (2016) 042514, http://dx.doi.org/10.1103/PhysRevA.94.042514.

[38] J. Morton, J. Barnes, B. Schueler, K. Mauersberger, Laboratory studies of heavy ozone, J. Geophys. Res. -Atmosp. 95 (1990) 901-907.

[39] M.H. Thiemens, T. Jackson, Pressure dependency for heavy isotope enhancement in ozone formation, Geophys. Res. Lett. 17 (1990) 717-719.

[40] J.A. Kaye, D.F. Strobel, Enhancement of heavy ozone in the earth's atmosphere? J. Geophys. Res. 88 (1983) 8447-8452.

[41] K.L. Feilberg, A.A. Wiegel, K.A. Boering, Probing the unusual isotope effects in ozone formation: Bath gas and pressure dependence of the non-mass-dependent isotope enrichments in ozone, Chem. Phys. Lett. 556 (2013) 1-8.

[42] S. Chakraborty, Anomalous Fractionation of Oxygen Isotopes in Photochemical Reactions (Ph.D. thesis), Ahmedabad, India, 2004.

[43] G.-R. Plattner, I. Sick, Coherence, Interference and the Pauli principle: Coulomb scattering of carbon from carbon, Eur. J. Phys. 2 (1981) 109-113.

[44] L. Verlet, Computer experiments on classical fluids. I. Thermodynamical properties of lennard-jones molecules, Phys. Rev. 159 (1967) 98-103.

[45] A. Scemama, IRPF90: a programming environment for high performance computing, 2009, arXiv:0909.5012v1.

[46] see for instance, W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in FORTRAN, Cambridge University Press, Cambridge (USA, 1992, p. 119.