G Remus

R M Morard

A Bolis

Benuzzi-Mounaix

J Bouchet

F Bottin

V Recoules

F Remus

G Morard

R M Bolis

A Benuzzi-Mounaix

Ab initio calculations of the B1-B2 phase transition in MgO

Ab initio calculations of the B1-B2 phase transition in MgO

I. INTRODUCTION

As one of the major component of the Earth's lower mantle, MgO has been extensively studied experimentally and theoretically. Since the NaCl (B1) structure is stable up to several megabars and thousand degrees Kelvin, MgO is also considered as an internal pressure standard and considerable efforts have been made to build its P -V -T equation of state (EOS). With the discovery of super-Earths it is also crucial to extend the thermodynamic data of MgO up to pressures and temperatures expected in these massive exoplanets 1 , but also for the rocky core of gas and ice giants 2,3 . It has been theoretically predicted for a long time that MgO undergoes a phase transition between the rocksalt structure to the caesium chloride (B2) structure in the 400-600 GPa range [4][5][6][7][8][9][10][11][12] . Experimentally the transition has only been directly observed recently by dynamic X-ray diffraction measurements 13 , and its probable signature in laserdriven decaying 14,15 and steady shocks experiments 12,16 has been evidenced.

At 0 K, the more recent ab initio calculations are in close agreement and therefore they have narrower the range of the transition pressure between 475-510 GPa [9][10][11][12] , see Fig. 1. The small differences in the results can be incriminated to the flavor of the functional used to describe the exchange correlation term : LDA 7,9,17 or GGA (PBE 10,18 ,AM05 11,12,19 ). At high temperature, the disagreement on the location of the phase boundary dramatically increases, from 7% at 0 K, up to 40 % at 8000 K between the lowest transition pressure (∼ 310 GPa 12 ) and the highest one (∼ 440 GPa 10 ), see Fig. 1. One of the goal of this paper is to understand these discrepancies and to propose a well-established phase transition boundary.

All the methods used to compute the phase transition line between the B1 and B2 structures are based on a decomposition of the Helmholtz free energy between the 12 (black), A. B. Belonoshko et al. 9 (blue) and D. Cebulla et al. 11 (red). The black and red lines include anharmonic corrections beyond the quasiharmonic approximation. The green squares are the results of B. Boates et al. 10 using velocity autocorrelation function and the green line a Kechin fit to these data. The inset shows the renormalisation of the thermal transition pressures by the athermal one. static part, or cold curve, obtained at 0 K and the ionic or vibrational contribution coming from the atomic motions in temperature. Since the methods are in good agreement for the transition pressure at 0 K, the differences are due to the way of obtaining the vibrational part.

The key ingredient of this contribution is the phonon spectra, or rather the number of phonon states by energy, the phonon density of states (PDOS) that can be used to calculate the vibrational free energy. The phonon spectra can be obtained at 0 K using either the density functional perturbation theory [20][21][22] (DFPT) or the finite displacements method 23 . Based on these calculations performed at several volumes, the so called quasiharmonic approximation 24 (QHA) introduces a temperature dependence of the phonon frequencies through thermal expansion that alters the volume of the solids. Applied to MgO 6,9 , the QHA predicts the lowest phase boundary and favors the B2 phase over the B1 phase at high temperature with a close to zero Clausius-Clapeyron slope above 8000 K, see Fig. 1. Despite its usefulness and wide application, the QHA suffers of two major flaws: (i) it cannot be applied if the crystal structure is dynamically unstable at 0 K, as the bcc phase of Zr 25 , Li 26 , or the actinides 27,28 and (ii) it neglects the intrinsic thermal effects in opposition to the extrinsic volume dependence.

To test the robustness of the QHA, A. Glensk et al. 29 have recently studied the anharmonicity of several fcc metals using upsampled thermodynamic integration 30 . They show the importance of the anharmonic effects, even at temperatures far below the melting point. More specifically, the validity domain of the QHA for MgO has been studied by A. Erba et al. 31 and Z. Wuet al. 32 by a carefull comparison with experimental data for several thermodynamic quantities. At room pressure, they estimate that the QHA is valid up to about 1000 K, and up to 3000 K at 100 GPa (around 30-40% of the melting temperature). This domain of validity of the QHA has to be taken cautiously, particularly when looking at phase transitions where differences of a few meV in the Gibbs free energy can induce a difference of several GPa on the transition pressure. Basically the QHA will fail if the phonon frequencies have a strong temperature dependence, i.e when atomic vibrations increase in magnitude compared to the atomic distance, as for example close to the melting point but also at a phase transition where soft modes can appear.

To correct the QHA, D. Cebulla et al. 11 and S. Root et al. 12 have used ab initio molecular dynamics (AIMD) but their methods still relies on the QHA to deal with temperatures where anharmonicity has to be included in the estimation of the free energy. Nonetheless, their results show that the QHA overestimates the stability domain of the B2 phase in temperature and that the anharmonic effects drive to higher transition pressures at higher temperatures, see Fig. 1. Using Fourier analyzing of the velocity autocorrelation function from AIMD B. Boates et al. 10 observe a similar effect but with a stronger intensity at 10000 K. Although this method should in principles include the whole anharmonicity present in the AIMD it requires long simulations and large supercells to avoid finite-size effects, as when employed with classical molecular dynamics. In Ref. [10], they used supercells of 64 and 128 atoms to simulate the B1 and B2 phases, rather small numbers compared to the several thousand atoms used in classical MD.

In the past ten years, strong efforts have been made to take into account explicit temperature effects and major advances have been obtained. New methods capturing the thermal properties of solids at non-zero temperature are now available and can be applied with ab initio calculations. These approaches combine ideas including finite large displacements, molecular dynamics sampling, self consistent harmonic theories, and different force fitting schemes. The most widely used methods include Self-Consistent Ab Initio Lattice Dynamics (SCAILD) 25 , Stochastic Self-Consistent Harmonic Approximation (SSCHA) 33,34 , Temperature Dependent Effective Potential (TDEP) 26,35,36 , Anharmonic LAttice MODEl (ALAMODE) 37 , Compressive Sensing Lattice Dynamics 38 . Other methods obtain anharmonic contributions via a derivation of the Gibbs energy 29 , or a series expansion of the interatomic forces constants 39 . A large number of new phenomena, intrinsically temperature dependent, can now be captured: the modification of the PDOS and free energy, the (T,P) phase transition boundaries, the evolution of elastic constants or Grüneisen coefficients, the phonon lifetimes, the thermal conductivity...

Here we have used the TDEP technique developed by O. Hellman and coworkers 26,35,40 to extract the vibrational frequencies from AIMD simulations, with an original implementation in the ABINIT package [START_REF]The ABINIT code is a common project of the Catholic University of Louvain (Belgium), Corning Incorporated, CEA (France) and other collaborators[END_REF][START_REF] Gonze | [END_REF] . This method takes roots in the pionner work of Hooton 43 and K. Esfarjani et al. 39 and has already been successful to describe phase transitions at high temperature as for Ti and Zr 26,44 , 4 He 35 or U and Pu 27,28 . In TDEP, a Taylor expansion on the atomic displacements is used to fit the Born-Oppenheimer molecular dynamics potential energy surface at finite temperature and to obtain the interatomic force constants (IFC). The second order IFCs give the phonon frequencies while the higher order terms can be used to calculate phonon lifetimes, thermal conductivity or the Grüneisen parameter.

The paper is structured as follows. In Sec. II we give the details of the AIMD simulations and the construction of the free energy. In Sec. III we give our results for the B1-B2 transition curve and we compare our results with previous calculations and recent experiments. This is followed by the presentation of several thermodynamics quantitites including thermoelasticity, Grüneisen parameter and wave velocities.

II. SIMULATION DETAILS

Simulations were performed using the ABINIT package [START_REF]The ABINIT code is a common project of the Catholic University of Louvain (Belgium), Corning Incorporated, CEA (France) and other collaborators[END_REF][START_REF] Gonze | [END_REF] in the framework of the Projector Augmented Wave (PAW) method 45,46 and by means of the local density approximation (LDA) 17 . Using ATOM-PAW [47][START_REF]ATOMPAW is a general license public code developed at Wake Forest University. Some of its capabilities have been developed at the Commissariat à l'énergie atomique[END_REF][START_REF] Dewaele | [END_REF] , we generated a PAW atomic data with a radius r PAW equals to 1.0 Å, with 2s, 2p, 3s, and 3p states as valence electrons (10 electrons) for Mg and a radius r PAW equals to 0.6 Å, with 2s and 2p states as valence electrons (6 electrons) for O. The cutoff energy chosen for the plane wave set along the simulations is 816 eV. To determine the B1-B2 phase boundary we need to calculate the free enthalpy G of both phases, with G(T, P ) = F (T, P ) + P V (T, P ), (1) where the free energy F is given by

F (T, V ) = E 0 (V ) + F vib (V, T ) + F el (V, T ). (2)
Here E 0 (V ) is the energy of the static lattice, F vib the vibrational free energy due to the atomic motions and F el the electronic free energy. For the static calculations, at T=0K, we use the unit cells of the B1 and the B2 phases with two atoms and a 8 × 8 × 8 Monkhorst-Pack mesh.

To obtain the vibrational part we use two methods, the QHA and the TDEP. They differ in the way of calculating the interatomic force constants (IFCs). The QHA is based on static calculations, so at 0 K, involving only the unit cells and the DFPT. The dynamical matrices are calculated on a grid of q points and then Fourier transformed, yielding to the IFCs. With this method, the phonon frequencies are only dependent of the volume. Our calculations were performed on 4 × 4 × 4 grid of q points in the Brillouin zone. The second method, the TDEP 26,35 , calculates directly the IFC in real space from sets of forces and displacements obtained from AIMD. For this purpose, we use supercells of 128 atoms for the B1 and B2 phases, corresponding to 4 × 4 × 4 rhombohedral B1 unit cells and cubic B2 unit cells (we also use larger and smaller supercells for convergence studies, see below Fig. 4). Note that contrary to previous calculations [10][11][12] , also based on AIMD, we use the same number of atoms for the two structures. Our AIMD were performed with a 2 × 2 × 2 Monkhorst-Pack mesh leading to the inclusion of 4 special k-points. This mesh was checked against a 3 × 3 × 3 Monkhorst-Pack mesh to ensure the required accuracy, less than 1 meV on the vibrational free energy. Simulations were performed in the N V T ensemble (constant number of particles, constant volume and temperature) and were run for about 5 ps using a time step (τ ) of 0.60 fs. Taking into account the whole symmetries of the system, we have to calculate 19 coefficients for the B1 phase (up to the 8 th shell of nearest neighbours) and 21 coefficients for the B2 phase (up to the 9th shell of nearest neighbours) to build the IFC matrices with the TDEP. So, contrary to the DFPT which uses unit cells, the TDEP requires AIMD on supercells of atoms to obtain the IFCs at a sufficiently large distance where they become negligible. It is therefore much more expensive in terms of computational time but goes beyond the QHA since the phonon frequencies are not only volume dependent but also explicitly temperature dependent.

Once the IFC are obtained, by either methods, a Fourier transform is performed to get the dynamical matrix at any q-point of the Brillouin zone and to calculate the PDOS. Then the free energy is given by :

F vib (V, T ) = k B T ∞ 0 dωg(ω(V, T )) ln[2 sinh(hω(V, T )/2k B T )].
(3) where the temperature dependence of the phonon frequencies, ω, and the PDOS, g(ω), is only valid for the TDEP method. Once the free energies are calculated with Eq. 2, they are fitted by a Vinet EOS to obtain their pressure dependence and then the free enthalpy with Eq. 1.

III. B1-B2 PHASE TRANSITION

A. QHA versus TDEP method

We show in Fig. 2 the variation of the PDOS as a function of temperature obtained with TDEP for the B1 and B2 phases at a fixed density of 6.75g/cm 3 . For both structures, the whole spectra softens with temperature. This softening is a signature of explicit thermal effects, and since in QHA the spectrum is independent of the temperature at fixed volume, it will induce lower vibrational free energies in TDEP. This is directly seen in Fig. 3 where we compare the free energies obtained with both methods. Up to about 2000 K, both methods give the same free energies, but at higher temperatures the free energies obtained with TDEP deviate from the QHA. More importantly, the free energies difference, ∆F vib = F B1 vib -F B2 vib is smaller in TDEP, as shown in the inset of Fig. 3. This means that the QHA overestimates the stability domain of the B2 phase compared to the TDEP method. In other words, the intrinsic anharmonic effects will move the phase boundary at higher pressure at high temperature, in agreement with previous works [10][11][12] .

We have also applied the TDEP method on several supercell sizes: 64, 128 and 216 atoms for the B1 phase and 54, 128 and 250 atoms for the B2 phase. The vibrational free energy differences between the B1 and B2 phases, ∆F vib , are shown on Fig. 4 as a function of the supercell sizes. At 2000 K, the effect of the supercell size is negligible, and the TDEP method is close to the QHA for all the volumes considered. At 8000 K we observe a size effect. The smaller supercell sizes overestimate the free energy differences, favoring the B2 phase, while the calculations with 128 atoms and more are in agreement.

We also observe a larger difference between TDEP and QHA for smaller volumes (larger pressure). This seems in contradiction with the belief that anharmonic effects should decrease with pressure along an isotherm. In fact, it is not. For both phases, the difference between QHA and TDEP decreases with the volume, but not in the same ratio : F vib of the B2 phase becomes closer to the QHA value than for the B1 phase.

Our results for the phase transition boundary are presented in Fig. 5 in comparison with the previous the- oretical data of Fig. 1. First, our QHA results are in agreement with the work of A. B. Belonoshko et al. 9 and D. Cebulla et al. 11 with similar slopes of phase boundary (see the inset of Fig. 5). We therefore observe the same discrepancy with the QHA data of S. Root et al. 12 . The phase boundary predicted with the TDEP method is steeper, increasing the stability domain of the B1 phase at high temperature. Consequently, our curve is in agreement with the one obtained by B. Boates et al. 10 using the VACF method and which, in principle, should take into account the full anharmonicity of the simulated system. The previous attempts to go beyond the QHA 11,12 and to include the intrinsic anharmonicity, predicted also an extension of the stability domain of the B1 phase, but to a lesser extent than the one that we obtain with the TDEP. In the next session we analyze these discrepancies and we will try to reconcile, if possible, the different methods.

B. Comparison with previous calculations

D. Cebulla et al. 11 , propose to take into account the anharmonicity by expressing the internal anharmonic energy U anh with :

U anh = U AIMD + E KIN -U vib -E 0 (4)
where U AIMD and E KIN are respectively the internal and kinetic energy of the AIMD, E 0 the electronic ground state energy and U vib is the harmonic internal energy given by : Then the anharmonic free energy at a temperature T f is obtained via a thermodynamic integration from an initial temperature T i :

U vib (V, T ) = ∞ 0 dωg(ω(V, T ))hω(T )/2 coth(hω(V, T )/2k B T ). (5)
F anh (T f ) = F anh (T i ) T f T i -T f T f Ti U anh (T ) T 2 dT (6)
Since F anh (T i ) is not known, the idea is to infer a temperature where this energy is negligible and to eliminate the first term of the right hand side of Eq. 6. This temperature has also to be high enough to neglect the zero point motion energy contained in U vib but not taken into account in U AIMD since the atomic motions follow a classical trajectory. In the same way, the temperature should be high enough so that the kinetic part contained in U vib reaches the classical value of 3k B T given by the AIMD. D. Cebulla et al. 11 chose a value of T i = 3000K. At this temperature and above, since the QHA and the TDEP are based on the harmonic approximation, U vib will be equal in both models to the classical value of 3k B T , and obviously, Eq. 6 will give the same correction. In other terms, this method is completely independent of the shape of the PDOS (see Fig. 2). So it clearly neglects the fact that the phonon frequencies are temperature dependent and it will not correct F vib : this intrinsic anharmonicity will still be missing in F anh . Therefore, the QHA and the TDEP cannot be reconciled with this method. But it can be used to estimate the anharmonic part not included in the TDEP.

To avoid the problem of the zero point motion and the kinetic energy we can use the harmonic potential energy 35 directly calculated with the IFCs obtained with the TDEP to calculate the anharmonic internal energy :

U anh = U AIMD -1/2 ij Φ ij u i u j -E 0 (7)
where the Φ ij are the IFCs and u i the atomic displacements. We present in Fig. 6 the internal anharmonic energy as a function of temperature and for a density of 6.75g/cm 3 . At low temperature U anh converges towards zero, the atomic motions in both phases are purely harmonic. As the temperature increases, anharmonic effects arise and strongly depend of the structure. For the B2 phase, U anh is negative and close to zero, while it strongly increases with temperature for the B1 phase. This shows again that anharmonic effects are more important in the B1 phase than in the B2 one and that even at high pressure, these effects cannot be neglected. Following Eq. 6, we can calculate F anh , see the inset of Fig. 6. Compared to the QHA, the intrinsic anharmonic effects described by the TDEP lower the vibrational free energy difference between the B1 and B2 phases, see Fig. 3, pushing the phase boundary to higher pressure at high temperature.

If we add F anh to F vib obtained with TDEP, it will have a cumulative effect and the phase boundary will steepen even further. For the B1 phase, F anh counts for about 20 % of the anharmonic correction to the QHA, the rest being given by the TDEP, while for the B2 phase, F anh is positive, so slightly brings back the TDEP toward the QHA value by 5%. S. Root et al. 12 propose an other approach to go beyond the QHA based on a thermodynamic integration (TI) to calculate the change in entropy along an isochore :

∆S = T f Ti 1 T ∂U AIMD ∂T V dT (8)
The main drawback of this method is obviously the number of simulations required to obtain the entropy at high temperature. The authors used an increment of 250 K, which means that to have the entropy at 8000 K, starting from 1000 K, about 30 AIMD simulations are necessary. Multiplied by the number of isochores, the computational cost can rapidly become too cumbersome.

We have compared this method with the TDEP one on several isochores, using the same supercell sizes as in Ref. 12, 64 and 54 atoms in the B1 and B2 phases respectively. An example of the entropy differences between the two phases is given in fig. 7 for ρ = 7.75g/cm 3 . In QHA this difference is a straight line, proportional to the temperature. Note that we have neglected the zero point motion and other quantum effects not included in eq. 8, but they have a small impact and only at low temperatures. The agreement between TDEP, eq.8 and QHA is excellent up to 2000 K where we start to observe a deviation. As the temperature increases, the entropy difference between B1 and B2 diminishes compared to the QHA, as already shown in fig. 3. The effect is stronger with eq.8 than with TDEP, meaning that it will push the phase boundary to higher pressures than the ones obtained with TDEP. A result in contradiction with the lower pressures reported by S. Root et al. 12 , see Fig. 5. To understand these discrepancies, we have firstly checked the influence of the exchange correlation functional, since S. Root et al. 12 use the Armiento-Mattson (AM05) 19 while we use the LDA. The vibrational free energies obtained with TDEP are slightly different (few meV), but in favor of the B1 phase with AM05 compared to LDA. Therefore, the discrepancies that we observe compared to Ref. 12 can not be attributed to the choice of the exchange correlation functional. S. Root et al. 12 mentioned deviation with the QHA around 5000 K and 400 GPa, while we already observed an effect at 2000 K. If we use this value of 5000 K as the starting temperature for the TI, we underestimate the entropy differences, as shown by the open circles in fig. 7. The consequence is of course a phase transition at pressures closer to the ones predicted by QHA and therefore by S. Root et al. 12 . We have also performed calculations with bigger supercells, 128 and 216 atoms for the B1 phase, 128 and 250 atoms for the B2 phase. At 2000K the internal energy is converged to 1 meV as a function of the supercell size but at 8000 K we observe differences up to 5 meV. This supercell size effect is similar to the one that we observe for the TDEP method, see Fig. 4, where we have shown that a small supercell favors the B2 phase.

Since, Eq. 8 is based on the internal energy coming from AIMD simulations with no further approximation as in TDEP, it should in principle contained the whole anharmonicity of the simulated system. The difference that we observe between TI and TDEP is close to F anh obtained with the method of D. Cebulla et al. 11 and if we add this correction to the entropy found in TDEP we almost recover the TI, see Fig. 7. This confirms that the anharmonic terms that we do not take into account with TDEP correct to a larger extend the QHA.

B. Boates et al. 10 used the Fourier transforms of the velocity autocorrelation functions (VACF) to calculate the vibrational entropies. They obtained a transition pressure of 405 GPa at 10000 K, close to the value obtained with the TDEP method, see Fig. 5. This method is usually used with classical MD, since it requires large supercells and considerably long simulation time to obtain the PDOS. B. Boates et al. 10 used 64 and 128 atoms only in their supercells for the B1 and B2 phases respectively. We have performed simulations up to 20 ps for several ρ, T points to obtain the VACF and to calculate the vibrational free energies. This lead us to think that the supercell sizes used in Ref. 10 are not large enough to converge the free energy. At 2000 K and 7.75g/cm 3 for example, a temperature where anharmonicity should be negligible, we do not recover the QHA result, due to a strong underestimation of the vibrational free energy of the B1 phase (100 meV), partially compensated by a similar effect on the B2 phase (40 meV), see Fig. 7. At 8000 K, it is slightly better, but again, certainly due to a cancelation of errors between the two phases.

To sum up, we are now able to explain the discrepancies between our calculations and the previous ones. These differences come from the treatment of the intrinsic anharmonic part for D. Cebulla et al. 11 or by an overestimation of the temperature at which the anharmonic effects have to be taken into account for S. Root et al. 12 . Concerning the work of B. Boates et al. 10 , we do not recover their results.

C. Comparison with experiments

Fig. 8 shows the experimental P -T curves with our calculated B1-B2 transition line and Hugoniot. We plot also recent ab initio predictions of the melting curve 9,12 . Recently, two decaying shock experiments 14,15 have studied the phase diagram of MgO. In the first one 14 , a clear bump and a slight slope change were observed in the P -T curve along the Hugoniot. The strong bump at 450 GPa, was interpreted as the B1-B2 phase transition and the second one at 650GPa, much softer, as the B2melt. The second experiment 15 , reported only the first bump, at similar thermodynamic conditions, but did not observe the slope change at higher pressure. Contrary to Ref. 14, they interpreted this feature as the B2-melt and not the B1-B2 solid phase transition, similarly to interpretation of decaying shock experiments performed on SiO2 50 or MgSiO3 51 . In addition, with impact shock compression experiments, S. Root et al. 12 observed only one breack in the Hugoniot around 360 GPa. Note that if a bump along the Hugoniot can be attributed to a phase transition the nature of this phase transition cannot be determined in the absence of microscopic in situ diagnos-tics. Can our calculations help to solve this controversy? Our predicted B1-B2 phase boundary is close to the bump observe in the decaying shock experiments at 450 GPa. Knowing that anharmonic effects not included in our calculations should increase the stability domain of the B1 phase (see the discussion in the previous section), a nice correlation between the bump observed in decaying shock and the B1-B2 phase transition could be deduced. However, different arguments are on the side of a melting observation during decaying shock experiments. Firstly, the low intensity or absence of a second bump related to melting signature in S. McWilliams et al 14 is ambiguous, as a large thermodynamic domain is investigated (up to 1 TPa and 30000 K) far above the predicted melting line 9,10,12,16,52 . Secondly, there is a strong disagreement between the calculated Hugoniot for the B2 phase and the experimental points above 450 GPa, as already noted by D. Cebulla et al. 11 . The theoretical Hugoniot in the B1 domain is consistent with the lower part of the experimental results (below 450 GPa). After the bump, the experimental points are close to our calculated liquid Hugoniot. This comparison between our theoretical and experimental Hugoniots seems therefore to indicate a direct transition between the B1 and the liquid phase along the Hugoniot, close to the triple point, with a melting point around 450 GPa and 10000 K as proposed in Ref [15]. This melting point is close to the computed melting line of the B1 phase of N. de Koker et al. 53 , using a comparison of the Gibbs free energies of the solid and the liquid and based on AIMD, but also with extrapolation of the B1 melting line obtained with the Z method by A. B. Belonoshko et al. 9,54 and with the two phase calculations of S. Root et al. 12 and R. Musella et al. 55 . But in this region, the calculated melting line of the B2 phase is above the results for the B1 phase 12 , see Fig. 8.

Ab-initio calculations of melting curves are difficult and their uncertainties can be around several hundred kelvins. These uncertainties are cumulative when comparing melting lines obtained with different structures. In our simulations, at 6.0 g/cm 3 and 10000 K, corresponding to a pressure of 330 GPa, the B2 structure melts spontaneously (see the red cross on Fig. 8), contrary to the B1 structure, while at the same density and 9000 K (283 GPa), the B2 structure stays solid. This melting of the perfect crystal along an isochore is known as the heat until it melts method and generally induces a superheating that results in an overestimation of the melting temperature 56 . New simulations, with larger supercells and with new methods to determine the melting are necessary to clarify this debate. In parallel, shock or diamond anvil cells experiments, coupled with X-ray diffraction, are strongly needed. Purple, magenta and orange dots : R. M. Bolis et al. 15 , green dots: R. S. McWilliams et al. 14 . The black solid line with filled circles is the B1-B2 boundary obtained with TDEP and the blue line the calculated Hugoniot. The dashed blue line is an extension of our Hugoniot curve for the liquid phase in the stability domain of the solid phase. The red curve is the calculated melt line of S. Root et al. 12 for the B2 structures and the red dashed one is the calculated melt line of N. de Koker et al. 53 for the B1 structure. The red cross indicates the melting of an AIMD with B2 as the starting structure of the simulation.

IV. THERMODYNAMICS PROPERTIES

One of the main advantage of using the TDEP method is that we can not only calculate a temperature dependent free energy but also extract from the AIMD various thermodynamical quantities as the thermoelasticity, the Grüneisen parameters or the thermal expansion.

A. Thermoelasticity

The elastic constants can be obtained from the secondorder force constants following [START_REF] Liebfried | Solid State Physics[END_REF][START_REF] Wallace | Thermodynamics of Crystals[END_REF][START_REF] Shulumba | [END_REF] :

C αβγδ = A αγβδ + A βγαδ -A αβγδ (9)
with,

A αβγδ = 1 2V ij Φ αβ ij d γ ij d δ ij ( 10
)
where d α ij is the αth component of the distance between the atom i and j. The full elastic constants matrix, C αβγδ , can be reduced following the symmetries and the Voigt notation. The temperature effects are then directly introduced in the elastic constants by the temperature variation of the IFCs. As shown by N. Shulumba et al. [START_REF] Shulumba | [END_REF] , this method is much faster and practical than the one using deformation matrix which needs several AIMD simulations or based on the QHA to calculate the free energies from the phonon dispersions of the strained lattices 60 .

The main drawback of this method is that since the elastic constants are related to phonon frequencies in the long wavelenght limit, it requires large supercells to converge the sum in Eq. 10. To circumvent this difficulty, we only retain the temperature dependence of Eq.10 that we apply to the static calculations at 0 K 59 .

We present on Fig. 9, the pressure and temperature evolution of the elastic constants (C 11 , C 12 and C 44 ), the isotropic bulk (K T ) and shear (G) moduli of MgO for the B1 and B2 phases, calculated using the TDEP method. K T and G were obtained using the Voigt-Reuss-Hill average. For our 0 K results, we use the stress-strain relations and a deformation matrix applied on the cubic unit cell of the B1 and B2 phases. We also plot on Fig. 9 the results of Karki et al. 60,61 at 0 K and at 2000 K using the QHA (for the B1 phase). For both phases we observe a softening of C 11 and C 44 while C 12 hardens with temperature. For the B1 phase, C 44 shows almost no pressure dependence above 150 GPa and the effect of the temperature is also very weak. In the B2 phase, C 44 has a different behavior, it increases with pressure and softens with temperature. Since K T is given by (C 11 + 2C 12 )/3, the effect of temperature on these two elastic constants cancel each other to result in a small decrease of the bulk modulus in temperature for both phases. The shear modulus G evolves in temperature in a similar way than C 44 .

B. Grüneisen parameter

The Grüneisen parameter γ provides a bridge between P and T and is a key quantity for the Earth's interior where P as a function of depth is well constrained, thanks to seismological data, but T is still poorly known. The evolution of γ with density, but also with temperature is therefore essential to build reliable EOS for planetary modeling 62,63 . The mode Grüneisen parameters are defined by the volume derivatives of the phonon-mode frequencies :

γ s = - V ω s ∂ω s ∂V T (11)
They can be obtained from the third-order force constants Ψ ijk following 40 :

γ s (q) = - 1 6ω 2 s (q) ijkαβγ Ψ αβγ ijk X α is (q)X β js (q) M i M j r γ k exp [iq.R j ] (12
) where M i and r α i are the mass and the αth component of the vector position of atom i, ω s the frequency of mode s, the X's are the eigenvectors and R i the lattice vector of the unitcell of atom i.

From the γ s 's the thermodynamical Grüneisen parameter is obtained by :

γ = 3Na s=1 γ s C V,s C V (13)
where the C V,s are the mode heat capacities, and C V the specific heat.

In EOS, a commonly used expression for the volume dependence of the Grüneisen parameter has been given by Al'tshuler et al. 64 :

γ = γ ∞ + (γ 0 -γ ∞ )x β (14)
where γ 0 is the Grüneisen parameter at ambient conditions, γ ∞ is the Grüneisen parameter at infinite compression (x = V /V 0 = 0) and β is a fitted parameter. Since MgO has a long history as an internal pressure standard at extremely high pressure and temperature there is a large collection of parameters for the volume dependance of its Grüneisen parameter. The Mie-Grüneisen model have been used to reproduce experimental P -V -T data and obtain the parameters entering Eq. 14.

We present in Fig. 10 several versions of Eq. 14 based on experimental data 8,65-68 and the results of A. Oganov et al. 69 based on theoretical P -V -T points obtained with AIMD, while the corresponding fitting parameters are given in Table I. We also report our values obtained with Eq. 13 for the B1 and B2 phases at 2000 and 8000 K. We have performed an AIMD at 300 K and at the equilibrium volume of the B1 phase. We obtain a value 66 . For the B2 phase we find larger values for γ, about 25 %, compared to the B1 phase. We observe a similar difference in the thermal pressure extracted from AIMD simulations. For the B2 phase, the volume dependence of γ is rather weak. We performed an AIMD simulation at 9 g/cm 3 (V /V 0 = 0.4) and we obtain a value of 1.08.

In the Mie-Grüneisen EOS, the temperature dependence of γ at constant V is neglected to directly obtain the thermal pressure in terms of the thermal energy : P TH = γE TH /V . O. L. Anderson 62 has shown that this assumption does not hold for MgO. Using experimental data, he obtained, at constant V , a value of 1.39 at 1800 K, so 10 % lower than at room temperature. At 8000 K, and for larger compression, we also observed a decrease of γ with temperature, see Fig. 10 and 11. The effect is larger for the B1 phase and is almost constant with volume, while for the B2 phase the temperature dependence diminishes with the compression (Fig. 11). This shows the limitation of the Mie-Grüneisen approach to obtain the thermal pressure.

C. Wave velocities

From the Grüneisen parameter we can also obtain the volume thermal expansion coefficient defined as:

α V = γC V K T V (15)
Then, knowing γ and α, we can obtain the adiabatic incompressibility K S = (1 + αγT )K T and with the shear modulus G, the compressional and shear wave velocities, V P and V S following : V p = K S + 4/3G ρ (16)

V s = G ρ (17)
Our results are presented in Fig. 12. In both phases, V P has a similar pressure and temperature dependence while V S shows a larger pressure and temperature dependence in the B2 phase that are mainly due to the differences observed in the evolution of the elastic constant C 44 between B1 and B2, see for V P and V S in iron at 6000 K (before the premelting effect 71 ) and 360 GPa while we found 4.5 and 9.7 % for the same conditions.

V. CONCLUSIONS Using a recent method to calculate the vibrational contribution to the free energy we have studied the phase transition between the B1 and B2 phase in MgO. Compared to previous theoretical works, we obtain a stronger correction to the quasiharmonic approximation at high temperature which pushes the transition line to higher pressures. We show that anharmonic effects cannot be neglected and have to be carefully calculated. We have also estimated the anharmonic part not described by our method and we show that its effect is to further increase the transition pressure. Superimposed on recent shock experiments, our results indicate a direct transition between the B1 structure towards the liquid along the Hugoniot close to the triple point B1-B2-liquid. New experiments using x-ray diffraction are necessary to confirm this assertion. Using the full capabilities of the TDEP we also derive several thermodynamic quantitites as the elastic constants, the Grüneisen parameters and the sound velocities that are crucial quantities to build reliable high pressure high temperature EOS to describe planetary interiors.

FIG. 1 .

 1 FIG. 1. (Color online) Theoretical B1-B2 phase boundary of MgO. Dotted lines : quasiharmonic approximation of S. Root et al. 12 (black), A. B. Belonoshko et al. 9 (blue) and D. Cebulla et al. 11 (red). The black and red lines include anharmonic corrections beyond the quasiharmonic approximation. The green squares are the results of B. Boates et al. 10 using velocity autocorrelation function and the green line a Kechin fit to these data. The inset shows the renormalisation of the thermal transition pressures by the athermal one.

FIG. 2 .

 2 FIG. 2. (Color online) PDOS of the B1 and B2 phases of MgO at 6.75g/cm 3 and for temperatures of 1000, 5000 and 8000 K.

FIG. 3 .

 3 FIG. 3. (Color online)Comparison of the vibrational free energies of the B1 (in red) and B2 (in blue) phases of MgO obtained with DFPT (lines) and TDEP (filled circles) at 7.75g/cm 3 . Inset : difference between the free energies of the B1 and the B2 phases in QHA (solid line) and with the TDEP method (filled circles and dotted line)

FIG. 5 .

 5 FIG. 5. (Color online) Same as Fig. 1, with our QHA results (purple open circles and dotted line) and TDEP results (purple filled circles and solid line).

FIG. 6 .

 6 FIG.6. (Color online) Internal anharmonic energy U anh for the B1 (in black) and B2 (in red) phases as a function of temperature obtained with Eq. 7 at 6.75g/cm 3 . Inset : F anh as a function of temperature.

FIG. 7 .

 7 FIG. 7. (Color online) Entropy difference between the B1and B2 phases at 7.75g/cm 3 as a function of the temperature. The solid line corresponds to the QHA, the filled black circles to the TI of Eq. 8 and the red squares to TDEP. The blue squares are similar to the red ones with the correction F anh of Eq. 6 while the open black circles corresponds to an TI starting at 5000 K. The green triangles are the results of the VACF at 2000 and 8000 K. The error bars have been estimated with respect to calculation with supercell sizes of 216 atoms and 250 atoms for the B1 and the B2 structures respectively.

FIG. 8 .

 8 FIG. 8. (Color online)Experimental shock temperature data in comparison with our predicted B1-B2 phase boundary. Purple, magenta and orange dots : R. M. Bolis et al.15 , green dots: R. S. McWilliams et al.14 . The black solid line with filled circles is the B1-B2 boundary obtained with TDEP and the blue line the calculated Hugoniot. The dashed blue line is an extension of our Hugoniot curve for the liquid phase in the stability domain of the solid phase. The red curve is the calculated melt line of S. Root et al.12 for the B2 structures and the red dashed one is the calculated melt line of N. de Koker et al.53 for the B1 structure. The red cross indicates the melting of an AIMD with B2 as the starting structure of the simulation.

FIG. 9 .

 9 FIG. 9. (Color online) Elastic constants, C11, C12 and C44, isotropic bulk (KT ) and shear (G) moduli of MgO of the B1 and B2 phases as a function of pressure and temperature. The dashed lines are the results of Karki et al. 60,61 at 0 K (black) and at 2000 K (red) using the QHA.

3 FIG. 11 .

 311 FIG. 10. (Color online) Grüneisen parameter of MgO. The black open symbols are the effective Grüneisen parameter of A. R. Oganov et al.69 . The lines are Grüneisen parameters used in semi-empirical equation of state8,65,66,68

Fig. 9 .

 9 Compared to iron, temperature has a weaker effects on the wave velocities. B. Martorell et al. 70 predict a decrease of 8.6 and 29.6%

FIG. 12

 12 FIG.12. (Color online) Compressional and shear wave velocities of MgO as a function of pressure at 0 (black), 2000 (red), 5000 (blue) and 8000 K (green). The kink in the curves indicates the phase transition between B1 and B2.

TABLE I .

 I Fitting parameters of function 14 of the Grüneisen parameter of the B1 phase of MgO. of 1.44 for the Grüneisen parameter, close to the experimental value of 1.54. Using this value as γ 0 we have fitted our results for the B1 phase at 2000 K to obtain γ ∞ and β, see TableI. Our direct calculations of γ are close to the results of A. Oganov et al.69 obtained by fitting the thermal pressure extracted of AIMD simulations performed with GGA. Our results are also in good agreement with the semi empirical EOS proposed by P. I. Dorogokupets et al.

	EOS	γ0	γ∞	β
	This work	1.440	0.800	3.255
	PAW 69	1.455	0.841	3.057
	Dorogokupets et al. 66	1.50	0.75	2.96
	Wu et al. 8	1.520	0.606	1.406
	Tange et al. 68	1.431	1.016	3.50
	Speziale et al. 65	1.524	1.325	11.8

VI. ACKNOWLEDGEMENTS

We would like to thank X. Gonze for helpful discussions. We also thank M. Desjarlais for sharing with us details of the method used in Ref. 12. This research was supported by the PlanetLab program of the Agence Nationale de la Recherche (ANR) grant ANR-12-BS04-0015-04.


François Robert

Lambert Baraut-Guinet

Pierre Cartigny

Peter Reinhardt

An experimental test for the mass independent isotopic fractionation mechanism proposed for ozone

Keywords: Ozone, Oxygen Isotopes, Mass Independent Fractionation, Experimental Observations, Theoretical Predictions

Ozone exhibits large and mass independent isotopic fractionations (MIF) in oxygen isotope ratios relative to molecular di-oxygen (O2) from which it is formed. An interpretation of this effect was proposed based on the behavior of the indistinguishable isotopes 16 O in scattering processes. We report here an experiment aimed at testing one of the predictions of this model.

O3 was formed by high frequency discharge in O2 with pressures ranging between 1.6 and 38 Torr. The isotopic evolution of the closed O2 reservoir was monitored during its distillation taking place during the continuous removal of ozone by condensation. Its composition evolves from a mass independent to a mass dependent fractionation along with the decrease in pressure.

The isotopic pathways defined by this evolution in the 3 isotopes diagram are in quantitative agreement with the theoretical prediction of the disappearance of MIF with the increase of the complex lifetime stabilized as ozone.

Introduction

The discovery of a mass independent oxygen isotope fractionation (MIF) during the synthesis of ozone was found to contain a considerable potential for applications in Earth and Atmospheric Sciences [START_REF] Thiemens | The mass-independent fractionation of oxygen -a novel isotope effect and its possible cosmochemical implications[END_REF][START_REF] Thiemens | Atmosphere science -Mass-independent isotope effects in planetary atmospheres and the early solar system[END_REF][START_REF] Thiemens | Introduction to Chemistry and Applications in Nature of Mass Indepedendent Isotope Effects Special Feature[END_REF] . A comprehensive set of experimental data illustrates this effect [START_REF] Thiemens | The observation of mass-independent fractionation of oxygen in an RF discharge[END_REF][START_REF] Heidenreich | A non-mass-dependent oxygen isotope effect in the production of ozone from molecular oxygen -the role of molecular symmetry in isotope chemistry[END_REF][START_REF] Bains-Sahota | Mass-Independent Oxygen Isotopic Fractionation in a Microwave Plasma[END_REF][START_REF] Anderson | Laboratory measurements of ozone isotopomers by tunable diode-laser absorption-spectroscopy[END_REF][START_REF] Morton | Oxygen fractionation of ozone isotopes O-48 through O-54[END_REF][START_REF] Morton | Laboratory studies of heavy ozone[END_REF][START_REF] Thiemens | Pressure dependency for heavy isotope enhancement in ozone formation[END_REF][START_REF] Mauersberger | Multi-isotope study of ozone -implications for the heavy ozone anomaly[END_REF][START_REF] Mauersberger | Ozone isotope enrichment: Isotopomer-specific rate coefficients[END_REF][START_REF] Janssen | Kinetic origin of the ozone isotope effect: a critical analysis of enrichments and rate coefficients[END_REF][START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF][START_REF] Mauersberger | Oxygen Isotope Processes and Transfer Reactions[END_REF][START_REF] Feilberg | Probing the unusual isotope effects in ozone formation: Bath gas and pressure dependence of the non-mass-dependent isotope enrichments in ozone[END_REF] .

The origin of this physical phenomenon has been extensively discussed [START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF][START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF][START_REF] Gao | Strange and unconventional isotope effects in ozone formation[END_REF][START_REF] Gao | On the theory of the strange and unconventional isotopic effects in ozone formation[END_REF][START_REF] Miklavc | Rates of formation of ozone isotopomers: a theoretical interpretation[END_REF][START_REF] Babikov | A quantum symmetry preserving semiclassical method[END_REF][START_REF] Baker | Classical and quasi-classical trajectory calculations of isotope exchange and ozone formation proceeding through O+O-2 collision complexes[END_REF][START_REF] Babikov | Formation of ozone: Metastable states and anomalous isotope effect[END_REF][START_REF] Fleurat-Lessard | Isotope dependence of the O+O2 exchange reaction: Experiment and theory[END_REF][START_REF] Marcus | Mass-independent isotope effect in the earliest processed solids in the solar system: A possible chemical mechanism[END_REF][START_REF] Pack | Some symmetry-induced isotope effects in the kinetics of recombination reactions[END_REF][START_REF] Schinke | Dynamical studies of the ozone isotope effect: A status report[END_REF][START_REF] Ivanov | On molecular origin of mass-independent fractionation of oxygen isotopes in the ozone forming recombination reaction[END_REF] but is still puzzling.

The isotopic composition of ozone is usually reported in relative variations expressed in per mil and noted in the d m O units:

d m O = {[( m O/ 16 O)O 3 / ( m O/ 16 O)O 2 ] -1 } 1000 (1)
with m standing for mass 17 or 18. Using these units, the isotopic compositions can be reported in a 3 isotopes diagram as d 17 O = f(d [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O) where quasi-linear correlations are defined by their slope p. For an infinite reservoir of O2, the isotopic composition of the reservoir (i.e.

( m O/ [START_REF] Feilberg | Probing the unusual isotope effects in ozone formation: Bath gas and pressure dependence of the non-mass-dependent isotope enrichments in ozone[END_REF] O)O 2 ) remains constant during the synthesis of ozone. In this condition, d m O stands also for the overall isotopic fractionation factor m a between O2 and O3:

d m O = ( m a -1) 1000 (2)
The variations in the two isotope ratios [START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF] O/ [START_REF] Feilberg | Probing the unusual isotope effects in ozone formation: Bath gas and pressure dependence of the non-mass-dependent isotope enrichments in ozone[END_REF] O and [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O/ [START_REF] Feilberg | Probing the unusual isotope effects in ozone formation: Bath gas and pressure dependence of the non-mass-dependent isotope enrichments in ozone[END_REF] O depend only on the difference in mass of the isotopes and, in this respect, are assumed to be mass dependent. In other terms, a change of 1‰ in d [START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF] O is accompanied by a change of ≈2‰ in d [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O; hence p=0.52 [START_REF] Urey | The thermodynamic properties of isotopic substances[END_REF][START_REF] Bigeleisen | The relative reaction velocities of isotopic molecules[END_REF][START_REF] Young | Kinetic and equilibrium mass-dependent isotope fractionation laws in nature and their geochemical and cosmochemical significance[END_REF] .

Contrary to this theoretical and experimental prediction, the MIF effect observed during the synthesis of ozone yields p=1.0.

If the isotope fractionation results from two processes having a mass dependent and a mass independent isotopic fractionation factors, m a,Mass Dep and h respectively, the overall isotopic fractionation factor m a can be expressed as:

m a = m a,Mass Dep . h (3)
with m a, Mass Dep. = 1+e(a)

and

h = 1+ e(h) (5)
At this stage, in order to simplify the notations, we designate [START_REF] Feilberg | Probing the unusual isotope effects in ozone formation: Bath gas and pressure dependence of the non-mass-dependent isotope enrichments in ozone[END_REF] O, [START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF] O and [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O by 6, 7 and 8, respectively. Using this formalism, the variations of p in the 3-isotope diagram can then be understood as follows: when ozone is formed by a mass dependent process, ( |e(h)| << |e(a)| ) the isotopic compositions of O3 are related by p≈0.52 and inversely, when ozone is formed purely by a mass independent process ( |e(h)| >> |e(a)| ) p≈1.00.

The question raised by ozone can then be re-formulated such as: what is the physical origin of e(h) [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] ? The solution of this question is complicated by the fact that the formation of ozone is a 3-body reaction and thus involves several reactions that yield isotopic selections. A simplified system describing the main reactions responsible for the formation of ozone is written as:

O + O2 à O3* (kf*) (6) O3* à O + O2 (kD) (7) O3* + M à O3 + M (kM) (8)
Reaction ( 6) describes the formation of the activated complex O3* (also referred to as metastable ozone) with the rate constant kf * (f for formation). Reaction (7) represents the spontaneous dissociation of the complex where the inverse of the rate constant kD (D for dissociation) stands for the lifetime of the complex: t = 1/ kD. Reaction (8) corresponds to the possible stabilization of O3* by a third body M absorbing the proper amount of internal energy (the density of M is noted [M]). M can be a molecule of the gas phase or a solid surface.

At low pressure (P<1 Torr), kD >> [M].kM ), the formation rate of ozone is written as:

d[O3]/dt = [O].[O2] - * / . -0 [1] -2 = [O].[O2] (𝑘 * 5 . 𝑘 1 [𝑀]. 𝜏) (9)
The isotopic fractionation factor m a is the ratio of the rate constants k*f, M.kM and t :

m a = m Rf x m RM x m R t (10)
with m Rf, m RM and m R t standing for the k*f, M.kM and t ratios involved in the formation of m66 and 666. Although no scientific consensus exists in the literature on the origin of the MIF effect [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF][START_REF] Gao | Strange and unconventional isotope effects in ozone formation[END_REF][START_REF] Gao | On the theory of the strange and unconventional isotopic effects in ozone formation[END_REF][START_REF] Miklavc | Rates of formation of ozone isotopomers: a theoretical interpretation[END_REF][START_REF] Babikov | A quantum symmetry preserving semiclassical method[END_REF][START_REF] Baker | Classical and quasi-classical trajectory calculations of isotope exchange and ozone formation proceeding through O+O-2 collision complexes[END_REF][START_REF] Pack | Some symmetry-induced isotope effects in the kinetics of recombination reactions[END_REF][START_REF] Ivanov | On molecular origin of mass-independent fractionation of oxygen isotopes in the ozone forming recombination reaction[END_REF] , it has been proposed [START_REF] Reinhardt | Mass independent isotope fractionation in ozone[END_REF][START_REF] Reinhardt | On the mass independent isotopic fractionation in ozone[END_REF] that m Rf and m RM are mass dependent contrary to the lifetime ratios m Rt of the different isotopomers of the complex: the lifetime t would be different if the complex is formed by distinguishable (such as m+66) or by indistinguishable isotopes (such as 6+66) [START_REF] Reinhardt | Mass independent isotope fractionation in ozone[END_REF][START_REF] Reinhardt | On the mass independent isotopic fractionation in ozone[END_REF] . The channels accounting for isotope exchange and non-exchange between O and O2 can be split when isotopes are distinguishable; but, this is not anymore possible when they are undistinguishable. As a consequence, the relative proportions of complexes formed by the exchange and non-exchange processes, cannot be calculated similarly when complexes formed via dis-or indistinguishable isotopes. Their average lifetime t is thus different. This effect takes place during the scattering of oxygen atoms by O2 and is independent of the interaction potential between O and O2. In other terms, MIF does not result from an isotope-exchange reaction, but reflects the isotopic composition of the fraction of the activated complexes that is sampled through their stabilization as ozone (reaction 8). The lifetime ratio m Rt is ascribed to the mass independent isotopic fractionation factor and is noted m h [START_REF] Gao | Strange and unconventional isotope effects in ozone formation[END_REF] .

m h has been modeled numerically in classical mechanics [START_REF] Reinhardt | Mass independent isotope fractionation in ozone[END_REF] , assuming that the quantum terms are negligible when rates are averaged in a thermal situation. This assumption has been numerically verified by [START_REF] Rao | Huge Quantum Symmetry Effect in the O + O2 Exchange Reaction[END_REF] . Several numerical predictions arising from these calculations are specific to this theory and can be used for practical purposes. In this paper, we report an experiment dedicated to test one of these predictions, namely: m h should vary with the lifetime of the complex hence with the pressure of the bath.

In the laboratory, ozone (O3) produced in a molecular oxygen plasma by electric or microwave discharge or by O2 photolysis, is condensed on the glass walls of the apparatus at liquid nitrogen temperature. At O2 pressures <10 Torr, a progressive transition has been observed [START_REF] Bains-Sahota | Mass-Independent Oxygen Isotopic Fractionation in a Microwave Plasma[END_REF][START_REF] Morton | Laboratory studies of heavy ozone[END_REF][START_REF] Thiemens | Pressure dependency for heavy isotope enhancement in ozone formation[END_REF][START_REF] Feilberg | Probing the unusual isotope effects in ozone formation: Bath gas and pressure dependence of the non-mass-dependent isotope enrichments in ozone[END_REF] between the two regimes of isotopic fractionation i.e. with p varying between 1.0 and 0.52 during the decrease in pressure. In addition, when p=1.0 or p=0.52 ozone is either enriched or depleted in the heavy isotopes, respectively. This effect was interpreted as the progressive transition between two distinct processes: a homogeneous (mass independent; p=1.0) and a heterogeneous (mass dependent; p=0.5) reaction taking place on the walls of the apparatus [START_REF] Bains-Sahota | Mass-Independent Oxygen Isotopic Fractionation in a Microwave Plasma[END_REF][START_REF] Morton | Laboratory studies of heavy ozone[END_REF][START_REF] Thiemens | Pressure dependency for heavy isotope enhancement in ozone formation[END_REF][START_REF] Feilberg | Probing the unusual isotope effects in ozone formation: Bath gas and pressure dependence of the non-mass-dependent isotope enrichments in ozone[END_REF] .

Bains-Sahota and Thiemens [START_REF] Bains-Sahota | Mass-Independent Oxygen Isotopic Fractionation in a Microwave Plasma[END_REF] ascribed this low pressure reaction to the formation of O2 by the recombination of two oxygen atoms while Janssen and Tuzson [START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF] called upon the fast isotopic equilibrium between O and O2: m + 66 « 6 + m6 (11) Ozone being formed by the reaction between oxygen atoms and O2 from the bath, the isotopic composition of ozone would reflect the isotopic composition of the oxygen atoms in (11). Since the equilibrium reaction (11) is much faster that the stabilization rate of the complex, Janssen and Tuzson [START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF] have proposed that the reverse reaction (7) overcomes the stabilization reaction (8). However, observations [START_REF] Mauersberger | Oxygen Isotope Processes and Transfer Reactions[END_REF] and theoretical considerations clearly indicate that the complex states reached isotopic equilibrium with the reactants in the gas phase at pressures well above 10 Torr, i.e. well above the domain in pressure where the progressive change of p is observed. Therefore, the theoretical origin of the p variations remains an open question.

We have studied the evolution of p for a finite reservoir of O2 where the isotopic composition of O2 changes along with the removal of O3 and with the change in pressure.

Based on these new results and using the theoretical prediction according to which m h is pressure dependent, we show that another interpretation is possible to account for the change of p at low pressure.

Experimental methods

Ozone is produced by a high-frequency discharge in natural oxygen (research grade purity)

and is continuously removed from the bath through its condensation on the glass walls of the apparatus at liquid N2 temperature. The oxygen gas is introduced in the reactor (about 3700 cm [START_REF] Thiemens | Introduction to Chemistry and Applications in Nature of Mass Indepedendent Isotope Effects Special Feature[END_REF] , cylinder-like shape with 5cm radius and 40cm length) through a Pyrex tube (Figure 1).

The isotopic composition of O2 is recorded as pressure steps i.e. when the reaction stops, the pressure of the remaining O2 is measured and a small aliquot of O2 is sampled and isotopically analyzed. The evolution of pressure is monitored by a piezo-sensitive gauge Pfeiffer APR26. is applied manually with a Tesla coil on one single metallic electrode (inox-made) located at the center of the Pyrex cylinder tube. A glow discharge is visible at the lowermost part of the electrode, where small light sparks can be seen. The lower part of the reactor, about 60% in volume, is immerged in liquid nitrogen (∼87 K).

The reaction and the mass spectrometric analyses are made through the following successive operations (cf. Figure 1): (i) the reactor is outgassed by opening the stopcocks (1) and (3) and similarly, for the admission of O2 in the reactor, (ii) the stopcock (3) is closed and the upper volume (20 cm 3 ) limited by the 3 stopcocks is continuously pumped out, (iii) when the discharge is stopped, the stopcock ( 1) is closed and an aliquot of O2 is sampled by opening the stopcock (3)

(ozone remaining trapped cryogenically on the reactor walls), (iv) after 30 sec. for pressure equilibration, the stopcock ( 3) is closed and an aliquot of O2 is injected in a Dual-Inlet Isotopic Ratio Mass-Spectrometer (DI-IRMS; see Figure 1) by opening the stopcock (2). This last operation yields a decrease of the pressure in the reaction vessel of the order of a few %.

The isotopic composition is measured with a mass spectrometer ThermoFischer Delta V DI-IR-MS. The analytical precision (1 std. dev.) on isotopic analyses is ±0.082 and ±0.040 ‰, for d 17 O and d 18 O, respectively. The mass resolution of the instrument is ≈150 (Dm/m). The amount of analyzed gas lies between 10 and 80 µmole STP.

Results

Results are reported in Table 1 and 5 experiments out of 8 are reported in Figure 2. O2 was introduced in the reaction vessel at an initial pressure Pi. In the Janssen and Tuzson experiment [START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF] the amount of synthetized O3 did not cause a noticeable change in the O2 pressure and thus, the isotopic fractionation factor m a was directly accessible by measuring d 17 O and d 18 O in O3. In contrast to this, the present experiment lead to a reservoir effect (hereafter referred to as a "distillation") where the pressure and the isotopic compositions of O2 diminish along with of the condensation of O3. The duration of the reaction, the initial and final pressures Pi and Pf are also reported in Table 1 for each step.

The evolution of the d 18 O with pressure is reported in the Figure 2 Note that p=0.71 is usually observed when the discharge takes place between two electrodes [START_REF] Bains-Sahota | Mass-Independent Oxygen Isotopic Fractionation in a Microwave Plasma[END_REF][START_REF] Thiemens | Pressure dependency for heavy isotope enhancement in ozone formation[END_REF][START_REF] Janssen | Kinetic origin of the ozone isotope effect: a critical analysis of enrichments and rate coefficients[END_REF][START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF] . We have also reproduced experimentally p=0.71 by using an electric discharge between two electrodes [START_REF] Baraut-Guinet | Unconventionnal mass-independent oxygen isotope effect in ozone by microwave discharge plasma[END_REF] . This departure of p from unity would reflect a mass dependent isotope fractionation occurring either during the dissociation of O3 by electron impact in the discharge or by a simple thermal dissociation of O3. In these conditions, ozone results from a two steps process: a formation mechanism along the slope 1.00 and a dissociation mechanism along the slope 0.52; hence a resulting slope of 0.71. In the course of this study, it has been verified that the slope p is exactly equal to 1.0 for a O2 pressure around 80 Torr and when using a single electrode to promote the discharge. Therefore, this dissociation mechanism is not taking place in the present experiment.

As shown in this Figure ,d 17 O and d 18 O in the remaining O2 becomes progressively depleted in [START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF][START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O when the pressure diminishes indicating that ozone is enriched in [START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF][START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O and lies on a slope p=1.0. The shift of slopes at low pressures from p =1.0 to p »0.52 is associated with an increase in d 17 O and d 18 O in the remaining O2 (p » 0.5; Sample #1), indicating that ozone is now depleted in [START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF][START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O in a mass dependent manner. Such a behavior is usually referred to as a cross-over in the isotopic fractionation factor.

This cross over is qualitatively in agreement with previous observations [START_REF] Bains-Sahota | Mass-Independent Oxygen Isotopic Fractionation in a Microwave Plasma[END_REF][START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF] . It was interpreted [START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF] as an evolution in the relative contributions of (i) ozone formed in the gas phase reaction with p≈0.71 corresponding to the maximum d As a summary, these new results are difficult to reconcile with a two end-member mixing model. They are compared in the next session with the predicted theoretical relation between m a and the lifetime of the complex

Interpretation

The present interpretation is based on the formalism developed by Reinhardt et Robert [START_REF] Reinhardt | Mass independent isotope fractionation in ozone[END_REF][START_REF] Reinhardt | On the mass independent isotopic fractionation in ozone[END_REF] .

Our purpose here is not to perform a complete physical treatment of oxygen chemistry in the plasma. Numerical simplifications can be made in the calculation of the isotopic fractionation m a if one considers only the ratio of isotopically substituted reaction rate constants [START_REF] Young | Kinetic and equilibrium mass-dependent isotope fractionation laws in nature and their geochemical and cosmochemical significance[END_REF] . Indeed, in most cases, the interaction potentials are not affected by the isotopic substitution and such ratios (as m Rf or m RM in (10)) can be supposed to be mass dependent. This assumption is still theoretically investigated in great details [START_REF] Teplukin | Several Levels of Theory for Description of Isotope Effects in Ozone: Symmetry Effect and Mass Effect[END_REF]38) .

In [START_REF] Reinhardt | On the mass independent isotopic fractionation in ozone[END_REF] , a simple model for the stabilization of the complex has been assumed: ozone would be stabilized above the minimum lifetime value tMin. Below tMin, the complex spontaneously dissociates and O2 returns to the gas phase via reaction (7). The relation between the weighted average lifetime 𝜏 and tMin has been computed and is reported in Figure 4 (see [START_REF] Reinhardt | On the mass independent isotopic fractionation in ozone[END_REF] for details).

In this Figure, the difference between 17 𝜏 and 18 𝜏 being negligible, m has been omitted in 𝜏 for simplicity. In the rest of the paper we will use 𝜏 in place of tMin, 𝜏 being an intensive parameter of the reaction forming ozone.

Considering only the natural isotopic abundances of oxygen (i.e. 8/6≈2x10 -3 and 7/6≈4x10 - 4 ), the 3 isotopomers 666 and m66 are the most abundant species. In italic, the knock-on atom of the molecule O2. Reaction (12) stands for the two reactions between m and 6. Note that reactions ( 12) and ( 14) correspond to the forward and backward reaction (11) for isotope exchange. The atom-molecule reactions producing [START_REF] Feilberg | Probing the unusual isotope effects in ozone formation: Bath gas and pressure dependence of the non-mass-dependent isotope enrichments in ozone[END_REF] O [START_REF] Feilberg | Probing the unusual isotope effects in ozone formation: Bath gas and pressure dependence of the non-mass-dependent isotope enrichments in ozone[END_REF] O [START_REF] Feilberg | Probing the unusual isotope effects in ozone formation: Bath gas and pressure dependence of the non-mass-dependent isotope enrichments in ozone[END_REF] O* is:

6 + 66 → 666* ( 15
)
The lifetime for the reaction ( 15) is noted 𝜏I with the subscript I for indistinguishable.

We define the lifetime ratio as the mass independent isotopic fractionation factor h:

h( 𝜏) = 𝜏 / 𝜏I (16)
Experimental data and numerical simulations [START_REF] Janssen | Kinetic origin of the ozone isotope effect: a critical analysis of enrichments and rate coefficients[END_REF][START_REF] Schinke | Isotope dependence of the lifetime of ozone complexes formed in O+O2 collisions[END_REF][START_REF] Schinke | The effect of zero-point energy differences on the isotope dependence of the formation of ozone : a classical trajectory study[END_REF] reveal that both mass dependent and mass independent isotopic mass fractionation takes place in the stabilization rate of O3* at the low pressure regime. For example, the average lifetime of the 666 complex is shorter than the lifetime of 688 but longer than the lifetime of 866 [START_REF] Schinke | Isotope dependence of the lifetime of ozone complexes formed in O+O2 collisions[END_REF] . This mass dependent effect is theoretically well understood and involves the quantum mechanical differences in the vibrational zero-point energies (DZPE) between entrance and exit channels in the O+O2 reaction [START_REF] Janssen | Kinetic origin of the ozone isotope effect: a critical analysis of enrichments and rate coefficients[END_REF][START_REF] Schinke | Isotope dependence of the lifetime of ozone complexes formed in O+O2 collisions[END_REF][START_REF] Schinke | The effect of zero-point energy differences on the isotope dependence of the formation of ozone : a classical trajectory study[END_REF] .

This vibrational contribution (noted n) is related to the difference in the DZPE between the reactants (O+O2) and the products (O3) and has been experimentally measured by [START_REF] Janssen | Kinetic origin of the ozone isotope effect: a critical analysis of enrichments and rate coefficients[END_REF] for several isotopomers of ozone. In the case of the natural abundances of [START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF] O and [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O, each reaction ( 12), ( 13) and ( 14) contributes for 1 3 @ in the overall isotopic fractionation factor [START_REF] Reinhardt | Mass independent isotope fractionation in ozone[END_REF] . This gives the mass independent isotopic fractionation factor m h( 𝜏) (i.e. the overall lifetime ratio) [START_REF] Reinhardt | On the mass independent isotopic fractionation in ozone[END_REF] :

m h( 𝜏) = 1 3 @ [h(𝜏). m n12] + 1 3 @ [h(𝜏). m n13] + 1 3 @ [h(𝜏). m n14] (17)
with the subscripts standing for the reactions ( 12), ( 13) and (14). The computed numerical n values are [START_REF] Reinhardt | On the mass independent isotopic fractionation in ozone[END_REF] :

17 n12 = 0.867; 17 n13 = 1.000 (DZPE=0); 17 n14 = 1.133 (18) 18 n12 = 0.746; 18 n13 = 1.000 (DZPE=0); 18 n14 = 1.254 (19) Note that the masses of the isotopes for the translation and rotation movements were also considered in the calculations of h( 𝜏) [START_REF] Reinhardt | Mass independent isotope fractionation in ozone[END_REF] ; these mass corrections are negligible and are not mentioned here for simplicity. Numerical results for 17 h( 𝜏) and 18 h( 𝜏) as a function of 𝜏 are reported in Figure 5. They exhibit a maximum around 𝜏 ≈ 30 ps which is in agreement with the maximum isotopic fractionation experimentally observed in ozone and in agreement with experimental determinations for 𝜏 [START_REF] Kaufman | M effect in the gas-phase recombination of O with O2[END_REF] . This gives at 𝜏 ≈ 30 ps, 17 h( 𝜏) = 1.114 and 18 h( 𝜏) = 1.116, giving p=1 in the three isotopes diagram.

In the gas phase, 𝜏 is constant and leads to the third order reaction (9). As a consequence, the ratio d 17 O/d [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O exhibits a plateau at pressure P≲100 Torr. However, if ozone is produced on the walls of the apparatus, the reaction should not depend anymore on the density of the third body [M] and the reaction ( 9) becomes second order. This situation is achieved if one supposes that 𝜏 varies as P -1 on the walls.

Looking now at Figure 5, it can be noted that h( 𝜏) varies rapidly with 𝜏. As shown in Figure 5, h(𝜏) → ≈ 1 when 𝜏 → ∞. 𝜏 can be related inversely to the pressure, as the occurrence of the two possible reactions (O3* with the wall at an average time constant and O3* in the gas stabilized by the encounter with M) tend to favor the reaction with the walls when the pressure diminishes. The average lifetime before a stabilizing collision will thus be shorter at high pressure than at low pressure (and, of course shorter than for a case without walls). Therefore, we assume that 𝜏 is inversely proportional to pressure.

𝜏 ∝ P -1 (20)
This hypothesis has two consequences: when the pressure diminishes (i) the reactions (9) should becomes progressively second order and (ii) since h(𝜏) → ≈ 1, the mass independent effect should disappear. This is in qualitative agreement with several observations reported in the literature [START_REF] Bains-Sahota | Mass-Independent Oxygen Isotopic Fractionation in a Microwave Plasma[END_REF][START_REF] Morton | Laboratory studies of heavy ozone[END_REF][START_REF] Thiemens | Pressure dependency for heavy isotope enhancement in ozone formation[END_REF][START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF] where the decrease in d 17 O-d [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O was observed below ≲30 Torr and ascribed to wall effects. Ascribing this effect to a change in the lifetime of the complex stabilized as ozone, the observed pathways in the d 17 O-d [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O diagram can be simulated numerically.

We define m f as the isotope exchange equilibrium constant of the reactions ( 12) and ( 14) (m+66 and 6+6m):

m f = ( m O/ 16 O)O / ( m O/ 16 O)O 2 (21)
The term m f was not considered in our previous theoretical approach [START_REF] Reinhardt | On the mass independent isotopic fractionation in ozone[END_REF] that was focused on the numerical evaluation -and origin -of the mass independent fractionation factor m h( 𝜏). In the present experiment where all isotopomers are scrambled, the overall isotopic fractionation factor m a( 𝜏) becomes:

m a( 𝜏) = m f . m h( 𝜏) (22)
Note that m f<1 while m h( 𝜏) >1 i.e. are of opposite signs in d units. Consequently, when the mass independent isotopic effect becomes negligible (i.e. m h(𝜏) ≈1), ozone should be depleted in [START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF][START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O relative to O2 with p≈0.5 while, when m h( 𝜏) > m f, ozone should be enriched in [START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF][START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O relative to O2 with p≈1. In the following calculations, we use the measured [START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF] (and calculated [START_REF] Kaye | Enhancement of Heavy Ozone in the Earth's Atmosphere?[END_REF] ) values for f: 18 f=0.965 (i.e. d 18 O= -35‰) and 17 f=0.982 (i.e. d 17 O= -18‰).

The isotopic fractionation of oxygen isotopes during the distillation of O2 via the condensation of O3 can be written as:

d[ m O] / d[ 16 O] = m a([ 16 O]) . {[ m O]/[ 16 O]} (23)
with m a([ 16 O]) being a function of the overall pressure:

m a([ 16 O]) ∝ m a(P -1 )) ∝ m a( 𝜏)) (24)
When m a is constant, equation ( 23) yields the so-called Rayleigh equation. Using the variations m a( 𝜏) reported in the Figure 5, the equation ( 23) has been solved numerically by adjusting t=30 ps for P =30 Torr. Results for d 17 O and d 18 O are reported in Figure 6. They reproduce the data with no additional free parameter. For the sake of clarity, the present calculation is still a first order approximation since the measured isotopic values stand for the function (23) integrated between the initial and final pressure of each step. We neglect this refinement.

Note that f values vary with the temperature of the bath. This temperature is difficult to estimate since there is a strong temperature gradient between the center of the plasma and the walls at Liq. N2 temperature. However, as long as (1-17 f).1000 = 0.52 (1-18 f ).1000 (i.e. 17 f mass dependently related to 18 f), the numerical simulation is (almost) insensitive to f values for 0.94< 18 f<0.98.

In experimental literature [START_REF] Bains-Sahota | Mass-Independent Oxygen Isotopic Fractionation in a Microwave Plasma[END_REF][START_REF] Morton | Laboratory studies of heavy ozone[END_REF][START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF] , the O2 reservoir is infinite i.e. O2 is continuously renewed in the reaction chamber. In such condition the mixing ratio between ozone formed in the gas phase and on the walls cannot be constrained by a distillation model but by the reaction rate ratio between ozone formed on the walls and formed in the gas. As shown by Morton et al. with ozone formed on the walls (p=0.52) resulting from a different reaction than ozone formed in the gas (p=1.0). This property of the evolution of the ozone d 17 O/d 18 O has been observed in all previously reported experiments [START_REF] Bains-Sahota | Mass-Independent Oxygen Isotopic Fractionation in a Microwave Plasma[END_REF][START_REF] Morton | Laboratory studies of heavy ozone[END_REF][START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF] . The change in h with the lifetime of the complex as reported in Figure 5, should also define a straight line in the d As mentioned previously, a cross over in the d 17 O-d [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O distillation pathway can also be predicted by previous experimental results [START_REF] Bains-Sahota | Mass-Independent Oxygen Isotopic Fractionation in a Microwave Plasma[END_REF][START_REF] Janssen | Isotope Evidence for Ozone Formation on Surfaces[END_REF] . We have therefore performed a similar distillation calculation with constant isotopic fractionation factors (i.e m a( τ) and m f). Contrary to the assumption where m a( τ) changes with the lifetime of the complex, these conditions do not reproduce the pathways defined by the present data.

The present observations confirm that glass walls trigger the stabilization of the complex but the numerical treatment of the distillation does not call upon other formation reactions than those involved in the gas phase (eq. 9). In other terms, the change in sign in the

d 17 O/d 18 O
values is not caused by the evolution of the mixing ratio between two ozone species but by the evolution of the overall isotopic fractionation factor. The stabilization on the walls yields a selection in the lifetime distribution of the complex.

This interpretation raised one issue: the transition between the two regimes of stabilization of the complex [O3]* should theoretically occur around 10 -2 Torr and not around 10 Torr as postulated here. This argument is also reinforced experimentally: in an afterglow of pulsed discharges [START_REF] Marinov | Direct observation of ozone formation on SiO2 surfaces in O2 discharges[END_REF] , it has been shown that the O3 surface production rate depends on the density of adsorbed O and O2 on solid surfaces and consequently, that ozone produced on the walls in the range 1-5 Torr should not contribute significantly to its overall production. This point remains an open issue.

Conclusions

The isotopic variations in d 17 O/d [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O observed experimentally during the condensation of ozone in a close reservoir of O2, are in agreement with variations with pressure of the mass independent isotopic fractionation factor h [START_REF] Reinhardt | On the mass independent isotopic fractionation in ozone[END_REF] . The cross over in the isotopic evolution of the remaining reservoir would not be caused by an additional reaction forming ozone but would be intrinsic to the overall isotopic fractionation factor.

In this model, the "anomalous" isotopic effect is not restricted to the formation of ozone nor to oxygen isotopes. If correct, this opens a wide range of possible applications to experimental or natural systems for other chemical elements.

(3)

Figure 1: Oxygen discharge experiment. After condensation of ozone at liquid nitrogen temperature on the walls of the vessel, an aliquot of the remaining O2 is sampled via the Aliquot set up and send directly to the mass spectrometer (see Experimental). 1). They define distillation pathways caused by the continuous removal of O 3 by its condensation at liquid nitrogen temperature on the walls of the vessel. The evolution of these pathways depends on the initial pressure of O 2 . The mass dependent (p=0.52) and mass independent (p=1) fractionation lines are shown for references. Note that the pathways do not define converging linear correlations.

  . In this figure, the pressure stands for the average value between Pi and Pf except for the first step where the pressure is the initial Pi value. The d 18 O shows first a decrease in the first pressure steps, followed by an increase for the last fractions of the remaining O2 reservoir. During the decrease in pressure, the d 18 O's go through minimum values which occur at different pressures among the experiments, depending on the initial pressure of the O2 reservoir. This minimum value corresponds to a 45-55% decrease in pressure relative to the initial pressure value.The evolutions of the isotopic compositions of O2 caused by the removal of O3 are also reported as pathways in the 3 isotopes diagram (Figure3). In several experiments, d 17 O and d[START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O show first a decrease along a slope p»1.0 followed by an increase along a slope1>p>0.5. The correlation slope p varies between »1.0 and »0.52 i.e. between a pure mass independent and a pure mass dependent fractionation.

17 O and d 18 OFigure 3 ,

 17183 Figure 3, it can be observed that the linear relations of the isotopic evolutions toward positive d 17 O and d 18 O values do not converge toward a unique composition that would stand for the component (ii).

( 9 )

 9 , the evolution of d 17 O/d 18 O with the pressure depends also on the size of the reservoir, demonstrating the role of the walls in the overall production rate of ozone. In agreement with the present results, Morton et al. (9) reported, for a volume of reaction (volume 35 mL, r»5cm) similar to the present experiment, a change in p beginning around 30 Torr. However, in conditions of an infinite reservoir, the evolution of ozone with pressure defines straight lines in the d 17 O/d 18 O diagram compatible with a two end-member mixing model i.e.

  [START_REF] Gellene | An explanation for symmetry-induced isotopic fractionation in ozone[END_REF] O-d[START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF] O diagram. Therefore, for experiments performed with an infinite reservoir of O2, it is impossible to decide if the evolution of the isotopic composition of ozone in the d 18 O-d[START_REF] Krankowsky | Atmospheric chemistry -Heavy ozone -A difficult puzzle to solve[END_REF] O diagram should be ascribed to h or to the evolution of the mixing ratio of ozone produced by two different reactions (in the gas and on the walls).

Figure 2 :

 2 Figure 2: The isotopic composition (d 18 O) of molecular oxygen is reported as a function of Pressure. Note that the minimum d 18 O values is reached at different pressure and is related to the initial pressure.

Figure 3 :

 3 Figure 3: The isotopic compositions of O 2 are reported in the 3 isotopes diagram in d m O units (data in Table1). They define distillation pathways caused by the continuous removal of O 3 by its condensation at liquid nitrogen temperature on the walls of the vessel. The evolution of these pathways depends on the initial pressure of O 2 . The mass dependent (p=0.52) and mass independent (p=1) fractionation lines are shown for references. Note that the pathways do not define converging linear correlations.

t

  Min(ps) t(ps)

Figure 4 :

 4 Figure 4: The relation between the weighted average lifetime ! and the minimum lifetime t Min . The difference between 17 ! and 18 ! being negligible at this scale, has been omitted.

Figure 6 .

 6 Figure 6. The computed isotopic compositions of O 2 in a distillation model are reported in the 3 oxygen isotopes diagram. The mass dependent (MDF) and mass independent fractionation (MIF) lines are shown as references. For comparison with experimental data, the dot symbols are used for calculated values while larger symbols are used for experimental data.

  The mass independent isotopic fraction factors m h (t) = ( m O/ 16 O) O3 / ( m O/ 16 O) O2 with m standing for mass 17 or 18 (cf. equation (17)) are reported as a function of the average lifetime ! of the complex O 3 *. The notations 6, 7 and 8 stand for 16 O, 17 O and 18 O. m h (t) is calculated for a scrambled mixing of the isotopomers 6-6m, 6-m-6 and m-66.

		1.130
		1.110 1.120		8 6 6
	)	1.100
	( t	1.090		7 6 6
	m a	1.080
		1.070
		1.060
		1.050
		1.040
		20.00	25.00	30.00	35.00	40.00	45.00
				t (ps)
	Figure 5:

Table 1 :

 1 In the successive columns, are reported: The samples in their chronological sampling order, the duration of the discharge, the pressure variation of the O2 reservoir and the isotopic compositions of the remaining O2 expressed in d m O (‰) units. (n.m. for not measured).

	Sample	Duration DPressure	d 17 O(‰)	d 18 O(‰)
		(min)	(Torr)
	#1,0	0	5.39	-0.08 ± 0.08	-0.14 ± 0.02
	#1,1	15	5.39, n.m.	2.73 ± 0.11	6.48 ± 0.05
	#1,2	30	3.22, 1.80	10.20 ± 0.21	21.89 ± 0.05
	#2,0	0	10.43	-0.08 ± 0.08	-0.14 ± 0.02
	#2,1	24	10.43, 6.5	-2.32 ± 0.05	2.25 ± 0.04
	#2,2	37	4.13, 2.58	-1.29 ± 0.16	4.43 ± 0.09
	#3,0	0	10.33	-0.08 ± 0.08	-0.14 ± 0.02
	#3,1	22	10.33, 5.83	-4.79 ± 0.07	-0.64 ± 0.04
	#3,2	37	4.9, 2.56	-1.48 ± 0.14	8.17 ± 0.05
	#3,3	48	2.06, 1.65	3.72 ± 0.09	18.41 ± 0.09
	#4,0	0	11.45	-0.08 ± 0.10	-0.12 ± 0.02
	#4,1	24	11.45, 6.30	-5.28 ± 0.08	-2.08 ± 0.03
	#4,2	47	5.05, 3.75	-3.00 ± 0.06	4.68 ± 0.04
	#5,0	0	18.05	-0.03 ± 0.05	-0.04 ± 0.03
	#5,1	22	18.05, 12.52	-7.96 ± 0.05	-6.79 ± 0.01
	#5,2	35	11.26, 9.21	-9.51 ± 0.04	-6.81 ± 0.03
	#5,3	51	8.50, 5.83	-9.83 ± 0.06	-4.32 ± 0.03
	#6,0	0	n.m.	-0.02 ± 0.08	-0.11 ± 0.30
	#6,1	19	26.52, 20.60	-8.51 ± 0.08	-8.27 ± 0.03
	#6,2	44	19.1, 11.99	-16.44 ± 0.03	-13.69 ± 0.03
	#6,3	68	10-79, 7.94	-15.59 ± 0.08	-10.13 ± 0.04
	#6,4	80	6.77, 4.97	-13.88 ± 0.13	-5.76 ± 0.03
	#6,5	97	4.14, 3.26	-7.57 ± 0.07	5.90 ± 0.04
	#7,0	0	n.m.	-0.02 ± 0.09	-0.09 ± 0.03
	#7,1	19	30.16, 23.03	-8.93 ± 0.06	-8.76 ± 0.04
	#7,2	44	21.49, 14.85	-18.80 ± 0.06	-16.77 ± 0.02
	#7,3	70	13.60, 9.01	-24.45 ± 0.04	-19.51 ± 0.03
	#7,4	80	8.18, 7.13	-24.60 ± 0.10	-18.28 ± 0.02
	#7,5	87	6.36, 6.12	-24.41 ± 0.07	-16.37 ± 0.02
	#7,6	95	5.51, 5.20	-25.19 ± 0.07	-15.35 ± 0.01
	#7,7	117	4.36, 2.46	-21.49 ± 0.02	-5.42 ± 0.03
	#8,0	0	37.56	-0.07 ± 0.05	-0.10 ± 0.01
	#8,1	28	37.56, 26.29	-11.85 ± 0.10	-11.68 ± 0.04
	#8,2	57	23.92, 16.72	-21.56 ± 0.06	-19.63 ± 0.04
	#8,3	73	15.03, 11.05	-26.67 ± 0.15	-22.69 ± 0.03

Acknowledgements

This work has been founded by the ERC Advanced Grant PaleoNanoLife (PI: F. Robert; 161764). We wish to thank Mark Thiemens, Dmitry Babikov, Fabien Gatti and Christof Janssen for their time in long debates. We also thank Guillaume Lombardi and Karim Hassouni for their expertise in plasma physics