The giant impact and the Earth and Moon formation

Periodic Reporting for period 4 - IMPACT (The giant impact and the Earth and Moon formation)

Reporting period: 2021-03-01 to 2022-02-28

The early age of the solar system has seen numerous giant impacts, on a scale energetic enough to un-make and make proto-planets, planetesimals and moons. The Earth did not escape this pattern; if early giant impacts could have formed ponds or even seas of magma, their signature cannot be directly retrieved today, but only inferred. But the last giant impact was large enough to result in a total breakup of the initial bodies and to transform the proto-Earth and the impactor, Theia, into a disk or a synestia. This object was formed of hot liquid droplets and gas. Large parts of it were probably in supercritical state. After its formation this starts to cool and condense.

The central accumulation of dense matter forms the Earth. The timescale for its accumulation is on the order of a few tens of hours. The central body is surrounded for a long time of an atmosphere. Models predict either the formation of a distinct central body and outer atmosphere or of a continuous transitional state with no sharp distinction, like a surface. The first case happens if the liquid-gas

The two scenarios have tremendous implications, for explaining the observed geochemical similarities between the Earth and the Moon. In the first case the disk must have been extremely homogenous from the beginning, for example because the impactor and the proto-Earth would have been very similar in composition. In the second case this similarity constrain is relaxed, as a persistent continuous object would be better mixed.

Choosing the right scenario depends on many things, like the energy of the impact, the masses, the angular momentum of the system. But also it largely depends on the thermodynamic relations of the different mineral components. As a function of the energy of the impact the supercritical state can be reached. And more importantly, depending on the position of the critical point the supercritical state can be maintained for a shorter or longer time. (Figure 1).

However, some of the most important thermodynamic parameters that are necessary in the modelling of such an impact are missing today: the position of the critical points, the boundaries of the critical state for multicomponent systems, the liquid-vapour equilibria, and the equation of state for both liquids and the supercritical fluids.

The estimations of the critical points currently in use have been based on extrapolations of the thermodynamic parameters from ambient conditions up to the high temperatures and low densities specific to silicate critical points. Only very recently more through approaches based on ab initio simulations have started to address simple relevant mineral systems, like silica and forsterite.

In our project we are computing these missing parts of the phase diagrams for the representative systems of the protolunar disk. So far, during these first 30 months we have addressed a wide set of chemical compositions, both pure idealized phases, i.e. iron, feldspars, brucite, and silica, as well as a realistic average composition of the Earth’s mantle, the Bulk Silicate Earth. We have also considered a series of volatiles, like CO, CO2, He, H2O, and Zn, whose behaviour during the Impact and its aftermath we are tracking.

We obtain the liquid spinodal points from which we infer the position of the critical points. We compute the equation of states of the liquids and of the supercritical fluids, we determine the fluid structure, the chemical speciation, and we calculate the transport properties.

The good realization of the project required a considerable methodological development. In the beginning in parallel with producing the raw data by simulations, we spent a large amount of time building the tools necessary for the interpretation of these raw data. Today we have available a powerful library for analysing molecular-dynamics runs. We are currently in the processing of beautification of the code and performing the last tests. The entire package will be released as open source to the entire community during 2019, and a paper will be submitted to an appropriate journal with the code description.
In the following we describe the methodology and then list the status of the simulations and interpretation of the results for the various mineralogical systems we are working on.
Methodology
For single component systems, the schematic phase relations are represented in Figure 1. (For multicomponent system this represents a cross-section through the compositional space). Solids are typically stable at high density and low temperature, liquids at high density and high temperature and gases at low density. All the grey areas in the Figure are two-phase regions. At high-enough temperature there is no transition between gases and liquids. This is the realm of the supercritical state, above the critical point.
In order to find the critical point, the closure point of the gas-liquid two-phase region (dome), we calculate the pressure and the internal energy as a function of density at different isotherms. Upon decompression, the pressure on the melt decreases and the behaviour of the melts follows the Maxwell construction of the liquid - vapor equilibrium. At high density from high pressures down to the saturation vapor pressure, the melts are stable. Below this, they become metastable with respect to the mixture of vapor and liquid. During this interval the pressure continues to decrease and the melt is more and more stretched, such as void spaces may appear inside the liquid at the atomic level. The minimum of the pressure marks the liquid spinodal. At densities smaller than that of the spinodal, the melts become unstable with respect to a mixture of gas and liquid, and bubbles filled with gas start to differentiate even at the scale of our simulation box. If the volume of the simulation box continue to increase the total pressure on the cell starts to increase, as more and more atoms fill in the gas bubbles.
We compute the variation of the pressure (as well as energy) as a function of density using first-principles (FP) molecular-dynamics simulations (MD). In MD the atoms move according to the Newtonian mechanism under the action of the interatomic forces. We compute these forces using the Projection-Augmented Wavefunction (PAW) implementation of the density-functional theory as in the VASP package. We employ the PBE GGA functional for the exchange-correlation part. Our simulations are performed in NVT ensembles, i.e. the Number of particles and the Volume of the cell are kept fixed. The Temperature is controlled by the use of external thermostats to oscillate around a desired average value.

In-house software
We developed a post-processing utility, described at point 1.2 which will be made available to the entire computational community, be it in geosciences, in physics, or in chemistry. We are currently in the final stage of source development and will be able to release a first complete version and submit a publication before the summer of 2019.
The UMD package consists of:
A parser between various ab initio molecular dynamics codes (VASP, QBox, abinit, more to come at a later stage) and the UMD format file
A tool to compute the pair distribution function
A tool to compute the speciation and clusterization based on an interconnectivity matrix (the atoms are bonded if they lie inside their corresponding first coordination sphere)
A tool to distinguish “gas”-like and “liquid”-like components
A tool to compute the gas vs liquid volume ratio
A tool to compute the atomic mean square deviation, from which diffusion coefficients are obtained
A tool to compute the diffusion coefficients of various chemical species identified during the speciation analysis
A tool to compute self-correlation properties (atomic velocities to build the vibrational spectrum, stresses to extract the viscosity)

Mineralogical systems

For each of these systems we build the Maxwell construction, obtain the liquid spinodal as a function of temperature and estimate the position of the critical point. Then we analyse the structural properties of the fluids: interatomic bonding, atomic coordination, clustering, and polymerization, liquid vs gas phase. Second we look at transport properties: mean square displacements and diffusion coefficients. Finally we try to obtain the entropy and free energies of the fluid Depending on the systems we also analyse the magnetic state, and the electrical conductivity.

The bulk silicate Earth
We consider the composition of the bulk silicate Earth as proposed by McDonough and Sun (1995), which we call hereinafter “pyrolite”. We cover the 0.75 – 7.5 g/cm3 density range and 2000 – 10000 K temperature range. We find the critical point above 6000K at densities below
At high density, the liquid is highly polymerized and viscous. At low density and low temperatures, in the 2000 to 4000 K range, we capture the nucleation of bubbles. These bubbles contain a low-density gas phase dominated by atomic species, alkaline and calc-alkaline cations, and SiOx groups. When volatiles are present such molecular species are the first ones to evaporate into these bubbles.
At high temperature, we identify the supercritical region characterized by one homogeneous fluid, rich in ionic species with short lifetimes. Its chemical speciation is very different from the one at ambient pressure conditions.
At low pressures all 4 Fe atoms are in a high-spin state almost the entire simulation. With increasing pressure, there is a statistical quantity of instances of lower Fe atoms. This population steadily increases with increasing pressure. At the highest pressure points the melt can be considered in a low spin state; at conditions typical of the magma ocean, we can safely assume the presence of the residual magnetic moment on the Fe atoms in the silicate melt.
The nucleation of the bubbles corresponds to the spinodal instability of the melt phase in the van der Waals description of the liquid-gas equilibrium. The spinodals are reached consistently regardless of the thermodynamic path we chose to obtain the low densities. The critical curves, inferred from the spinodals, are necessary to understand the separation and degassing of volatiles during the recovery from a giant impact. As such we obtain that the largest part of the disk was in the supercritical state for a long time. The condensed central Earth inside synestia doesn’t have a “surface”.

Pure iron
We performed a full calculation for the critical point of iron. Our lowest temperature is set as 3000 K, close to the boiling point at 1 bar. The highest temperature is 12000 K. At ambient pressure conditions, i.e. P = 0 GPa, and 3000K, the liquid density is 6.2 g/cm3, compared with experimental values of 6.22 g/cm3. At the same pressure but 4000K, the liquid density decreases to 5.7 g/cm3, only 0.2 g/cm3 larger than the experimental value of 5.5 g/cm3.
For iron, we observe the existence of pressure minima for all isotherms up to 9000 K. Along this isotherm, the minimum pressure is obtained at a density of about 2.25 g/cm3, and the maximum of the pressure at a density of 2.05 g/cm3. Starting with the 9350 K isotherm and above, with decreasing density, the pressure monotonically decreases; this is characteristic to the supercritical state. Therefore, the critical density must be found between 2.05 g/cm3 and 2.25 g/cm3 and between 9000 K and 9350 K.
The position of the critical point of iron that we find in our simulations is low enough in density for the Earth's core to be in a supercritical state or in a molten state throughout the giant impact and the condensation of the protolunar disk. Only iron ejected in the outermost parts of the disk, like the one originating from the impactor's core can reach the liquid-vapor dome. Indeed, the temperatures in these outer regions of the disk or of the synestia are low enough not to exceed the critical temperature.
Then we compute the entropy along the spinodal and the Hugoniot lines. The energy - pressure - volume data allow us to compute the Hugoniot lines via Rankine-Hugoniot equations. The warm Hugoniot line with initial conditions at 1GPa and 1500 K intercepts the iron melting curve at 130 GPa. The hot Hugoniot line is calculated with an initial conditions at 40 GPa and 4000K, which may be a representative condition for Mars-size impactor. For peak shock conditions at 15000K, we estimate that the entropy can reach up to 19.1 kB/atom and 18.6 kB/atom along the warm and the hot Hugoniot lines respectively. Our results suggest that above 7500 K, the entropy along the warm Hugoniot line is less than that along the hot Hugoniot line. We note that the entropy difference between these two Hugoniot lines is relatively small (0.5 kB/atom) from 7500 K to 15000K. Based on the computed entropy along the warm Hugoniot line, the shock pressure required to reach the entropy necessary for the onset of vaporization upon release and cooling is 312 GPa.
The entropies that we obtain allowed us to conclude that in most impacts with planetesimals their core would undergo partial vaporization. A manuscript is submitted to Science Advances with these results.

Feldspars
We consider the end-member terms of the feldspar class of minerals: KAlSi3O8, NaAlSi3O8, and CaAl2Si2O8. We sampled numerically a wide range of temperatures and pressures: from 2000 to 7000 K and from 0.5 to 6.4 g/cm3, which cover the liquid side of the liquid-vapor dome and the position of the critical point.
The cations are the most diffusive elements at low densities. The self diffusivities of Al, O, and Si are similar along each isotherm. The difference between isotherms is reduced when the temperature increases. At high temperature and low density, the self diffusivity of every elements tends toward 1-2.10−7m2s−1.
At 3000 K, species with an even number of oxygen (SiO4, SiO6, AlO4, AlO6, AlO8) reach a maximum proportion bigger than species with an odd number of oxygen (SiO5,244AlO5, AlO7). The peak for SiO4 and AlO4 is located around 2.2g cm−3for CaAl2Si2O8 and KAlSi3O8 and around 2.1g cm−3 for NaAlSi3O8. For lower densities the proportion of 4 and 5 fold coordination decrease whereas the proportion of SiO3 and AlO3increase. At 6000 K the general behavior is the same but no species reach a proportion higher than 0.5. Species with an odd number of oxygen are as major as species with an even number of oxygen. Moreover, when the density decreases, very small species as SiO, AlO, SiO2 and AlO2appear, whereas at 3000 K we only saw SiO3 and AlO3.
At 4000 K, a few number of species appear in the gas phase (populating the bubbles), with clusters living in average between 60 and 300 fs. The major species are isolated Na, SiO, SiO2 and NaSiO3. At 6000 K we see a lot more species, especially with 4 to 13 atoms in their formula. We see less isolated Na but more isolated O, along with a276huge quantity of the species seen at 4000 K. We also see O2 molecules living up to 550 fs.
We find the critical points below 0.8 g/cm3, at 4000-4500 K for K- and Na- terms and at 6000-6500 K for Ca term, thus confirming the refractory character of Ca.

Volatiles in the bulk silicate Earth
We consider different series of simulations where we add He, CO, CO2, H2O, CO2+He, CO2+H2O, CO+H2O, Zn, and In to the pyrolite. The volatiles can either be added as molecules or as replacement reactions. In the second case we have:
Mg ó 2H+
Mg ó Zn
Al ó In
Just as for the pure dry pyrolite we monitor the position of the spinodals. But additionally we look for the speciation of the volatiles in the melt, and the formation of the bubbles.
The simulations with pyrolite + CO and pyrolite + CO2 revealed a far more complex chemical behaviour of carbon than expected. At 0 GPa and 3000–5000 K, carbon exists exclusively in unpolymerized forms (e.g. CO, CO2, and CO3) while at higher pressures, carbon becomes increasingly polymerized (i.e. higher fraction of carbon atoms bonded to another carbon atom). We found that the more reduced pyrolite melts have a larger quantity of polymerized carbon with respect to the more oxidized pyrolite melts with similar carbon concentrations. Most of the polymerized carbon exists in C2O, C2O2, and C2O3 molecules at ~65 GPa, whereas at 115 GPa, carbon exists predominantly as C2O2–6 and C3O4–6 molecules. The oxidation state has a significant effect on the types of species present. The more oxidized pyrolitic melts only contain C2O4, C2O5, C2O6, and C2O7, with nearly no C3Ox species forming while the more reduced melts have a wide range of CxOy species present. These species are long-lived. Chemically they represent organic radicals. We can say that we observed the formation of oxalates in the deep magma ocean.
Moreover we observed strong clustering between carbon and iron, which increases with increasing pressure. For all the melt compositions, nearly all carbon clusters with three or more carbon atoms contain iron above 50 GPa. Our results also confirmed that carbon bonds directly to silicon, substituting oxygen. The degree of bonding between carbon and silicon seems to be largely independent of temperature and carbon concentration.
As we decrease the density, we observe the formation of bubbles or cavities, which contain volatile species. At 3000-4000 K, the cavities contain almost exclusively CO and CO2 species whereas at 5000 K, a rich variety of species, such as O, O2, SiO, Mg and Na enter the cavities.
Most of the carbon is retained in the melt at a density of 2 g/cc; however, as the density approaches 1 g/cc, about half of the carbon exists in the cavities rather than the molten silicate phase. The fact that a significant fraction of carbon is retained in the melt implies that carbon may not be as volatile at the densities and temperatures of the proto-lunar disk as previously thought. Consequently, Earth may have not lost all of its carbon during the impact. Instead, the affinity of carbon for iron implies that in the large magma oceans of the proto-Earth and Early Earth (prior to and after the Moon-forming impact) carbon would most probably follow iron during segregation of an iron metallic liquid, and a majority of carbon present in the deep mantle may have been transported into the iron-rich core. After the moon-forming impact, metal iron and silicate melt would have been well-mixed and subsequent core segregation would have resulted in siderophile elements being sequestered into the Earth's core.

* We have refined the Maxwell construction to obtain critical points for a variety of mineralogical systems
* We have implemented the method of the Gibbs ensemble to determine the vapor -liquid equilibrium. According to this procedure we define two simulation boxes, one at high density, i.e. the liquid, and one at low density, i.e. the vapor. We allow transfer of atoms between the two boxes with the Boltzmann probability based on the energy of the systems. At equilibrium the probability of the transfer is the same in the two directions. Given appropriate computational resources, this procedure can be used also to find the chemical equilibrium and the element partitioning between two geochemical reservoirs. This is a methodological development that we adopted during the project.
* We expect to obtain the critical points of a large majority of mineralogical systems, in agreement with what was initially planned for the project. But we will also apply the Gibbs ensemble method to a series of systems, including multi component silicate melts, like pyrolite, to compute element partitioning between vapor and liquids, or between two liquids
* We will fully characterize the structural and the transport properties of a large series of melts and supercritical fluids that are relevant for the early Earth

Schematic 3-steps history of disk condensation. Initial state, shock state, isoentropic condensation

Periodic Reporting for period 4 - IMPACT (The giant impact and the Earth and Moon formation)

Share this page

Download