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

HAL
open science

# Diffusion, chemical bonding, and kinetic fractionation of noble gases in the primordial magma ocean

Adrien Saurety, Ozge Ozgurel, Chris Mohn, Razvan Caracas

► To cite this version:

Adrien Saurety, Ozge Ozgurel, Chris Mohn, Razvan Caracas. Diffusion, chemical bonding, and kinetic fractionation of noble gases in the primordial magma ocean. Geochimica et Cosmochimica Acta, 2024, 378, pp.144-152. 10.1016/j.gca.2024.06.021. hal-04754196

HAL Id: hal-04754196

https://hal.science/hal-04754196v1

Submitted on 27 Oct 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

# Diffusion, chemical bonding, and kinetic fractionation of noble gases in the primordial magma ocean

Adrien Saurety<sup>a,b,*</sup>, Ozge Ozgurel<sup>b</sup>, Chris Mohn<sup>b</sup>, Razvan Caracas<sup>a,b</sup>

## ARTICLE INFO

Associate editor: Manuel Moreira

**Keywords:**

Ab initio molecular dynamics

Magma ocean

Noble gases

## ABSTRACT

Noble gases record in their radiogenic/non-radiogenic isotope ratios the traces of accretion history and degassing events from the early Earth, including giant impacts and the magma ocean stage, to large-scale recent volcanism. The diffusional behavior of noble gases in silicate magmas is one of the essential parameters needed to correctly model devolatilization from the magma ocean and subsequent volcanic eruptions. Here, we study the diffusivity of He, Ne, Ar, Kr, and Xe in molten pyrolite using ab initio molecular dynamics. We find that the diffusivity of the noble gases decreases with increasing atomic size. The diffusion of Ar is comparable to that of Mg. We compute the partial molar volumes and find that at high pressure, He is almost one order of magnitude more compressible than the heavier noble gases. The analysis of the speciation of the noble gases in the melt indicates that the lifetime of the coordination polyhedra formed with O is not sensitive to the mass of the noble gases. Finally, we compute the electron localization function along the Xe–O interatomic vectors and show that there is no chemical bonding between the Xe and surrounding O atoms from the hosting silicate melt.

## 1. Introduction

Noble gases are inert volatiles that record information about long-term mantle evolution and Earth's degassing history in their radiogenic vs. non-radiogenic isotopes ratio. The stable non-radiogenic isotopes are of primordial origin, formed before the solar system, and accreted to the early Earth. Since they are continuously degassed throughout geological time, only reservoirs entirely isolated from the surface may contain significant amounts of primordial isotopes. By contrast, radiogenic isotopes are by-products of common nuclear decay reactions continuously replenished over geological time.

A series of recent first principles simulations (Xiong et al., 2021; Yuan and Steinle-Neumann, 2021; Li et al., 2022; Wang et al., 2022; Ozgurel and Caracas, 2023) have shown that the noble gases are preferentially partitioned into the magma ocean (MO) during the core formation stage of the early Earth. The partition coefficient between the liquid iron and the molten silicate is $10^{-4}$ to $10^{-2}$, depending on the compositions of the two liquids, the temperature, and the pressure. The magma ocean must therefore have stored almost the entire amount of the noble gases available during core formation.

As the cooling of the global MO advanced, crystals started to separate

into a mushy layer, which accumulated roughly at mid-depths in the MO due to neutral buoyancy (Labrosse et al., 2007; Ballmer et al., 2017; Caracas et al., 2019). Once the rheological transition was reached (Abe, 1995; Vigneresse et al., 1996), the global MO separated into a shallow MO (SMO) and a basal MO (BMO). This led to different degassing behaviors of the two developing reservoirs.

Since noble gases are chemically inert, only physical processes such as diffusion and behavior under compression can affect their isotopic and elemental fractionation. As the BMO was separated from the surface, depending on the crystallization dynamics, these elements could have been trapped in the crystallizing bridgmanite and magnesiowüstite or have remained in the magma. In the former case, the noble gases could have been dispersed by solid convection throughout the solid lower and upper mantle. In the latter case, they could have accumulated at the base of the mantle, in isolated reservoirs, like the large-low shear velocity provinces (Trønnes et al., 2019; Burke et al., 2008) or in the ultra-low velocity zones (Trieloff and Kunz, 2005; Class and Goldstein, 2005; Ozgurel and Caracas, 2023).

In the SMO, strong stirring followed by vigorous convection dominated the dynamics of the devolatilization. During this stage, large amounts of gas could have degassed if able to reach the surface and be in

<footnote><sup>a</sup> Institut de Physique du Globe de Paris, Université Paris Cité, CNRS, 1 Rue Jussieu, 75005 Paris, France</footnote>

<footnote><sup>b</sup> The Centre for Earth Evolution and Dynamics (CEED) and The Center for Planetary Habitability (PHAB), University of Oslo, 0371 Oslo, Norway</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: saurety@ipgp.fr (A. Saurety).</footnote>
contact with the atmosphere (Salvador and Samuel, 2023; Ozgurel and Caracas, 2023). One of the limiting factors of these processes is related to the atomic diffusion of the noble gases in the melt.

Hence, atomic diffusion in the ultramafic liquids of the MO is a fundamental parameter, which needs to be considered in magma convection models. If large-scale mixing occurs due to convection cells, local-scale element homogenization and devolatilization of the melts in the magma ocean is determined by diffusion. However, the experimental and theoretical data for diffusion in molten silicates are limited to low pressures (<5GPa) and relatively low temperatures (<2000 K) (Nowak et al., 2004; Lux, 1987; Guillot and Sator, 2012). Instead, with ab initio simulations we can access conditions typical of the entire magma ocean and compute various physical/chemical properties, including the diffusion coefficients.

Hence, here we study the physical properties of noble gases in the magma ocean at conditions consistent with the temperatures and pressures of the magma ocean of the early Earth. For this, we employ ab initio molecular dynamics. We analyze in detail the structural changes of the pyrolite melt induced by the presence of noble gases and calculate their partial molar volumes. Then, we determine the self-diffusion coefficients of all noble gases in the molten silicate. Lastly, we investigate the possibility of chemical bonds occurring between Xe and O, and we find no proof thereof. We conclude that the behavior of noble gases in the magma ocean is ruled by physical aspects, like mobility and compaction.

## 2. Methods

We perform ab initio molecular dynamics (AIMD) simulations using the VASP software package (Kresse and Furthmüller, 1996; Kresse and Furthmüller, 1996; Kresse and Hafner, 1993). The interatomic forces are obtained from density-functional theory (Kohn and Sham, 1965) in the projector augmented wavefunction (PAW) approach. The kinetic energy cut-off of the plane waves outside of the PAW spheres is set to 550 eV; the cut-off for the Hänkel functions, which describe the augmented charge inside the PAW spheres, is set to 800 eV. The exchange-correlation term is approximated by the Perdew-Burke-Erzernhof formalism of the generalized gradient approximation (Perdew et al., 1996). The electronic temperature of the system is treated using the Fermi distribution via the smearing width as (k_B T) (Mermin, 1965). The Brillouin zone is sampled at the Gamma point.

We perform the simulations in the canonical (NVT) ensemble, where the number of particles (N) and the volume (V) are fixed. The temperature (T) is set by an external Nose-Hoover thermostat with a mass of 4. The time steps are 1 fs long.

We use pyrolite (Caracas et al., 2019; Solomatova and Caracas, 2019) as a model to approximate the composition of the bulk silicate Earth (McDonough and Sun, 1995). The dry pyrolite melt is modeled using a simulation box with 153 atoms and NaCa₂Al₃Fe₄Mg₃₀Si₂₄O₈₉ molar composition. This approximates the silicate Earth composition (McDonough and Sun, 1995) to better than 1 wt%. To minimize the interactions between the noble gases, the simulations are performed by inserting only two noble gas atoms into the simulation cell, one noble gas at a time. These yield noble gas concentrations of 0.25 wt%, 1.22 wt%, 2.42 wt%, 4.95 wt%, and 7.52 wt% for the He, Ne, Ar, Kr, and Xe, respectively. These values are significantly larger than those of the actual mantle and fluids, estimated to be on the order of hundreds of ppm (Tucker et al., 2018; Kendrick et al., 2011). However, these concentrations are sufficiently small to avoid noble gas interactions between periodic images. Based on our previous experience (Solomatova and Caracas, 2019), we expect the bulk properties of volatiles in terrestrial magma oceans to be correctly described, even if the concentrations are somewhat different than the real ones.

For each simulation, we allow at least one picosecond for thermalization. This ensures even distribution of the atomic velocities and equilibration of the energy and temperature. The production simulations last for at least six picoseconds, with the higher-pressure

simulations being longer. The number of atoms in the simulation box and the mass of the thermostat induce fluctuations of the pressure on the order of 3 GPa and of temperature on the order of 300 K, with a wavelength of about 30–40 fs. Thus, even in the shortest of the simulations, the studied systems perform at least 200 oscillations of the different thermodynamic parameters around an average value, which confer reliable bulk properties of the melt. The post-processing analysis is done with the UMD package (Caracas et al., 2021).

## 3. Results

### 3.1. The structure of the melt

The first atomic coordination sphere in liquids is revealed by analyzing the pair distribution function (PDF). This provides a measure of the number of coordinating atoms as a function of distance to a central atom. The PDFs are always averaged with respect to all the atoms and over the entire simulations cell and finally normalized asymptotically to 1.0 at large distances.

The PDF starts with an exclusion zone, where the coordination is zero as the interatomic repulsion is strong. At the end of the exclusion zone, the PDF starts to increase, and the first peak yields the most likely bond distance. This is often confounded with the average bond distance. The first minimum of the PDF defines the radius of the first coordination sphere, the second minimum defines the second radius of the second coordination sphere, and so on. Due to the limited size of the simulation cells, we only analyze the first coordination shell in our calculations. The integral of the PDF up to the radius of the first coordination sphere gives the average coordination number for that pair of atoms.

Fig. 1 shows the PDFs for all the noble gases in pyrolite as a function of pressure, along the 5000 K isotherm, while Figure S1 presents a few other examples at 4000 K. Over the 0–110 GPa pressure range, the interatomic distances show reductions of up to 19% for He–O and up to 15% for the rest of the noble gases. The coordination radii over the same pressure range decrease by about 13% for He–O, and by about 17 to 16% for the heavier noble gases. This reduction in the radius of the first coordination sphere shows that the noble gases lie closer and closer to the surrounding oxygen when the pressure increases. Compared to the van der Waals radius of the noble gas atoms in ideal gas conditions, even the slightest compression induces a large decrease in the atomic radius. This is also clearly seen in the partial molar volumes, discussed in detail below.

Further, using the radius of the first coordination spheres, we can construct the connectivity matrix between all the atoms at all simulation steps. From there, it is just a small step to track the behavior of all atoms, reconstruct their trajectories, and analyze the time evolution of the structure. Fig. 2 shows such a time analysis: it depicts the lifetime of all atomic polyhedra having a noble gas as the central atom and oxygen as a ligand. The image is obtained from the simulations at 100 GPa and 5000 K.

The sizes of the coordinating polyhedra increase with increasing atomic size of the central atoms, as expected. The larger the noble gas, the larger its interaction range; thus, a higher number of O lies in its coordination sphere. At these pressures and temperatures, He has between 6 and 10 coordinating oxygen atoms, Ar has between 8 and 14, and Xe has as many as 10 to 16. But comparing the lifetimes of these atomic clusters yields remarkable similarities between the different noble gases. The maximum lifetimes are between 50 and 100 fs, and the average lifetimes are around a dozen femtoseconds. Even large gases, like Xe, do not form XeOx clusters that can last for picoseconds or even for hundreds of femtoseconds. This represents a first marker for the fast diffusion of Xe and the lack of bonding between Xe and O.

### 3.2. Partial molar volumes

The partial molar volume of a chemical species in the melt is
<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Pair distribution function (PDF) between noble gas atoms and oxygen atoms as a function of the distance between the two atoms. Simulations were performed at 5000 K, at 10 GPa and 110 GPa. The first maximum is the most likely interatomic distance, and the first minimum marks the radius of the first coordination sphere. The decrease in the interatomic distance and the radius of the first coordination sphere with increasing pressure is highlighted by the colored arrows.</figcaption></figure>

computed from the change in the equation of state (EOS) of the melt induced by the presence of this species. For this, we first obtain the EOS by fitting the pressure–density points along each of the three isotherms with a third-order Birch–Murnaghan EOS for the pure pyrolite melt. Then, we fit an EOS for each simulation with pyrolite melt and dissolved noble gases at each temperature. The difference in volume between the two EOS (with and without noble gas) at a given pressure and temperature yields the partial molar volume of that noble gas at those conditions. This sequence of partial molar volumes allows us to smoothly interpolate the results over the entire range of magma ocean pressures and temperatures.

Fig. 3a shows the increase of pressure due to the noble gas as a function of the pressure of the pure pyrolite. As expected, the presence of the noble gases in the same volume increases the pressure.

The partial molar volumes of the different noble gases increase with atomic size and decrease with pressure. There is no explicit dependence on temperature. Among all the noble gases, only He has an atomic volume smaller than the average atomic volume of the dry melt (defined as the volume divided by the total number of atoms). The other noble gases (Ne, Ar, Kr, and Xe) have a partial atomic volume equivalent to or superior to the average atomic volume of the pyrolite (Fig. 3b).

We can also express the compressibility of the gases in the molten silicate by comparing their partial molar volumes to the corresponding van der Waals volumes at standard conditions. Fig. 3c shows the compressibility of each noble gas along 5000 K isotherm. Up to about 40

GPa, the He volume relative to its van der Waals volume is reduced by almost 50 %, while the volumes of the other noble gases are reduced by about one-quarter of their initial value. The difference might reside in the electronic structure of the noble gases. As He is small compared to the other atoms, it can be more easily accommodated in the void spaces between the other atoms in the melt. Therefore, He does not interfere with the atomic arrangements of the system when pressure increases. On the other hand, the large size of the heavier noble gases prevents an efficient atomic reorganization under compression. Consequently, the larger compressibility of He compared to the other gases likely reflects the compressibility of the entire system rather than capturing the compressibility of He alone.

### 3.3. Xe–O chemical bonding

While X-ray diffraction studies (Leroy et al., 2018; Sanloup et al., 2005) find large amounts of electrons around and between the Xe and O atoms, the absolute value of the electronic density measured in diffraction is no proof of the existence of a chemical bond. Rather, it is either the localization of the electrons between atoms (which is not visible in X-ray diffraction) or the charge transfer between them that marks a chemical bond. Consequently, we investigate the chemical bonding between the Xe atoms and their coordinating O atoms by analyzing the Electron Localization Function (ELF) (Becke and Edgecombe, 1990; Gibbs et al., 2013; Koumpouras and Larsson, 2020; Silvi and Savin, 1994) and performing the Bader charge analysis (Bader, 1994).

The ELF gives the probability density of finding an electron in the neighborhood of another electron with the same spin. In ELF, ionic bonding is characterized by the localization of the electrons around the anions, and covalent bonding by their localization in the interatomic space.

Fig. 4 shows the ELF projected on three planes intersecting the melt along the Xe–O distance for one snapshot of the melt at 5000 K and 110 GPa. The ELF fails to show any electrons localized along the Xe–O distance in either of the three cross-sections. Instead, they all show the electronic clouds accumulated around each of the Xe and O atoms (Fig. 4c and 4d). There is little electron localization between the atoms, suggesting the absence of covalent bonding.

To quantify the potential charge transfer that might occur during an ionic bonding, we compute the Bader charges (Bader, 1994; Sanville et al., 2007; Tang et al., 2009) of all atoms in the pyrolytic melt, including Xe and the other noble gases. In the Bader analysis, attraction basins of the electronic density are identified around each atom. The saddle points between the basins define the boundaries of the atoms. The electronic density is integrated within the atomic boundaries to calculate the negative atomic charge. The difference with the nuclear charge defines the total atomic charge. Because of the utilization of pseudopotentials, we need to add the core electrons to the total electronic density before performing this analysis. The results of this analysis are presented in Table 1.

As the chemical bonding for real materials is usually a mixture of ionic and covalent bonding, as expected, all the charges on all the atoms that we calculated for pyrolite are fractional. Thus, all cations in the pyrolitic melt (Na, Ca, Mg, Fe, Al, Si) are positively charged, and the oxygen atoms are negatively charged. This is a clear indication of partial ionic bonds between the cations and the neighboring oxygen atoms. Except for iron, the Bader charges are weakly anomalous. On the contrary, iron exhibits positive charges between 0.38 and 1.21. These mark a complex Fe–O bonding in the silicate melt, which might be worth exploring further in a future study.

The charges on noble gases are very weak, either positive or negative, and very close to the nominal neutral values. Even for Xe, which exhibits the largest charges, these have a value of only 0.15. This makes the presence of ionic bonding between oxygen and any of the noble gases, including Xe, highly unlikely.
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Fig. Lifetime of $NgO_x$ cluster at 5000 K and 110 GPa. Ng stands for the noble gases (He, Ne, Ar, Kr or Xe). Each bar of the histogram stands for a cluster, and the height of the bar stands for the lifetime of this cluster. The average lifetime of a chemical compound is written under the compound's name.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** (a) The pressure increases of the silicate magma due to adding two moles of noble gas to one mole of $NaCa_2Al_3Fe_4Mg_{30}Si_{24}O_{89}$ pyrolite. (b) The partial atomic volume of the noble gases in pyrolitic melt as a function of pressure. The dashed line represents the average volume of an atom of the volatile-free pyrolite melt. (c) The compressibility of the noble gases is expressed relative to the van der Waals atomic volume at ambient conditions. The shadowed zones are estimates of the cumulated errors ($1\sigma$) in the equations of state.</figcaption></figure>

Hence, with no electron localization and no charge transfer, we can safely conclude that Xe, and, in general the other noble gases, exhibit no chemical bonding with any of the surrounding atoms, including oxygen.

The Xe–O relation is that of geometrical interaction due to compaction, and not an actual chemical bond, and the formation of the short-lived $XeO_x$ clusters in the melt is entirely due to geometrical effects.
<figure><img src="image_8.png" /><figcaption>Fig. 4. The Electron Localization Function (ELF) in the Xe-bearing pyrolite melt at 120 GPa and 5000 K for the first (a) and second (b) Xe atoms. The three cross-sections in each diagram hold the different Xe–O directions. The thin black lines are the limits of the simulation box. Details of one slice from each of the two diagrams (respectively slices 1.1 and 2.2) are shown in (c) and (d). All diagrams show that there is no localization of electrons between the xenon and the oxygen atoms.</figcaption></figure>

**Table 1**

Computed Bader charges compared to the nominal charges. The range covers all the charges of all the atoms of the same type. Values are computed for one snapshot at 110 GPa and 5000 K. All the cations, as well as the O atoms, show a clear charge transfer. The noble gases have Bader charges close to 0, clearly showing no charge transfer from and to the noble gases and, thus, no participation in ionic bonding.

<table><thead><tr><th>Atomic type</th><th>Bader electronic charge</th><th>Nominal charge</th></tr></thead><tbody><tr><td>Na</td><td>+0.67</td><td>+1</td></tr><tr><td>Ca</td><td>+[1.38.. 1.42]</td><td>+2</td></tr><tr><td>Mg</td><td>+[1.53.. 1.66]</td><td>+2</td></tr><tr><td>Fe</td><td>+[0.38.. 1.21]</td><td>+2</td></tr><tr><td>Al</td><td>+[2.40.. 2.48]</td><td>+3</td></tr><tr><td>Si</td><td>+[2.38.. 3.22]</td><td>+4</td></tr><tr><td>O</td><td>−[0.80.. 1.72]</td><td>−2</td></tr><tr><td>He</td><td>−[0.067.. 0.070]</td><td>0</td></tr><tr><td>Ne</td><td>−[0.016.. 0.075]</td><td>0</td></tr><tr><td>Ar</td><td>−[0.034.. 0.167]</td><td>0</td></tr><tr><td>Kr</td><td>+[0.014.. 0.079]</td><td>0</td></tr><tr><td>Xe</td><td>+[0.062.. 0.332]</td><td>0</td></tr></tbody></table>

### 3.4. Diffusion coefficients

The variation of the mean square displacements (MSD) of the atoms over time yields the diffusion coefficients. Fig. 5 shows several examples of MSD. As is standard for MSD, there is a change of slope: the first part of the trajectories represents the ballistic regime, and the second part is the diffusive regime. The difference between the two regimes is better seen in a log–log plot (Figure S2). In the ballistic regime, the atoms move according to short-term interatomic forces and collisions. The time required to leave the ballistic region may vary depending on the pressure, temperature, and nature of the atom diffusing. For the noble gases, the ballistic regime lasts between 100 and 1000 fs. In the diffusive regime, the collisions are dissipated, and the self-diffusion coefficients can be calculated as the slope of the MSD with respect to time, according

to the Einstein relation:

$$D = \text{MSD}/2nt \quad (1)$$

where D is the diffusion coefficient, t is the time and n is the number of degrees of freedom (3 in our case).

The fit quality increases with the length of the simulation, i.e., the distance traveled, and with the number of atoms. This is particularly obvious when comparing He, which covers the longest distances due to its low mass, to Kr, which travels the shortest distances. For the same reasons, Kr shows more deviations from linearity than He. However, the missing statistics can be compensated with longer simulation times, and Fig. 5 shows that even Kr covers sufficient distances in our simulations to allow us to obtain reliable self-diffusion coefficients.

In general, the diffusion of noble gases decreases with increasing atomic size, which is expected based on atomic size considerations. He diffuses around one order of magnitude faster than Mg, and the diffusion of Ar is comparable to that of Mg.

Due to its small size and the lack of chemical bonding interactions between the He and the ions in the silicate, the light He atoms move by “bouncing” back and forth between larger atoms. The helium trajectory is similar to a (random) mean free path between collisions.

At very high pressure, the diffusion of large noble gases, like Ne, Ar, Kr, and Xe, is similar to that of Mg. It was previously suggested that Ar has similar diffusivities as CO₂ (Nowak et al., 2004; Spickenbom et al., 2010), and indeed we find comparable diffusion coefficients between the two gases (Vuilleumier et al., 2015; Solomatova and Caracas, 2019).

Using the diffusion coefficients calculated at various pressures and temperatures, we can fit an Arrhenius law as:

$$D = D_0 \exp(-(E_0 + P\Delta V)/RT) \quad (2)$$

where D is the calculated diffusion coefficient, D₀ is a pre-exponential factor, E₀ is the activation energy, P is the pressure, ΔV is the activation volume, R is the ideal gas constant, and T is the temperature. The
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Mean square displacements of selected atoms in the volatile-bearing pyrolite melts. The exemplified noble gases are Ar (a) and Kr (b) at 120 GPa and 5000 K, and He (c) at 70 GPa and 4000 K. The slope of the linear regression is proportional to the self-diffusion coefficient D. All the diffusivities computed in this study are shown in panel (d).</figcaption></figure>

equation yields a line in log space in the D-T reference, or a plane in the D-P-T space. We fit our D and T values to Arrhenius equation in log scale to determine E₀ (Fig. 6). The results for the fitted parameters (D₀ and E₀) are listed in Table 2. With these values, we can compute the diffusivity of each noble gas in molten pyrolite in the 0–120 GPa pressure range, and in the 2000 to 5000 K temperature range.

The activation energy, E₀, is a purely macroscopic concept. At the microscopic level, the change in energy associated with the displacement of a given diffusive atom is “drowned” in the energy fluctuations of the entire melt (Fig. 7). It is impossible to distinguish the effect of the jump of the noble gases on the internal energy among all the thermal activities. Therefore, the concept of activation energy at the microscopic scale cannot be associated with one single event like an atomic jump, but with an average likelihood of hopping or moving events.

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Log of the diffusivity plot against 10000/T, with T the Temperature. The lines are fitted using the Arrhenius equation. The slope of the lines yield the E₀ + PΔV terms.</figcaption></figure>

<caption>**Table 2**</caption>

The parameters of the Arrhenius law fitted on the self-diffusion results of noble gases in a pyrolitic melt. The pressure and temperature range where we fit those parameters is from 3000 to 5000 K and 5 to 120 GPa.

<table><thead><tr><th>Noble gas</th><th>ΔV (cm³ mol⁻¹)</th><th>E₀ (kJ mol⁻¹)</th><th>D₀ (×10⁻² cm² s⁻¹)</th></tr></thead><tbody><tr><td>He</td><td>1.17</td><td>85.39</td><td>1.01</td></tr><tr><td>Ne</td><td>0.75</td><td>127.78</td><td>1.13</td></tr><tr><td>Ar</td><td>1.08</td><td>84.85</td><td>0.25</td></tr><tr><td>Kr</td><td>0.67</td><td>118.67</td><td>0.31</td></tr><tr><td>Xe</td><td>1.15</td><td>224.79</td><td>6.67</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 7.** The internal energy as a function of the simulation time. The red square highlights the moment when an Ar atom jumps from one atomic cage to another. There is no visible energy barrier.</figcaption></figure>
### 3.5. Kinetic devolatilization of noble gases

Degassing from diffusion close to the surface dominates the time that the magma ocean has an exposed large-scale surface. Even the formation of bubbles implies a previous stage where atoms diffuse, along the chemical gradients, towards the walls of the bubbles, which help them escape from the melt.

With the computed diffusion coefficients (Fig. 5d), we can estimate the distance traveled by atoms in each time interval as a function of pressure and temperature. To smooth out our results obtained at discrete conditions over the entire magma ocean, we employ the Arrhenius law fit on Table 2 to estimate the diffusion rates of the noble gases in the magma ocean. We consider the moment of its formation and after the formation of the first lid on the surface (about 4000 K according to (Miyazaki and Korenaga, 2019) and we plot the diffusivity of noble gases in the early magma ocean (Fig. 8).

The diffusivity of He is around 10 to 20 times higher than that of Xe. Using the diffusion coefficients and realistic thermodynamic conditions within the magma ocean, we estimate that He atoms near the surface need an average of 20 min to travel 1 cm, which can be seen as a characteristic distance between nucleating bubbles in a mafic magma at about 40 kms depth In contrast, under identical conditions, Xe atoms need an average of 3 h to cover the same distance. Consequently, of all the atoms present in the studied systems, He travels the longest distances due to its low mass, small size, and chemical inertness. Ne moves fast at high temperatures; it is five to ten times faster than Xe in the SMO. Ar and Kr have a similar diffusivity. These differences in diffusivity will bear their imprint on fractional degassing only if the convection time-scale, i.e., the residence time of the magma close to the surface, is shorter than these intervals. Above 100 GPa Ar, Kr, and Xe travel comparable distances to Mg at high pressure.

## 4. Discussion

Close to the surface, the formation of bubbles of volatile compounds such as CO2 or H2O will help noble gases to cross hundreds of meters near the surface. This phenomenon plays a major role in the noble gas degassing of the SMO. Indeed, the noble gases can join the gas phase of these bubbles and travel with them much more efficiently than diffusing in the liquid phase (Mukhopadhyay and Parai, 2019).

However, despite the importance of the vesiculation, at the centi-metric scale (Jambon et al., 1986), the mobility of the atoms, i.e. the diffusivity, should also be considered alongside the concentration of gases to estimate the actual outfluxes of gases from the melt. There is a competition between the solubility and the diffusivity of the noble gases to escape the silicate melt (Mukhopadhyay and Parai, 2019); these yield the actual flux of volatiles that determines the degassing via the Fick laws of diffusion. Solubility data from (Shibata et al., 1998; Iacono-Marziano et al., 2010) suggest that Ne is 5 times more soluble in

silicate melt compared to Ar and Kr, and 20–25 times more soluble than Xe. While Ar and Kr solubility are of the same order of magnitude, the diffusion of Ne at the MO surface conditions is found to be five to ten times faster than that of Xe. Our diffusivity results suggest that, for every Xe atom that reaches the surface and degasses, almost one order of magnitude more Ne atoms undergo the same process. Although the differences in solubilities will tend to extract more Xe than Ne, the slow diffusion of the heavy noble gases tends to counterbalance this fractionation. This difference in the diffusivity could significantly affect the expected proportion of noble gases such as Ne and Xe.

This effect could even be stronger for He and Ne fractionation. Indeed, the solubility of He in silicate at low pressure might up to 3 times higher than that of Ne (Lux, 1987), while He seems to diffuse twice as fast as Ne under the conditions of shallow magma ocean. This means that, under identical conditions, the outflux of He would be six times greater than that of Ne due to the combined effects of higher diffusivity and greater solubility. Though closer to the surface, as the concentration of bubbles increases and the distance separating them decreases, the diffusion has a smaller effect; at these conditions, the outflux is dominated by the solubility.

Nevertheless, our estimations of outflux are closer to present-day MORB values than the estimated values of outflux from the primitive mantle reservoir (Tucker and Mukhopadhyay, 2014; Iacono-Marziano et al., 2010). If these numbers are correct, then this model might help explaining the actual ³He/²²Ne ratio of MORB. Considering the diffusivity of He, Ne, and Ar that we calculate here, we estimate that for residence times of the magma close to the surface (less than 2 m from the surface) of less than 4.6 h, the He/Ne fractionation will be affected at least partly by the difference in diffusivity. Past this time, the He/Ne ratio of the magma will be due to the difference in solubility (Jambon et al., 1986) because of efficient degassing. These time estimates become 5.5 h for the He/Ar ratio.

Other events might alter these gas ratios and need to be taken into account when rebuilding the big picture of the noble gases fractionation. They could include, for example, atmosphere-removing impacts that allow new degassing events, amplifying the fractionation between He and Ne (Tucker and Mukhopadhyay, 2014), different entrapment in hidden reservoirs (Mukhopadhyay and Parai, 2019), or fractionation during the crystallization of the magma ocean.

But the mixing in the magma ocean occurs at different scales. Convection displaces large volumes of melt and ensures stirring over a large scale; it depends on viscosity, density, and temperature difference, amongst other parameters (Salvador and Samuel, 2023). We can safely assume that adding the different noble gases does not change the properties of the pyrolite melt dramatically because of their chemical inertness and low concentrations in the melt. Those convective movements allow noble gases to travel efficiently over hundreds of kilometres in the magma ocean.

As the rheology of the deep magma ocean was marked by the

<figure><img src="image_12.png" /><figcaption>**Fig. 8.** The self-diffusion coefficient of each noble gas in the magma ocean as a function of pressure. Estimated diffusivity of the noble gases immediately after the formation of the magma ocean (a) and after the formation of a surface lid (b). The diffusivity is computed using our computed values, and the magma ocean adiabat comes from Miyazaki and Korenaga (Miyazaki and Korenaga, 2019), which assume the lid formation when the surface reaches 4000 K.</figcaption></figure>
rheological transition occurring at mid-depths (Abe, 1995), this led to the separation of two magma oceans, a shallow one and a basal one (Labrosse et al., 2007; Ballmer et al., 2017; Caracas et al., 2019). As a result, the global magma convection must have been separated into smaller cells, and stirring and turbulence must have been considerably reduced. In this regime, the crystal growth rate was estimated at $10^{-5}$, $10^{-4}$ cm s$^{-1}$ (Solomatov and Stevenson, 1993). For comparison, diffusion rates of the noble gases imply travel velocities on the order of $10^{-4}$. $10^{-5}$ cm s$^{-1}$ close to the surface of the magma ocean and on the order of $10^{-6}$–$10^{-5}$ cm s$^{-1}$ at the bottom of the magma ocean. Therefore, chemical incompatibilities suggest that the noble gases should remain in the melt, while the diffusion mechanics might be slower than the advancement of the crystallization front. In the shallow magma ocean, the gas escape would be highly efficient for all the noble gases. This would also lead to their continuous depletion, as large-scale volcanic eruptions of the early Earth would contribute to massive degassing of this liquid reservoir (Saal et al., 2002). In the basal magma ocean, there is competition between the advancing crystallizing front and the diffusion of the gases. With high diffusion rates, He remains in the melt, and would keep accumulating in the melt until the last droplets of the magma ocean solidify; thus, it would most probably concentrate in the last crystallizing parts of the mantle, like the D'' layer and the ultra-low velocity zones (Coltice et al., 2011; Gleeson et al., 2021; Deschamps et al., 2012). The heavier noble gases, like Xe, could be partially trapped in the freshly forming bridgmanite crystals; Ne would be on the edge of crystal entrapment. Consequently, Ar, Kr, and Xe, trapped in bridgmanite, would be more evenly spread throughout the solid mantle, while He would accumulate in specific deep reservoirs. While the actual distribution of the noble gases in the mantle today is the result of the long evolution of the mantle, involving convection and global circulation (Tackley, 2000), our molecular dynamics simulations of the magma ocean suggest a possible effect of this diffusion on the elemental fractionation taking place already in the beginning and during the crystallization of the magma ocean.

This trapping of Xe in the bridgmanite crystal could be a potential mechanism for retaining Xe in the deep mantle and may provide a plausible explanation for the missing Xenon paradox (Pepin, 1991; Tolstikhin and O'Nions, 1994), while the partial entrapment of Ne could contribute towards the observed $ ^{3}He/$$^{22}Ne ratio.

## 5. Conclusions

Ab initio molecular dynamics simulations of the magma ocean at pressures and temperatures conditions encompassing the magma ocean adiabat reveal the chemical inertness and kinetic fractionation trends of the noble gases. Our results show He as the fastest diffusing species in the pyrolite melt. The other noble gases, from Ne to Xe, have diffusion coefficients that are similar to the ones of Mg, and almost one order of magnitude smaller than He. These differences in diffusion are a major source of kinetic fractionation. Thus, we suggest that the devolatilization of He at the surface of the magma ocean is considerably more efficient than the devolatilization rate of the other noble gases. In the deep magma ocean, Ne and Xe can be partially trapped by the crystallizing front, while He easily escapes it. This leads to potentially different accumulation rates of the noble gases in the silicate melts and the subsequent geochemical reservoirs of the lowermost mantle.

We also obtained the partial molar volumes of the noble gases and observed that the compressibility relative to ambient conditions is larger for He than for the other noble gases. Thus, the small size and large compressibility of He suggest that it can more easily fit into the interstitial space of the pyrolytic melt than the larger noble gases. Therefore, the presence of He does not affect the compaction of the pyrolite melt, whereas incorporating Ne, Ar, Kr, and Xe in melt solution leads to less compressible melts. However, the actual concentration of these gases in a terrestrial setting is too small for any of these effects to be observed.

Finally, we explored the chemical interactions at high pressures and

high temperatures. We computed the Bader charges for all the pyrolite atoms and the noble gases and analyzed the ELF along the Xe–O distances. We found no evidence of charge transfer from or to the noble gases, ruling out ionic bonds to the surrounding atoms. The ELF analysis along the Xe–O directions dismisses possible charge localization between the atom, ruling out covalent bonds between Xe and the surrounding O atoms. Consequently, there is no proof of chemical bonding between the noble gases and the pyrolite, including Xe and O. This is an important result, suggesting that the incorporation of Xe and other noble gases in silicate melts, and most probably in silicate solids, must be a purely geometrical construction resulting from compaction and lacking a chemical basis.

## Author contributions

AS performed the simulations, RC designed the study. All the authors participated in the analysis of the simulations and the writing of the manuscript.

## CRediT authorship contribution statement

**Adrien Saurety:** Writing – review & editing, Writing – original draft, Investigation, Formal analysis, Data curation. **Ozge Ozgurel:** Writing – review & editing, Writing – original draft, Formal analysis. **Chris Mohn:** Writing – review & editing, Writing – original draft, Formal analysis. **Razvan Caracas:** Writing – review & editing, Writing – original draft, Formal analysis, Funding, Design the study.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data are available through Mendeley Data at https://doi.org/10.17632/kzjkdbt27m.1.

## Acknowledgments

RC acknowledges support from the European Research Council under EU Horizon 2020 research and innovation program (grant agreement 681818 – IMPACT), from the Research Council of Norway through project HIDDEN 325567, and access to supercomputing facilities via eDARI stl2816 grants, PRACE RA4947 grant, and Uninet2 NN9697K grant. This work was also supported by the Research Council of Norway through its Centres of Excellence funding scheme, projects number 223272 (CEED) and 332523 (PHAB). The simulations were performed on the betzy supercomputer under Uninett Sigma2 grant NN9697K, and on the irene supercomputer of the CCRT supercomputing center under e-DARI stl2816 and PRACE RA9497 grants. AS acknowledges Davide Ceresoli for his Python code for the equation of states, Yoshinori Miyazaki for his data on the magma ocean’s adiabat, and Guillaume Avice for the fruitful discussions about noble gases cosmochemistry.

## Supplementary material

Supplementary material to this article can be found online at https://doi.org/10.1016/j.gca.2024.06.021.

## References

* Abe, Y., 1995. The Earth's central part: its structure and dynamics, chapter basic equations for evolution of partially molten mantle and core. Terra Scientific 215–230.

* Bader, R.F.W., 1994. Atoms in molecules: a quantum theory. The International series of monographs on chemistry. Oxford University Press, Clarendon Press.
Ballmer, M.D., Lourenço, D.L., Hirose, K., Caracas, R., Nomura, R., 2017. Reconciling magma-ocean crystallization models with the present-day structure of the Earth's mantle. *Geochim. Geophys. Geosys.* 18, 2785–2806.

Becke, A.D., Edgecombe, K.E., 1990. A simple measure of electron localization in atomic and molecular systems. *J. Chem. Phys.* 92, 5397–5403.

Burke, K., Steinberger, B., Torsvik, T.H., Smethurst, M.A., 2008. Plume Generation Zones at the margins of Large Low Shear Velocity Provinces on the core–mantle boundary. *Earth Planet. Sci. Lett.* 265, 49–60.

Caracas, R., Hirose, K., Nomura, R., Ballmer, M.D., 2019. Melt-crystal density crossover in a deep magma ocean. *Earth Planet. Sci. Lett.* 516, 202–211.

Caracas, R., Kobsch, A., Solomatova, N.V., Li, Z., Soubiran, F., Hernandez, J.-A., 2021. Analyzing melts and fluids from ab initio molecular dynamics simulations with the UMD package. *J. Visual. Exp.* e61534.

Class, C., Goldstein, S.L., 2005. Evolution of helium isotopes in the Earth's mantle. *Nature* 436, 1107–1112.

Coltice, N., Moreira, M., Hernlund, J., Labrosse, S., 2011. Crystallization of a basal magma ocean recorded by Helium and Neon. *Earth Planet. Sci. Lett.* 308, 193–199.

Deschamps, F., Cobden, L., Tackley, P.J., 2012. The primitive nature of large low shear-wave velocity provinces. *Earth and Planetary Science Letters* 349–350, 198–208.

Gibbs, G., Cox, D., Boisen Jr., M., et al., 2013. The electron localization function: a tool for locating favorable proton docking sites in the silica polymorphs. *Phys. Chem. Minerals* 30, 305–316.

Gleeson, M., Soderman, C., Matthews, S., Cottaar, S., Gibson, S., 2021. Geochemical constraints on the structure of the Earth's deep mantle and the origin of the LLSVPs. *Geochim. Geophys. Geosyst* 22.

Guillot, B., Sator, N., 2012. Noble gases in high-pressure silicate liquids: a computer simulation study. *Geochim. Cosmochim. Acta* 80, 51–69.

Iacono-Marziano, G., Paonita, A., Rizzo, A., Scaillet, B., Gaillard, F., 2010. Noble gas solubilities in silicate melts: new experimental results and a comprehensive model of the effects of liquid composition, temperature and pressure. *Chem. Geol.* 279 (3–4), 145–157.

Jambon, A., Weber, H.W., Braun, O., 1986. Solubility of He, Ne, Ar, Kr, Xe in a basalt melt in the range of 1250–1600°C: geochemical implications. *Geochim. Cosmochim. Acta.* 58, 3519–3525.

Kendrick, M.A., Scambelluri, M., Honda, M., Phillips, D., 2011. High abundances of noble gas and chlorine delivered to the mantle by serpentinite subduction. *Nature Geosci* 4, 807–812.

Kohn, W., Sham, L.J., 1965. Self-consistent equations including exchange and correlation effects. *Phys. Rev.* 140, A1133–A1138.

Koumpouras, K., Larsson, J.A., 2020. Distinguishing between chemical bonding and physical binding using electron localization function (ELF). *J. Phys.: Condens. Matter* 32, 315502.

Kresse, G., Furthmüller, J., 1996. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. *Comput. Mater. Sci.* 6, 15–50.

Kresse, G., Furthmüller, J., 1996. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. *Phys. Rev. B* 54, 11169–11186.

Kresse, G., Hafner, J., 1993. *Ab initio* molecular dynamics for liquid metals. *Phys. Rev. B* 47, 558–561.

Labrosse, S., Hernlund, J., Coltice, N., 2007. A crystallizing dense magma ocean at the base of the Earth's mantle. *Nature* 450, 866–869.

Leroy, C., Sanloup, C., Bureau, H., Schmidt, B.C., Konópková, Z., Raepsaet, C., 2018. Bonding of xenon to oxygen in magmas at depth. *Earth Planet. Sci. Lett.* 484, 103–110.

Li, Y., VoCADLO, L., Ballentine, C., Brodholt, J.P., 2022. Primitive noble gases sampled from Ocean Island basalts cannot be from the Earth's core. *Nat. Commun.* 13, 3770.

Lux, G., 1987. The behavior of noble gases in silicate liquids: Solution, diffusion, bubbles and surface effects, with applications to natural samples. *Geochim. Cosmochim. Acta* 51, 1549–1560.

McDonough, W.F., Sun, S., 1995. The composition of the Earth. *Chem. Geol.* 120.

Mermin, D., 1965. Thermal Properties of the Inhomogeneous Electron Gas. *Phys. Rev.* 137.

Miyazaki, Y., Korenaga, J., 2019. On the timescale of magma ocean solidification and its chemical consequences: 2. Compositional differentiation under crystal accumulation and matrix compaction. *J. Geophys. Res. Solid Earth* 124, 3399–3419.

Mukhopadhyay, S., Parai, R., 2019, Noble Gases: A Record of Earth's Evolution and Mantle Dynamics, Ann. Rev. Earth Planet. Sci. 47:1, 389-419.

Nowak, M., Schreen, D., Spickenbom, K., 2004. Argon and CO2 on the race track in silicate melts: a tool for the development of a CO2 speciation and diffusion model. *Geochim. Cosmochim. Acta* 68, 5127–5138.

Ozgurel, O., Caracas, R., 2023. The magma ocean was a huge helium reservoir in the early earth. *Geochim. Perspect. Lett.* 25, 46–50.

Pepin, R.O., 1991. On the origin and early evolution of terrestrial planet atmospheres and meteoritic volatiles. *Icarus* 92, 2–79.

Perdew, J.P., Burke, K., Ernzerhof, M., 1996. Generalized gradient approximation made simple. *Phys. Rev. Lett.* 77, 3865–3868.

Saal, A.E., Hauri, E.H., Langmuir, C.H., Perfit, M.R., 2002. Vapour undersaturation in primitive mid-ocean-ridge basalt and the volatile content of Earth's upper mantle. *Nature* 419, 451–455.

Salvador, A., Samuel, H., 2023. Convective outgassing efficiency in planetary magma oceans: Insights from computational fluid dynamics. *Icarus* 390, 115265.

Sanloup, C., Schmidt, B.C., Perez, E.M.C., Jambon, A., Gregoryanz, E., Mezouar, M., 2005. Retention of xenon in quartz and earth's missing xenon. *Science* 310, 1174–1177.

Sanville, E., Kenny, S.D., Smith, R., Henkelman, G., 2007 Apr 15. Improved grid-based algorithm for Bader charge allocation. *J Comput Chem.* 28 (5), 899–908.

Shibata, T., Takahashi, E., Matsuda, J., 1998. Solubility of neon, argon, krypton, and xenon in binary and ternary silicate systems: a new view on noble gas solubility. *Geochim. Cosmochim. Acta* 62 (7), 1241–1253.

Silvi, B., Savin, A., 1994. Classification of chemical bonds based on topological analysis of electron localization functions. *Nature* 371, 683–686.

Solomatov, V.S., Stevenson, D.J., 1993. Nonfractional crystallization of a terrestrial magma ocean. *J. Geophys. Res. Planets* 98, 5391–5406.

Solomatova, N.V., Caracas, R., 2019. Pressure-induced coordination changes in a pyrolitic silicate melt from ab initio molecular dynamics simulations. *J. Geophys. Res. Solid Earth* 124, 11232–11250.

Spickenbom, K., Sierralta, M., Nowak, M., 2010. Carbon dioxide and argon diffusion in silicate melts: Insights into the CO2 speciation in magmas. *Geochim. Cosmochim. Acta* 74, 6541–6564.

Tackley, P.J., 2000. Mantle convection and plate tectonics: toward an integrated physical and chemical theory. *Science* 288, 2002–2007.

Tang, W., Sanville, E., Henkelman, G., 2009. A grid-based Bader analysis algorithm without lattice bias. *J. Phys.: Condens. Matter* 21, 084204.

Tolstikhin, I.N., O’Nions, R.K., 1994. The Earth’s missing xenon: a combination of early degassing and of rare gas loss from the atmosphere. *Chem. Geol.* 115, 1–6.

Trieloff, M., Kunz, J., 2005. Isotope systematics of noble gases in the Earth's mantle: possible sources of primordial isotopes and implications for mantle structure. *Phys. Earth Planet. In.* 148, 13–38.

Trønnnes, R.G., Baron, M.A., Eigenmann, K.R., Guren, M.G., Heyn, B.H., Løken, A., Mohn, C.E., 2019. Core formation, mantle differentiation and core-mantle interaction within Earth and the terrestrial planets. *Tectonophysics* 760, 165–198.

Tucker, J.M., Mukhopadhyay, S., 2014. Evidence for multiple magma ocean outgassing and atmospheric loss episodes from mantle noble gases. *Earth Planet. Sci. Lett.* 393, 254–265.

Tucker, J.M., Mukhopadhyay, S., Gonnermann, H.M., 2018. Reconstructing mantle carbon and noble gas contents from degassed mid-ocean ridge basalts. *Earth Planet. Sci. Lett.* 496, 108–119.

Vigneresse, J.L., Barbey, P., Cuney, M., 1996. Rheological transitions during partial melting and crystallization with application to felsic magma segregation and transfer. *J. Petrol.* 37 (6), 1579–1600.

Vuilleumier, R., Seitsonen, A.P., Sator, N., Guillot, B., 2015. Carbon dioxide in silicate melts at upper mantle conditions: Insights from atomistic simulations. *Chem. Geol.* 418, 77–88.

Wang, K., Xiancai, J., Xiandong, L., Mengzi, Z., Kun, Y., 2022. Partitioning of noble gases (He, Ne, Ar, Kr, Xe) during earth's core segregation: a possible core reservoir for primordial noble gases. *Geochim. Cosmochim. Acta* 321, 329–342.

Xiong, Z., Tsuchiya, T., Van Orman, J.A., 2021. Helium and argon partitioning between liquid iron and silicate melt at high pressure. *Geophys. Res. Lett.* 48 e2020GL090769.

Yuan, L., Steinle-Neumann, G., 2021. The helium elemental and isotopic compositions of the Earth's core based on ab initio simulations. *J. Geophys. Res. Solid Earth* 126 e2021JB023106.