Thierry Alboussière

Eric Debayle

Isabelle Daniel

Anaïs

Alistair Gabriel

Elham Ali Anastasis

Tim Zhi

The Moon Impact formation theory was first proposed in 1975 by Hartmann and Davis. To this day, it is still the prevalent theory to explain the unique aspects concerning the Earth-Moon system. Although well accepted in the scientific community, this theory is still subject of debate regarding the conditions of the impact and the resulting scenario. It is still challenging to encompass in a single solution the problematic aspects of chemical equilibration that took place post-impact. Different models of impact have been proposed over the years in order to address the aforementioned issues. These models often rely on equations of state to describe the behavior of the materials present in the Earth and the Moon. Equations of state (EOS) describe the distribution of materials' phases and they are used on collision models based in hydrodynamic simulations to predict the final composition. The use of precise EOS lead to more correct post-impact chemical models where phases can be estimated correctly.

The majority of minerals present on rocky bodies in the Solar System are SiO 2 based and they are the building blocks of the Earth and the Moon. In terms of composition, SiO 2 represents more than 44% of lunar mare basalts and highlands, and about 45% of Earth's primitive mantle. However similar in parts of their composition, Moon and Earth differ considerably when it comes to the presence of volatile elements, where H 2 O is arguably the most relevant. The impart of energy given by the large impact is sufficient to melt and vaporize silicate minerals. It produces a silicate atmosphere, where physicochemical properties and geochemical signatures of the resulting scenario depend on the fractionation of liquid-vapor. The presence of volatiles can change the behavior of typical equations of state of silicate materials and eventually be responsible for increasing the amount of vapor and supercritical material in such conditions.

We aim to construct equations of state for the silica-water binary system and subsequently study structural and transport properties of such system, comparing the behavior of the main elements under different conditions. We employ ab initio molecular dynamics (AIMD) methods where forces and energy are calculated with density functional theory as implemented in the VASP ® package. We start by building a supercell containing 72 molecules of SiO 2 and describe its behavior under temperatures from 4000 to 7000 K and densities from 0.2 to 2.33 g/cm 3 . We follow by inserting 9, 18, 36 and 72 H 2 O molecules inside the SiO 2 supercell to obtain four different water-silica ratios and compare changes in the thermodynamic behavior as well as potential effects in elemental structural and transport properties. We study the four systems at temperatures of 2000, 3000, 4000 and 5000 K and densities of 0.34 to 2.77 g/cm 3 depending on the system. We calculate the critical point of SiO 2 in between 5000 and 5500 K, from 0.6 to 0.85 g/cm 3 and 0.15 to 0.25 GPa. The presence of water has a direct effect on the placement of the critical point, reducing it in more than 2000 K at specific conditions. Our results show that supercritical state materials are often underestimated in moon forming hydrodynamic models, which may impact elemental mixing in the outcome of the Giant Impact.

A thesis is the result of the collective work of many people.

Résumé

La théorie de la formation de la Lune par un impact géant a été proposée pour la première fois en 1975 par Hartmann et Davis. À ce jour, c'est toujours la théorie qui prévaut pour expliquer les aspects uniques concernant le système Terre-Lune. Bien qu'elle soit bien acceptée dans la communauté scientifique, cette théorie fait encore débat quant aux conditions de l'impact et au scénario qui en résulte. Il est toujours difficile d'intégrer au sein une solution unique les aspects problématiques de l'équilibrage chimique qui ont eu lieu après l'impact. Différents modèles d'impact ont été proposés au fil des années afin de résoudre les problèmes susmentionnés. Ces modèles s'appuient souvent sur des équations d'état (EOS) pour décrire le comportement des matériaux présents sur la Terre et la Lune. Les EOS décrivent la distribution des phases des matériaux et sont utilisées sur des modèles de collision basés sur des simulations hydrodynamiques pour prédire une composition finale. L'utilisation d'EOS précises conduit à des modèles chimiques post-impact plus exacts pour lesquels l'estimation des phases est plus rigoureuse.

La majorité des minéraux présents sur les corps rocheux du Système solaire sont à base de SiO 2 et sont les éléments constitutifs de la Terre et de la Lune. En termes de composition, SiO 2 représente plus de 44% des mare basaltes et des hauts plateaux lunaire, et environ 45% du manteau primitif de la Terre. Bien que similaires sur certains aspects de leur composition, la Lune et la Terre diffèrent considérablement en ce qui concerne la présence d'éléments volatils, où H 2 O est sans doute le plus pertinent. L'apport d'énergie donné par le grand impact est suffisant pour faire fondre et vaporiser les minéraux silicatés. Il produit une atmosphère silicatée, dont les propriétés physico-chimiques et les signatures géochimiques du scénario résultant dépendent du fractionnement liquide-vapeur. La présence de volatils peut modifier le comportement des équations d'état typiques des matériaux silicatés et éventuellement être responsable de l'augmentation de la quantité de vapeur et de matériau supercritique dans de telles conditions.

Notre objectif est de construire des équations d'état pour le système binaire silice-eau afin d'en étudier les propriétés structurelles et de transport, en comparant le comportement des principaux éléments dans différentes conditions. Nous utilisons des méthodes de dynamique moléculaire ab initio (AIMD) où les forces et l'énergie sont calculées selon la théorie fonctionnelle de la densité telle qu'implémentée dans le package VASP. Nous commençons par construire une supercellule contenant 72 molécules de SiO 2 et décrivons son comportement sous des températures s'étendant de 4000 à 7000 K et des densités allant de 0,2 à 2,33 g/cm 3 . Nous insérons 9, 18, 36 et 72 molécules H 2 O à l'intérieur de la supercellule SiO 2 pour obtenir quatre rapports eau-silice différents et comparer les changements dans le comportement thermodynamique ainsi que les effets potentiels sur les propriétés structurelles et de transport élémentaire. Nous étudions les quatre systèmes à des températures de 2000, 3000, 4000 et v 5000 K et des densités allant de 0,34 à 2,77 g/cm 3 selon le système. Nous calculons le point critique de SiO 2 entre 5000 et 5500 K, de 0,6 à 0,85 g/cm 3 et de 0,15 à 0,25 GPa. La présence d'eau a un effet direct sur le placement du point critique, le réduisant à plus de 2 000 K dans des conditions spécifiques. Nos résultats montrent que les matériaux à l'état supercritique sont souvent sous-estimés dans les modèles hydrodynamiques de formation de lune, ce qui peut avoir un impact sur le mélange élémentaire dans le résultat de l'impact géant.

To my mother. May her laughter echo forever.

Solar System Overview

The Solar System is largely formed by hydrogen and helium, where all the other elements represent only 2% of the total elemental abundances. The solar mass represents nearly 99% of the total mass and represents the mean composition of the Solar System [START_REF] Lodders | Solar elemental abundances[END_REF]. The remaining ∼1% of the mass is distributed in the form of planets, satellites, asteroids, comets and other diverse bodies. The diversity of objects in the Solar System seems in contrast with its actual composition, especially considering that they are assumed to originate from the same initial composition [START_REF] Taylor | Planetary crusts: their composition, origin and evolution[END_REF]. However, a different theory based on titanium isotopes suggests that thermal processing in the disk before the period of accretion lead to condensation from a non-solar composition [START_REF] Trinquier | Origin of nucleosynthetic isotope heterogeneity in the solar protoplanetary disk[END_REF]. Nonetheless, the differences in size, structure and composition observed throughout the Solar System are a result of formation and evolutionary processes in the protoplanetary disk.

The understanding of processes involved during the formation of planetary bodies and their satellites are still cause of debate. In the protoplanetary disk, radial variations in physicochemical conditions affected the stability of the elements and altered the composition of the disk with increasing solar distance [START_REF] Russell | Timescales of the solar protoplanetary disk[END_REF], and references therein), as shown in Figure 1.1. Far from the radiation and high temperatures of the Sun, the outer part of the disk had cryogenic temperatures. Temperature estimates for this region can vary from 50 to 150 K [START_REF] Boss | Temperatures in protoplanetary disks[END_REF]. In these conditions, volatile elements condensed and resulted in the formation of icy planets and moons [START_REF] Chambers | Planet Formation[END_REF]. However, closer to the Sun, the inner part of the disk had higher temperatures, reportedly varying from 500 to 1500 K [START_REF] Boss | Temperatures in protoplanetary disks[END_REF], where evidence suggests silicates to be partially vaporized [START_REF] Chambers | Planet Formation[END_REF][START_REF] Messenger | The population of starting materials available for Solar System construction[END_REF]. In this region, the protoplanetary disk was rich in silicate material and refractory elements, originating essentially rocky bodies. Accretion in the inner part of the disk happened in several steps. The direct condensation of micrometer sized grains from the gas, together with the later processing of these grains, lead to the formation of the oldest solids known, such as calcium-aluminum inclusions (CAIs), ameboid-olivine aggregates (AOA) and chondrules [START_REF] Russell | Timescales of the solar protoplanetary disk[END_REF], which are all present in carbonaceous chondrites. Condensation and evaporation were dominant processes in the first stages of accretion and happened much more rapid than in later stages, where larger bodies and planetesimals were already formed. Planetary embryos underwent a period of stochastic events of accumulation, often involving giant impacts. At the same time, these bodies were subject to secondary processes in the form of aqueous alteration and thermal metamorphism, which involved transformation of the mineralogical phases already present. Furthermore, they also underwent a certain degree of melting and differentiation, resulting in the formation of cores and layered structures.

Models of accretion still struggle to yield an accurate description of the current structure of the Solar System and also depend on types of collisions considered as a mechanism in place during this period [START_REF] Chambers | Late-stage planetary accretion including hit-and-run collisions and fragmentation[END_REF]. Nevertheless, it is currently understood that planets, asteroids and the more than 200 satellites formed as a result of surrounding conditions regarding their positioning in the protoplanetary disk and the posterior processes that took place. The unique combination of such processes resulted not only in different sizes and structures, but, most importantly, in different chemical signatures for every single parent body in the Solar System. The study of the imprinted signatures shed light into the evolution of planetary bodies and relied strongly on materials available for in situ analysis, which for many years were made essentially by meteorites, micrometeorites and interplanetary dust particles [START_REF] Grady | Types of extraterrestrial material available for study[END_REF].

With the beginning of the Space Age in 1957 and subsequent space missions up to the present-day, samples could be directly gathered from the parent body and returned to Earth. The Apollo program was responsible for returning more than 300 kg of lunar material and lead to remarkable discoveries on the isotopic signature of the Moon and to breakthroughs regarding the Earth-Moon system. The analysis of lunar oxygen isotopic composition revealed to be surprisingly nearly identical to that of the Earth, calling for new models of formation of the Earth-Moon system that could explain such similarity [START_REF] O'neil | Oxygen isotope analyses of selected Apollo 11 materials[END_REF][START_REF] Clayton | Oxygen isotopic compositions of Apollo 15, 16, and 17 samples, and their bearing on lunar origin and petrogenesis[END_REF]. Questions on the nature of the processes involved during the formation of the Moon remain unclear to this day.

The Earth-Moon system

Inner planets and satellites of the Solar System are largely composed of silicate materials, where SiO 2 is the major building block. The silicate portion often envelopes a metallic core made essentially of Fe-Ni alloy [START_REF] Li | Experimental constraints on core composition[END_REF]. The formation of the core is a result of the silicate-metal fractionation during the period of accretion, when the energy imparted by impacts and decay of radioactive elements resulted in high temperatures that facilitated the formation of magma oceans (Elkins-Tanton, 2012; [START_REF] Day | Evaporative fractionation of volatile stable isotopes and their bearing on the origin of the Moon[END_REF]. At these conditions, denser materials displace towards the nucleus and lighter materials ascent towards the surface, resulting in a fully differentiate body. It is also at this stage that large-scale degassing can happen, leading to significant loss of volatile elements [START_REF] Elkins-Tanton | Magma oceans in the inner Solar System[END_REF]. The silicate fraction can be separated into a deeper portion, called mantle, and a more superficial part, called crust. Because of the increasing pressure, minerals in the mantle tend to present a tighter packing than those found in the crust. Furthermore, the crust is exposed to superficial secondary processes, such as cratering and, in systems that have an atmosphere, these processes can also include alteration by weathering [START_REF] Taylor | Planetary crusts: their composition, origin and evolution[END_REF].

The Earth and the Moon are differentiated bodies with crust, mantle and core. The primitive silicate composition of the Earth was established by [START_REF] Mcdonough | The composition of the Earth[END_REF] and is referred to as the Bulk Silicate Earth (BSE). This estimate defines the composition of the silicate portion before mantle-crust differentiation and after core formation [START_REF] Palme | Cosmochemical estimates of mantle composition[END_REF]. Formed shortly after accretion, Earth's core has more than 3000 km of radius and represents roughly 16% of the total volume [START_REF] Li | Experimental constraints on core composition[END_REF], as represented in Figure 1.2a. With respect to solar abundances, the Earth is depleted in moderate volatile and volatile elements [START_REF] Halliday | The origin and earliest history of the Earth, Treatise on Geochemistry[END_REF]. A diagram describing the variation on terrestrial elemental concentrations with respect to carbonaceous chondrites composition values in [START_REF] Anders | Abundances of the elements: Meteoritic and solar[END_REF] is shown in Figure 1. 3. Lunar composition is less accessible in comparison to the Earth and it depends on lunar samples and seismic data. Seismic velocity measurements show that the lunar core is signifi- Figure from [START_REF] Li | Experimental constraints on core composition[END_REF]. b) Estimate internal structure of the Moon. The solid inner core is represented in orange with radius of 240 km, predicted fluid outer core is indicated in yellow with a radius of 330 km, partial melt enveloping the core is in dark yellow with a radius of 480 km, mantle is indicated in green and crust in grey. Figure adapted from [START_REF] Weber | Seismic detection of the lunar core[END_REF].

cantly smaller than the terrestrial core [START_REF] Weber | Seismic detection of the lunar core[END_REF] and represents about 2 to 3% of its volume [START_REF] Taylor | Planetary crusts: their composition, origin and evolution[END_REF], as shown is Figure 1.2b. Such evidence when associated with the mean density estimate for the Earth (5.54 g/cm 3 ) and the Moon (3.55 g/cm 3 ), indicate a depletion on lunar contents of iron with respect to that of the Earth [START_REF] Taylor | Planetary crusts: their composition, origin and evolution[END_REF]. Furthermore, the Moon also shows depletion in moderate volatile and volatile elements when compared to terrestrial values (Lock et al., 2018, and references therein). Conversely, measurements of refractory elements on lunar samples show higher contents than in the BSE.

A compilation of bulk silicate Moon compositional values can be found in [START_REF] Melosh | New approaches to the Moon's isotopic crisis[END_REF]. Nevertheless, the satellite-to-parent mass ratio of the Earth-Moon system is the largest observed in the Solar System [START_REF] Taylor | Planetary crusts: their composition, origin and evolution[END_REF]. The only other system where orbiting satellite and parent body show a larger ratio is Pluto-Charon.

The uniqueness of the Moon does not come solely from being the only satellite orbiting the Earth. Besides the compositional differences observed, the Earth and the Moon share unique aspects that are not generally observed elsewhere in the Solar System. The roughly identical isotopic composition on many stable isotopes is a key feature that constraint models of lunar formation and need to be addressed when modeling different scenarios.

Isotopic composition

First in situ measurements of oxygen isotopes in lunar samples revealed a virtually identical composition with respect to terrestrial values [START_REF] O'neil | Oxygen isotope analyses of selected Apollo 11 materials[END_REF][START_REF] Clayton | Oxygen isotopic compositions of Apollo 15, 16, and 17 samples, and their bearing on lunar origin and petrogenesis[END_REF] and similar values continue to be reported by following authors [START_REF] Clayton | Genetic relations between the Moon and meteorites[END_REF][START_REF] Wiechert | Oxygen isotopes and the Moon-forming giant impact[END_REF][START_REF] Hallis | The oxygen isotope composition, petrology and geochemistry of mare basalts: Evidence for large-scale compositional variation in the lunar mantle[END_REF]. A discrete deviation in ∆ 17 O values has been described in the literature, where authors argue that this divergence can originate from a contribution from late veneer impactors [START_REF] Greenwood | Oxygen isotopic evidence for accretion of Earth's water before a high-energy Moon-forming giant impact[END_REF]. Other authors argue that ∆ 17 O variations can be a result of condensation of SiO 2 in the presence of hydrogen from a silicate atmosphere in the aftermath of a giant impact [START_REF] Cano | Distinct oxygen isotope compositions of the Earth and Moon[END_REF]. However, this deviation is not in agreement with values previously reported in the literature [START_REF] Wiechert | Oxygen isotopes and the Moon-forming giant impact[END_REF]. A comprehensible review on the oxygen isotopic composition of Solar System bodies has been published by [START_REF] Ireland | Oxygen isotopes and sampling of the Solar System[END_REF], showing the contrasting isotopic signatures on analytically accessible bodies, such as meteorites, asteroids, satellites and planets.

Further analysis on lunar samples showed persistent identical values not only for oxygen isotopes, but on a range of stable nuclides, such as titanium [START_REF] Zhang | The proto-Earth as a significant source of lunar material[END_REF]), silicon (Armytage et al., 2012) and chromium [START_REF] Mougel | Chromium isotopic homogeneity between the Moon, the Earth and enstatite chondrites[END_REF]. As reported on titanium and chromium measurements, the exposure to cosmic rays generates secondary effects that can modify the isotopic composition and, thus, corrections on such effects can broaden the range of near identical values for the Earth-Moon isotopic composition. However, the similarity in the isotopic signature is not observed for some of the short-lived nuclides that are sensitive to differentiation and planetary formation processes, such as the 182 Hf-182 W system [START_REF] Kruijer | Tungsten isotopes and the origin of the moon[END_REF][START_REF] Jacobsen | Accretion and Early Differentiation History of the Earth Based on Extinct Radionuclides[END_REF][START_REF] Jacobsen | The Hf-W isotopic system and the origin of the Earth and Moon[END_REF]. Terrestrial mantellic heterogeneities have been observed in tungsten isotopes [START_REF] Touboul | 182 W Evidence for long-term preservation of early mantle differentiation products[END_REF] and new measurements on lunar vanadium isotopes indicate a chondritic component compared to the bulk silicate Earth [START_REF] Nielsen | Isotopic evidence for the formation of the Moon in a canonical giant impact[END_REF].

The isotopic similarities in the Earth-Moon composition require a model with a common origin for both bodies. Mechanisms that have been described in order to explain such observations consider either the inheritance of nuclides from an identical isotopic reservoir in the protoplanetary disk, or the mixing and equilibration as a result of a giant impact [START_REF] Pahlevan | Equilibration in the aftermath of the lunar-forming giant impact[END_REF][START_REF] Canup | Forming a Moon with an Earth-like composition via a giant impact[END_REF][START_REF] Cuk | Making the Moon from a fast-spinning Earth: A giant impact followed by resonant despinning[END_REF][START_REF] Lock | The origin of the Moon within a terrestrial synestia[END_REF]. Given the heterogeneity of the protoplanetary disk in terms of chemical distributions, the simultaneous accretion from an identical reservoir is unlikely. Therefore, models that consider the formation of the Earth-Moon system from a giant impact remain the best hypothesis to this day to explain the identical isotopic signatures.

Water

Regardless of the similarity in isotopic compositions, Moon and Earth differ considerably concerning the abundance of volatile species. Their contents can differ concerning sodium, potassium and others (Melosh, 2014, and references therein). Being the basis of life as we know, H 2 O is arguably the most relevant volatile specie to be consider. The presence of water affects properties of materials, such as viscosity and melting, but also can give evidence on the origin and formation of planetary bodies [START_REF] Day | Evaporative fractionation of volatile stable isotopes and their bearing on the origin of the Moon[END_REF]. Water is ubiquitous on Earth and measurements of oxygen isotopes on zircon crystals suggest the presence of a terrestrial water ocean as early as 4.4 Gyr [START_REF] Wilde | Evidence from detrital zircons for the existence of continental crust and oceans on the Earth 4.4 Gyr ago[END_REF][START_REF] Valley | A cool early Earth[END_REF]. The origin of the water on Earth is still not entirely understood. However, studies show that more than 70% of water present today was accreted prior to the Moon forming impact and up to 30% from later chondritic contribution during the late veneer [START_REF] Greenwood | Oxygen isotopic evidence for accretion of Earth's water before a high-energy Moon-forming giant impact[END_REF]. With the lack of visible liquid water or ice caps, the Moon was thought to be completely anhydrous for many years [START_REF] Taylor | Where has all the lunar water gone?[END_REF]. The improvement in instrumental accuracy made possible to detect H 2 O contents on lunar rocks, revealing a more hydrated Moon than previously thought [START_REF] Hauri | Water in the Moon's interior: Truth and consequences[END_REF]. Depending on the lithological type, lunar H 2 O can vary from 3 to 292 ppm [START_REF] Canup | Origin of the Moon[END_REF].

In general, volatile loss can be a result of degassing during volcanic eruptions or during a magma ocean stage. [START_REF] Saal | Volatile content of lunar volcanic glasses and the presence of water in the Moon's interior[END_REF] analyzed lunar glass beads in sample cross section, where results showed an increasing depletion in H 2 O contents and further volatiles species from the center towards the rims, indicating significant degassing during volcanic eruptions. Furthermore, the authors also argue that a deep layer of ilmenite may hold volatiles during the crystallization of the lunar magma ocean (LMO) and prevent a complete volatile loss. Conversely, [START_REF] Day | Evaporative fractionation of volatile stable isotopes and their bearing on the origin of the Moon[END_REF] argue that volcanic degassing is not sufficient to explain the water contents and that the volatile depletion can only be explained by degassing during the lunar magma ocean. Nevertheless, water contents can be expected to increase towards the end of crystallization of magma oceans (Figure 1.4;[START_REF] Elkins-Tanton | Formation of early water oceans on rocky planets[END_REF], meaning also an increase of its effects on the residual melt. The inclusion of water in models of lunar formation is important in order to provide a better understanding of the processes involved (Saal et al.,

Lunar formation models

Before Apollo missions, the Moon was said to be cold and primitive. As a general rule, satellites can accrete from an early debris disk, but they can also be captured objects. Early models of lunar formation proposed before Apollo missions were fission [START_REF] Darwin | XIII. on the precession of a viscous spheroid, and on the remote history of the Earth[END_REF], capture [START_REF] Urey | The capture hypothesis of the origin of the Moon, The Earth-Moon System[END_REF] and coaccretion [START_REF] Schmidt | A Theory of Earth's Origin[END_REF]. In the fission model, a rapidly spinning Earth would become so elongated that would eventually rupture, giving origin to the Moon from the separated fragment. However, it fails to explain the origin of such high angular momentum and what would be the processes involved to dissipate it afterwards. The model proposed by [START_REF] Urey | The capture hypothesis of the origin of the Moon, The Earth-Moon System[END_REF], considers that the Moon accreted in a different region of the protoplanetary disk and was later captured into the terrestrial orbit. The capture theory had difficulties conciliating the velocities involved and the mechanisms of capture. With the isotopic measurements in the following years, this theory was discarded. The leading most accepted theory was croaccretion [START_REF] Schmidt | A Theory of Earth's Origin[END_REF][START_REF] Wood | Moon over Mauna Loa: a review of hypotheses of formation of Earth's Moon[END_REF], where the Moon and the Earth were believed to have accreted from the same protoplanetary disk region. Even though this model would agree with posterior isotopic observations, it failed to explain the differences in bulk compositional values, angular momentum and evidences of existence of a magma ocean [START_REF] Wood | Moon over Mauna Loa: a review of hypotheses of formation of Earth's Moon[END_REF]). Apollo's returned samples shed light into the composition of the Moon and showed a differentiated body with igneous activity. Since the first oxygen isotopic measurements made on lunar returned samples, considerable progress was made to address the questions relative to the formation of the Moon. It became clear that the Moon and the Earth must had shared a common origin. Seismic data and further analytical measurements, continue to reveal characteristics that need to be considered and explained by any lunar formation model. Large impacts were already believed to be a mechanism of accretion occurring throughout the Solar System and were also responsible for melting processes after the extinction of 26 Al (Elkins-Tanton, 2012). To this day, the theory that seems to better address the singular aspects of the Earth-Moon system is the Giant Impact Theory.

Giant Impact Theories

In the Giant Impact Theory, the Moon was formed by accretion from a dust and gas disk formed in the aftermath of a large collision between the proto-Earth and an impactor. This theory was first proposed by [START_REF] Hartmann | Satellite-sized planetesimals and lunar origin[END_REF], followed by [START_REF] Cameron | The origin of the Moon[END_REF] and has continued to evolve ever since. Different scenarios mostly based on smoothed particles hydrodynamic simulations have been invoked over the years (e.g. [START_REF] Hartmann | Satellite-sized planetesimals and lunar origin[END_REF][START_REF] Canup | Forming a Moon with an Earth-like composition via a giant impact[END_REF][START_REF] Cuk | Making the Moon from a fast-spinning Earth: A giant impact followed by resonant despinning[END_REF][START_REF] Reufer | A hit-and-run giant impact scenario[END_REF][START_REF] Rufu | A multiple-impact origin for the Moon[END_REF][START_REF] Lock | The origin of the Moon within a terrestrial synestia[END_REF], challenged to explain the Moon's large mass, orbit inclination, isotopic composition, volatile budget and iron depletion within a single solution. They consider different initial settings which affect possible final elemental distributions in the protolunar disk. Many reviews on the Moon formation problematic have been published recently [START_REF] Asphaug | Impact origin of the moon?[END_REF][START_REF] Barr | On the origin of Earth's moon[END_REF][START_REF] Lock | Geochemical constraints on the origin of the moon and preservation of ancient terrestrial heterogeneities[END_REF][START_REF] Canup | Origin of the Moon[END_REF] with an accurate description of the solutions presented and also the potential inadequacies of each model. Here, we are going to describe briefly two current models of impact found in the literature: the canonical impact [START_REF] Canup | Forming a Moon with an Earth-like composition via a giant impact[END_REF] and the fast-spinning Earth model [START_REF] Cuk | Making the Moon from a fast-spinning Earth: A giant impact followed by resonant despinning[END_REF][START_REF] Lock | The origin of the Moon within a terrestrial synestia[END_REF].

Canonical Impact model

The canonical impact model generally considers a low-velocity oblique collision by a large impactor [START_REF] Canup | Lunar-forming impacts: processes and alternatives[END_REF]. It has been described extensively in the literature with varying initial simulation settings [START_REF] Canup | Origin of the moon in a giant impact near the end of the Earth's formation[END_REF][START_REF] Canup | Dynamics of lunar formation[END_REF][START_REF] Salmon | Lunar accretion from a roche-interior fluid disk[END_REF][START_REF] Canup | Forming a Moon with an Earth-like composition via a giant impact[END_REF]. The model proposed by [START_REF] Canup | Forming a Moon with an Earth-like composition via a giant impact[END_REF] starts from the collision between two differentiated bodies and consider a much larger impactor than previous studies, where the mass of the impactor is similar to the target's mass. The scenarios with most fitting results consider an impactor-to-total mass ratio between 0.4 and 0.5, where the impactor contributes in similar proportions to the disk and to the target. The collision creates a disk with 50 to 90% of its mass in the vapor phase where the material comes mostly from the impactor's mantle. The composition of the disk is homogeneous and depleted in iron. The impactor composition was taken to be similar to Mars and results showed that the final silicate composition diverges less than 1%. Temperatures in this model can vary roughly from 2000 K to more than 6400 K in an 80,000 km across disk. The total simulated time was 1 day and, once particles are in the disk, they can accumulate at a distance of 3.8 Earth's radius to form the Moon. In a concomitant work, the period necessary for the complete accretion of the Moon from the disk was estimated in about 100 years [START_REF] Salmon | Lunar accretion from a roche-interior fluid disk[END_REF]. Oxygen represents the biggest constraint in previous canonical impact models, since most disk material originates from the impactor [START_REF] Canup | Origin of the Moon[END_REF]. However, the large mass proposed in the work of [START_REF] Canup | Forming a Moon with an Earth-like composition via a giant impact[END_REF] allows the homogenization of the composition before accretion. Nevertheless, the model still struggles to resolve the excess in the angular momentum and also applies a model with simple Ab initio study of silica and hydrated silica during giant impacts composition, where it only considers iron for the core particle and a dunite mantle. Furthermore, recent studies indicate that the probability to conciliate the homogeneity observed in different isotopic systems as an outcome of the standard giant impact scenario (i.e. the canonical model) is of less than 0.5% [START_REF] Fischer | The origin of the moon's earth-like tungsten isotopic composition from dynamical and geochemical modeling[END_REF].

Fast-spinning Earth impact model

The fast-spinning Earth model of lunar formation considers a high energy collision between a small impactor and a fast-spinning proto-Earth. In the model proposed by [START_REF] Cuk | Making the Moon from a fast-spinning Earth: A giant impact followed by resonant despinning[END_REF], the average terrestrial period of rotation is 2.5h at the moment of impact, where posterior orbital resonance between the Sun and the Moon would reduce the angular momentum. In this model, most of disk material comes from the target. The high energy of the impact vaporizes silicates and forms a silica atmosphere and vapor rich disk. Although the mixing could resolve isotopic signatures of the Earth-Moon system, modeling of the chemical evolution is necessary to further investigate final compositions. The density profile of the debris disk varies from 20 g/cm 3 in the center to 0.01 g/cm 3 within the Roche limit and towards 10 -4 g/cm 3 in the outermost parts. [START_REF] Lock | The structure of terrestrial bodies: Impact heating, corotation limits and synestias[END_REF] proposed that an impact in a fast-spinning body would produce a new planetary object, which was named synestia. [START_REF] Lock | The origin of the Moon within a terrestrial synestia[END_REF] describe a Moon forming impact with the formation of a synestia, where moonlets form from condensation and equilibration with a BSE atmosphere in low pressure regions outside the Roche limit, where equilibration temperatures would be close to 4000 K. This model projects enough mass in the disk to account for the large lunar mass with mixing of impactor and target material. The authors consider that partitioning of volatiles in the vapor would explain the observed lunar budget, but the model is not consistent with current estimates of water contents [START_REF] Saal | Volatile content of lunar volcanic glasses and the presence of water in the Moon's interior[END_REF] in the Moon.

Equations of state and phase diagrams

Equations of state (EOS) summarize the thermodynamic properties of a given material, usually expressed as a functional relation connecting the parameters that describe the material in the form of P(p,T) [START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF]. They are an essential input in current smoothed-particle hydrodynamics (SPH) impact simulations because they predict pressure-temperature-velocity variations, ultimately governing the final outcome of any scenario. For instance, the distribution of the mass in the disk is affected not only by the choice of material, but also by the equation of state the describe the material [START_REF] Benz | The origin of the moon and the single-impact hypothesis III[END_REF], showing ultimately an influence in mass injection, elemental distribution and multiphase partitioning.

The EOS initially used in SPH simulations was the one developed by [START_REF] Tillotson | Metallic equations of state for hypervelocity impact[END_REF]. Observations on problematic representation of the thermodynamic behavior of multiphase states showed a significant misrepresentation of the effects of pressure gradients in the gas [START_REF] Benz | The origin of the moon and the single-impact hypothesis III[END_REF], favoring the utilization of a different equation of state, called ANEOS. ANEOS was developed by [START_REF] Thompson | Improvements in the CHART D radiationhydrodynamic code III: revised analytic equations of state[END_REF] and it was based on analytical approximations for application in metallic systems. It gave a better representation of mixed phases, but it treated the gas as only monoatomic particles. This assumption lead to a high-energy standard for vaporization, underestimating of post-impact gas formation. [START_REF] Melosh | A New and Improved Equation of State for Impact Computations[END_REF] corrected the previous assumption that the gas is monoatomic with a new formulation of the ANEOS, allowing the presence of clusters and, therefore, leading to a higher gas proportion in final simulations. In a comparison between the original and the modified version, the author demonstrated an increase of a factor of 3 in the amount of vapor that could be generated from impact [START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF]. However, the description of the gas only accommodates diatomic molecules, which can be problematic for materials presenting a gas phase largely dominated by molecules with more than two atoms, such as water. Furthermore, the improvement in ANEOS made by [START_REF] Melosh | A New and Improved Equation of State for Impact Computations[END_REF] was also reported to impact the amount of the total mass in the disk and to increase its contents of iron [START_REF] Canup | Dynamics of lunar formation[END_REF].

Goal of the thesis

Models of lunar formation rely strongly in equations of state that describe the behavior of materials. The main goal of this work is to provide accurate information concerning key materials involved in the formation of the Moon at conditions relevant to current models available. To do so, we specifically focus on the behavior of silica and silica-water systems at low-density and high-temperature regimes.

SiO 2 is one of the major building blocks of the Earth-Moon system and plays a central role in lunar formation processes. We aim to describe SiO 2 behavior in better detail focusing on the construction of a reliable phase diagram in the liquid-vapor region and in establishing an accurate critical point. Subsequently, we focus on the study of structural and transport properties under varying conditions, comparing the behavior of the main elements at every step. Furthermore, we aim to assess the influence of the presence of water in the silicate system, in terms of positioning of the critical point and variation of phases, where we also focus on the potential effects on structural and transport properties of individual particles. To do so, we construct four hydrated systems and then compare the variations with the anhydrous SiO 2 baseline.

Chapter 2

Methodology Contents

Multiscale modeling

For many years, the grounds for understanding the behavior of a particular system was mostly given by theoretical formulation and experimentation. In order to accurately describe the implications of certain conditions to a material, theory and empirical knowledge were always confronted with each other. With the exponential increase in computing power and numerical methods, modeling became an essential tool to predict behaviors of systems or to design and test potentially interesting new materials. Moreover, it can often complement experimental work by allowing modeling at extreme conditions. Particularly, planet forming scenarios often involve extremely high temperatures and low densities, conditions generally difficult to reach and that generate limited data points in most current experimental setups.

In the past few decades, it became possible to bridge between scales and describe larger systems from finer scales with the introduction of multiscale modeling. As shown schematically in Figure 2.1, with this technique, a system of interest can be modeled by combining methods of different complexities, given that the physics at finer scale is better understood than that at coarser scale [START_REF] Fish | Mesoscopic and multiscale modelling in materials[END_REF]. For instance, macroscale description of plastic and elastic behavior of a given material can be coupled with modeling of microscale dislocations in order to provide information concerning discrete phenomena that needs understanding in better detail, such as chemical reactions or local defects [START_REF] Weinan | Principles of multiscale modeling[END_REF]. By applying multiscale modeling, is also possible to combine models that provide information essentially regarding displacement of particles in space, such as molecular dynamics, with subatomic models that describe the behavior of the electrons, i.e. ab initio methods. The choice of the method depends on the aimed resolution and it must balance accuracy and complexity, since finer-scale calculations can rapidly become too complex to handle even by modern supercomputers.

Ab initio methods describe subatomic scales of a given material from a quantum mechanics view point. It consider electrons and nuclei as the basic particles and depend essentially on the physics of the atomic species present in the system to be simulated [START_REF] Lee | Computational materials science: an introduction[END_REF]. The combination of ab initio methods with molecular dynamics approaches allows the description of the trajectory of the particles while calculating force and energy from first principles. It provides accurate thermodynamic information and also allows the analysis of electronic and magnetic properties. Density functional theory (DFT) is among the most popular approaches of ab initio calculations [START_REF] Van Mourik | Density functional theory across chemistry, physics and biology[END_REF] with the goal to solve Schrödinger equation by recasting the problem and maintaining accuracy. Coupled with molecular dynamics methods (MD), it is a powerful tool to derive equations of state and estimate phase diagrams under various conditions [START_REF] Söderlind | Assessing density-functional theory for equation-of-state[END_REF]. The DFT-MD approach will be detailed in the next section. Further reading concerning multiscale modeling and current available methods can be found in [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF], [START_REF] Fish | Practical Multiscaling[END_REF], [START_REF] Fish | Mesoscopic and multiscale modelling in materials[END_REF] and references therein.

Ab initio molecular dynamics simulations

Until the mid-80s, available molecular dynamics methods relied on empirical interatomic potentials which did not allow access to any information regarding the electronic structure. The first combined approach of molecular dynamics with density functional theory was proposed by [START_REF] Car | Unified approach for molecular dynamics and densityfunctional theory[END_REF], where the electronic wavefunctions are introduced as dynamic variables and the update of the electronic structure is calculated on-the-fly with ionic motion [START_REF] Boero | Car-parrinello molecular dynamics, Encyclopedia of Nanotechnology[END_REF]. By coupling these techniques, it was possible for the first time to compute ground-state electronic properties of larger systems from first principles while having the ionic motions given by classical mechanics. In modern DFT, the Born-Oppenheimer approximation was introduced, where the motion of nuclei and electrons are decoupled and the ground-state energy is calculated at every timestep. In a nutshell, ab initio molecular dynamics (AIMD) aims to describe accurately dynamic systems from a quantum mechanics perspective by solving Schrödinger equation. The means to realistically do so requires a series of approximations that can be computationally demanding (Table 2.1).

Schrödinger equation

In quantum mechanics, the Schrödinger equation describes the time evolution of particles for a given quantum state of a system. For a one-dimensional space, the time-dependent equation can be described as:

ih ∂ ∂t ψ = Ĥψ (2.1)
where the wavefunction ψ contains all the system properties for every particle in n-dimensional space at a given time t. In a system that do not change with time, the Schrödinger timeindependent equation can be generally written as:

Eψ = Ĥψ (2.2)
where E is the energy and Ĥ is the Hamilton operator. The Hamiltonian Ĥ can be expressed as: Ĥ = Te + Tn + Ven + Vnn + Vee (2.3) where Te and Tn are the kinetic energy operator for the electrons and nuclei, respectively, and Ven , Vnn and Vee are the potential energy operators accounting for nucleus-electron, nucleusnucleus and electron-electron given via the Coulomb potential. Due to its high dimensionality with 3N coordinates for N particles, the complexity of this equation allows it only to be solved for atoms with small number of particles, such as hydrogen and helium. Solving Schrödinger equation for complex systems requires a number of approximations in order to address the manyparticle problem. Reformulations implemented in the DFT method applied in the calculations of this work are going to be briefly described in the following sections.

Born-Oppenheimer approximation

In the approach proposed by [START_REF] Born | Zur quantentheorie der molekeln[END_REF], the motion of nuclei and electrons is separated from each other based on the assumption that their large difference in mass results in a difference in motion rates, with the light electrons adjusting to the dynamic of the heavier nuclei. Considering nuclei as fixed in space, the electronic wavefunction, then, depends parametrically on the position of the nucleus and remain at their instantaneous groundstate, providing a potential energy surface in which the nucleus would move. By decoupling the movement of the particles, nuclear and electronic wavefunctions can be expressed as a product and treated separately, where the electronic Hamiltonian is then expressed as: Ĥe = Te + Ven + Vee (2.4) The Born-Oppenheimer approximation helped to develop electronic methods that do not explicitly depend on nuclear motion. With its implementation in ab initio molecular dynamics methods, the energy is minimized and forces are calculated at each DFT cycle, with minor restrictions concerning the length of the timestep [START_REF] Kühne | Second generation Car-Parrinello molecular dynamics[END_REF].

Density Functional Theory

Density functional theory (DFT) is the result of the work from Hohenberg, Kohn and Sham on a series of theorem and equations in the 1960s. [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF] determined in a two-folded theorem that the ground-state energy is a functional of the electronic density, ultimately replacing the wavefunction by the electron density n(r). [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] mapped the many-body interacting problem on an non-interacting manner, where the unknown terms are placed in the exchange-correlation functional and the minimizing density is established in terms of single-particle orbitals. DFT became more predictive in the with the development of more appropriate forms of the exchange-correlation functionals that increase the accuracy of the method. The theorems will be presented briefly herein.

Hohenberg and Kohn theorem

Inspired by the precedent Thomas-Fermi theory [START_REF] Fermi | Un metodo statistico per la determinazione di alcune prioprietà dell'atomo[END_REF][START_REF] Thomas | The calculation of atomic fields[END_REF][START_REF] Kohn | Nobel lecture: Electronic structure of matter-wave functions and density functionals[END_REF] that introduced the main idea of electron density n(r), [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF] applied the concept to the treatment of the many-body electron problem. In general terms, the two-folded theorem demonstrate that the energy is a functional of the electron density and, therefore, all ground-state observables depend only on the electron density. The relation that expresses the ground-state energy functional E HK [n] can be written as:

E HK [n] = F HK [n] + ⟨ Vext ⟩ (2.5)
where ⟨ Vext ⟩ is the external potential. The free energy F HK [n] is a functional dependent on the kinetic T and potential electron-electron interaction Vee operators, which can be express as:

F HK [n] = ⟨ T ⟩ + ⟨ Vee ⟩ (2.6)
The ground-state density is given by energy minimization through a variational approach with trial wavefunctions. In this approximation, the solution depends only on the electronic positions as a function of 3 dimensional coordinates and independently of N.

Kohn-Sham equations and exchange-correlation functional

The Kohn-Sham approach [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] mapped the interacting particle problem on a non-interacting manner, where the minimization of the energy is not given with respect to the many-body wavefunction, but in terms of single-particle orbitals [START_REF] Marx | An introduction to ab initio molecular dynamics simulations[END_REF]. The free energy contains only terms that can be calculated exactly where the remaining unknown quantities are incorporated in the exchange-correlation term. The energy can then be expressed as:

E KS [n] = T ip [n] + ⟨ Vext ⟩ + E H [n] + E xc [n] (2.7)
where T ip is the kinetic energy of independent particles and E H is the Hartree potential energy. The exchange-correlation energy is the relation between the known quantities with the approximated terms, given by:

E xc [n] = ⟨ T ⟩ -T ip [n] + ⟨ Vee ⟩ -E H [n] (2.8)
where ⟨ T ⟩ is the exact kinetic energy and ⟨ Vee ⟩ is the exact interaction energy term. The oneparticle Schrodinger equation can be written as:

- h2 2m + v ef f (r) ϕ i (r) = ϵ i ϕ i (r) (2.9)
where h is Planck's reduced constant, m is mass and the electron density, eigenvectors ϕ i and eigenvalues ϵ i are given by:

n(r) = ∞ n=1 f i |ϕ i (r)| 2 (2.10)
where f i is the integer occupation number of the i th ϕ i orbital. The minimum energy for a fixed number of particles yields the effective potential v ef f [START_REF] Martin | Interacting Electrons[END_REF], as in:

v ef f (r) = v ext (r) + δE H [n] δn(r) + δE xc [n] δn(r) (2.11)
These three equations (2.9 -2.11) enter a self-consistent cycle in order to be solved since the effective potential that yields the orbitals which determine the density is also a functional of the density. Furthermore, the quality of the approximation used for the exchange-correlation energy term sets the quality of the electronic structure calculation [START_REF] Iftimie | Ab initio molecular dynamics: Concepts, recent developments, and future trends[END_REF]. The local-density approximation (LDA) was introduced in the original KS work [START_REF] Martin | Interacting Electrons[END_REF]. Generalized-gradient approximation (GGA) was introduced later on and improved accuracy in structure and bonding estimations in comparison with LDA [START_REF] Martin | Interacting Electrons[END_REF]. Different exchange-correlation functionals have been introduced over the years and extensive research has been done comparing the effects of the functionals over properties of materials [START_REF] Fritsch | Climbing Jacob's ladder: A density functional theory case study for Ag 2 ZnSnSe 4 and Cu 2 ZnSnSe 4[END_REF]. Further explanation of the Kohn-Sham approximation can be found in [START_REF] Martin | Interacting Electrons[END_REF] and [START_REF] Woods | Computing the self-consistent field in Kohn-Sham density functional theory[END_REF].

Molecular dynamics

The molecular dynamics method allows the computation of equilibrium and transport properties of a many-body system using classical mechanics solving Newton's law of motion [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF]. In this method, atoms are treated as hard spheres with established initial positions and move by the integration of net forces over incremental timesteps that can vary roughly from femtoseconds (fs) to picoseconds (ps). In classical MD, the forces acting on atoms are calculated from interatomic potentials, whereas, coupled with DFT, forces are calculated from first principles, as explained previously in section 2.2.3.

In a MD cycle, every atom starts from unique initial positions and velocities. The interatomic forces acting on each atom are calculated based on the method of choice and the net value is then computed. From the net forces, the velocities are estimated from Newton's law of motion, based on the assumption that they are constant over short time increments. The velocities are ccontrolled with the help of a thermostat in order to assure the energy of the system remain stable. Once the new forces and velocities are established, the positions of the atoms are updated through integration algorithms (e.g. Verlet or Leap-Frog). The MD cycle restarts systematically until an appropriate preestablished number of iterations is achieved, where equilibrium is reached and the macroscopic properties can be calculated. Thermostat, integration algorithms and timestep balance the accuracy of the method, where a comprise with the speed of calculations should be made.

The constraints imposed to the collection of iterations from the MD cycle represents the statistical ensemble on which the simulation is performed [START_REF] Mazurek | A review on combination of ab initio molecular dynamics and NMR parameters calculations[END_REF]. These constraints can be chosen based on the properties of interest and typically concern energy (E), pressure (P), temperature (T) and number of particles (N). For instance, the canonical ensemble fixes the number of particles, the volume and the temperature and is, therefore, abbreviated as NVT. In such ensemble, the temperature is regulated with a thermostat, where it is allowed to oscillate around the nominal value. Examples of different possible ensembles are also the microcanonical (NVE) and the isobaric-isothermic (NPT). Further explanation can be found in [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF], [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF], [START_REF] Lee | Computational materials science: an introduction[END_REF] and references therein.

In this work, we employ the DFT-MD method as implemented in the Vienna Ab initio Simulation Package (VASP ® ; [START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF], which is schematically represented in Figure 2.2. The simulation cycle starts from the initial position of the ions, where the electron density is initially guessed and the free energy is calculated based on the single-particle orbital approximations derived by Kohn-Sham. The DFT loop is repeated until the energy is minimized to be then plugged-in into Newton's equation of motion for the nuclei, where the positions are updated. The properties of the system are computed at the end of each cycle.

Initial configurations

The initial positions of atoms used as input in our AIMD simulations were obtained from the cristobalite mineral structure from [START_REF] Downs | The pressure behavior of α cristobalite[END_REF], available in the American Mineralogist Crystal Structure Database. Cristobalite is one of the low-pressure and hightemperature polymorphs of SiO 2 , typically stable up to 2000 K, and with a tetragonal structure. Its unit cell contains 12 atoms with lattice parameters established at a = 4.97170, b = 4.97170 and c = 6.92230 Å. At this configuration, SiO 2 has 16 electrons and 2.33 g/cm 3 at ambient conditions. In order to obtain a better representativeness of the system, we build a supercell by placing additional unit cells adjacently to each other, observing that the final shape approaches a cubic structure to maintain the isotropy of forces acting on the system. The unit cell was multiplied in a 3x3x2 proportion along the aforementioned axis and resulted in a supercell with a = 14.91510, b = 14.91510 and c = 13.84460 Å as lattice parameters, summing 216 atoms and 1152 electrons. The new structure was, then, forced into a cubic structure with 15 Å of lattice (Figure 2.3) and served as the material for the subsequent pure SiO 2 simulations and as basis for the construction of the silica-water mixtures. The modification of the structure represents a difference of less than 10% in volume from the original supercell and does not have a significant impact on final results given that this study focus on equilibrated fluid phases.

The unit cell of cristobalite was used as basis for the construction of silica-water mixtures, where H 2 O units were manually placed in the intermolecular space of the silicate chain, in order to have a evenly distribution of molecular types. We placed one, two and four H 2 O molecules inside the unit cell and then continued to build the supercell in the same manner as previously done for pure SiO 2 , resulting in three mixtures with silica-water proportions of 4:1, 2:1 and 1:1. The fourth system was constructed by removing alternately half of the H 2 O molecules of the system, forming a 8:1 silica-water proportion. In terms of weight percent of H 2 O, these systems vary from 3.61to 23.07 wt%. The respective concentrations, together with the total number of particles and electrons of each system is indicated in Table 2.2.

Computational details

This work was carried using ab initio molecular dynamics methods (AIMD) applying densityfunctional theory (DFT) implemented in the Vienna Ab Initio Simulation Package (VASP; [START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF]. The exchange-correlation functional was treated based on Perdew-Burke-Erzenhof (PBE; [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] generalized gradient approximation (GGA), which considers the electron density and its gradient at a given point [START_REF] Lee | Computational materials science: an introduction[END_REF]. The wavefunctions are represented in the project-augmented wave (PAW; [START_REF] Blöchl | Projector augmented-wave method[END_REF] pseudopotentials, where the core electrons are frozen and only the valence electrons are explicitly considered in the calculations, with a cutoff energy to determine how many planewaves are used. The cutoff energy was set at 550 eV and the augmentation sphere at 800 eV, and the Brouillon zone was sampled at the gamma-point. Convergence tests can be found in the Appendix A.

The atoms were placed inside a cubic cell with periodic boundary conditions, according to the method described in the previous section. The simulations were performed in the canonical ensemble (NVT), where the temperature was controlled by the Nosé-Hoover thermostat [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Hoover | Canonical dynamics: Equilibrium phase-space distributions[END_REF]. The choice of density and temperature was based on available phase diagrams for silicate systems [START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF] along with predicted protolunar disk conditions (see Section 1), aiming specifically for the most likely region where the critical point is placed.

Simulations on anhydrous SiO 2 were performed on a cell containing 72 formula units, summing 216 atoms and 1152 electrons. The density varied from 0.2 to 2.33 g/cm 3 and temperature values were set at 4000, 4500, 5000, 5500, 6000 and 7000 K (Figure 2.5). Thermodynamic properties were calculated over time intervals of 1 fs and averaged over the complete time length, which varied from 10 to 70 ps depending on initial equilibration. Simulations on hydrated systems were performed on four different compositions, where the number of atoms and electrons of each individual system can be found on Table 2.2. Densities varied from 0.34 to 2.77 g/cm 3 and the temperature values were set to 2000, 3000, 4000 and 5000K, depending on the system (Figure 2.6). The time interval was set to 0.5 fs and the thermodynamic properties were averaged over 9 to 20 ps length, also depending on initial equilibration.

A total of 84 points were calculated for pure SiO 2 and 137 points for the hydrated systems dataset, summing a total of 221 computed systems. The results of the production runs were processed using the UMD package [START_REF] Caracas | Analyzing melts and fluids from ab initio molecular dynamics simulations with the UMD package[END_REF], where initial equilibration steps were discarded.

Equations of state and phase diagrams

The construction of the thermal equations of state P(ρ,T) is done by the compilation of the DFT-MD simulations dataset. Total pressure and internal energy are given by the sum of three quantities related to the ionic, electronic and interacting term contributions. The pressure at each density point is averaged from the loop of iterations and then computed to construct the isotherms. The ensemble of P(ρ,T) conditions defines the equations of state of the studied systems, on which we fit a third-order polynomial equation in order to approximate the behavior of phases that serve as basis for the construction of phase diagrams.

The first-derivative of the cubic polynomial fit gives the local minima and maxima of the respective isotherms. These local values aforementioned are defined as spinodal points, meaning that they delimit the thermodynamic stability of a phase-homogeneous fluid. The curve delimited by those points comprises the liquid-vapor coexistence dome and the metastability domain. To the left of the local maxima is the gas domain, whereas to the right of the local minima is the liquid domain. In the supercritical regime, liquid and gas behave as a continuous phase and, therefore, spinodals become undefined. The schematic phase-diagram as described above is illustrated in Figure 2.7

In the van der Waals equation of state, the critical constants relate to the van der Waals coefficients and the critical point is established at a flat inflection where the first and secondorder derivatives are zero [START_REF] Atkins | Atkins' Physical Chemistry[END_REF]. In order to estimate the critical constants of our systems, we use the linear regression of the second-order derivatives bracketed by the spinodal points to constraint the critical point pressure-density-temperature setting.

Structural properties

Atomic pair distribution function and coordination numbers

Fluid structure is derived from the atomic ordering analysis given by pair distribution functions (PDF). They describe quantitatively the interatomic distances for a given atomic pair from a central atom at a given radial distance to its neighbors, as shown in Figure E.1a and expressed by the equation below:

g ij (r) = V N (N -1) ⟨ i j̸ =i δ(r -r ij )⟩ (2.12)
where r is the radial distance between particles i and j in a system with N atoms with volume V , δ is the Dirac delta function and ⟨⟩ is the ensemble average. The first coordination sphere of a given particle is defined by the nearest neighboring atoms placed in a shell at distance r0 ( The structure of the atomic pair reflects directly the general structure of the fluid. The sharpness or broadness of the peaks indicate the positioning of the atoms in space and, therefore, it becomes possible to estimate the phase state of the system. PDF also provide information regarding the different molecular species present in the fluid, where the integral of the pair distribution function up to the first minimum yields the coordination number of the atomic pair (Figure E.1b), as in:

C (rm) = rm 0 g(r)4πN r 2 dr (2.13)
where C is the coordination number and r m is the radial distance at the first minimum trough of the PDF.

Polymerization analysis

Silicate systems are known for forming polymers of variable length depending on thermodynamic conditions. The degree of connectivity of the chain on such systems is associated with fluid properties, such as viscosity, self-diffusion coefficients and more [START_REF] De Clermont Gallerande | Quantification of non-bridging oxygens in silicates using x-ray raman scattering[END_REF], which are key to understand the behavior of silicate systems at larger scales. The analysis of the variation on coordination numbers proportions can indicate the extent of molecular dissociation in the fluid. Considering this, we propose a new measurement for estimating the degree of polymerization based on the evaluation of the relation between coordinations of silicon and oxygen atoms with respect to each other. The relation is given by the weighted sum of the coordination of silicon with respect to oxygen over the weighted sum of the coordination of oxygen with respect to silicon, expressed as:

V T,ρ = [SiO x ]x [OSi y ]y (2.14)
where the sum of the fractions of all Si-O coordination numbers x weighted by their index x over the fractions of all O-Si coordination numbers y weighted by their index y yields the polymerization degree V at a given temperature T and density ρ. The general evaluation of this relation allows the assessment of the Si-O-Si bridging and describes the structure of the fluid by providing the amount of overall bonding between the two atoms.

Transport properties

The motion of particles is measured directly from the calculated trajectories in terms of the mean-square displacement (MSD) and self-diffusion coefficient (D). The MSD scales linearly with time in the fluid and is given by the deviation of a initial position in respect with time, as described in the equation:

M SD = ⟨[r i (t 0 + t) -r i t 0 ] 2 ⟩ (2.15)
where r is the position of the particle i from a initial time t 0 to time t, over the ensemble average ⟨⟩. The slope of the MSD yields the self-diffusion coefficient according to the relation below:

D = 1 6 lim t→∞ M SD t t (2.16)
The diffusion coefficient is calculated for every simulated point and fitted with a linear function in order to describe the dynamic behavior of the atoms at fixed temperature or fixed density.

Phase diagram and the critical point

From our ab initio molecular dynamics calculations performed on anhydrous SiO 2 systems, we build our equations of state and present the results in the pressure-density phase diagram shown in Figure 3.1. A total of six temperatures were simulated, on a range varying from 4000 to 7000 K. For each temperature-density condition, the value of the pressure was computed and is represented by an individual point in the aforementioned diagram. A table with all discrete values of pressure is given in Appendix B.1. We started from the normal density of cristobalite (2.33 g/cm 3 ) at all isotherms and we sample the phase space by gradually increasing the cell volume in order to reach lower densities, down to 0.20 g/cm 3 . Reported data on SiO 2 [START_REF] Green | Bulk properties and near-critical behaviour of SiO 2 fluid[END_REF][START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF] showed that previous predictions of the critical point for this material should be generally between 5000 and 6000 K, with a most likely value approaching 5400 K. Considering this, we space isotherms every 500 K up to 6000 K and we increase the number of simulation points in order to obtain a better resolution in these conditions.

The dataset is fitted with a cubic function up to 2 g/cm 3 for each isotherm in order to represent the mean behavior of the system. The minima and maxima values are extracted by calculating the first derivatives of the fitted curve and represent the spinodal points, as explained in detail in section 2.5. The region bracketed by the spinodal point along with the linear regression of the second order derivatives indicates the most likely region for the critical point. We place the critical point of SiO 2 between 0.6 and 0.85 g/cm 3 , 5000 to 5500 K and 0.15 to 0.25 GPa. Above 5500 K, the spinodal points are indistinct and, therefore, we define SiO 2 in supercritical state. We define the liquid-vapor coexistence dome between the liquid and gas spinodals, defined by the gray area in Figure 3.1. The binodal curve defining the metastability region of the fluid is not represented in this phase diagram.

Critical points of rock-forming silicate minerals have been vastly documented [START_REF] Kobsch | Behavior of feldspars during the Giant Impact[END_REF][START_REF] Kobsch | The critical point and the supercritical state of alkali feldspars: Implications for the behavior of the crust during impacts[END_REF][START_REF] Xiao | Critical vaporization of MgSiO 3[END_REF][START_REF] Townsend | Liquid-vapor coexistence and critical point of Mg 2 SiO 4 from ab initio simulations[END_REF][START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF][START_REF] Green | Bulk properties and near-critical behaviour of SiO 2 fluid[END_REF]. We show in Figure 3.2 the temperature-density critical constraints reported in the literature on the different minerals, along with the critical point of SiO 2 as calculated in this study. Minerals containing refractory elements such as calcium, aluminum and magnesium, report critical points at higher temperatures, where anorthite (CaAl 2 Si 2 O 8 ) is considerably higher in comparison to the other minerals, placed in between 7000 and 7500 K [START_REF] Kobsch | Behavior of feldspars during the Giant Impact[END_REF], follow by forsterite between 6550 -6750 K (MgSiO 3 ; Xiao and Stixrude, 2018) and enstatite 6040 -6440 K (Mg 2 SiO 4 ; [START_REF] Townsend | Liquid-vapor coexistence and critical point of Mg 2 SiO 4 from ab initio simulations[END_REF]. Alkali feldspars form a solid solution with interchangeable cations and behave similarly with respect to their near [START_REF] Kobsch | Behavior of feldspars during the Giant Impact[END_REF] SiO 2 [START_REF] Green | Bulk properties and near-critical behaviour of SiO 2 fluid[END_REF] MgSiO 3 [START_REF] Xiao | Critical vaporization of MgSiO 3[END_REF] Mg 2 SiO 4 [START_REF] Townsend | Liquid-vapor coexistence and critical point of Mg 2 SiO 4 from ab initio simulations[END_REF] SiO 2 [START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF] NaAlSi 3 O 8 [START_REF] Kobsch | The critical point and the supercritical state of alkali feldspars: Implications for the behavior of the crust during impacts[END_REF] KAlSi 3 O 8 [START_REF] Kobsch | The critical point and the supercritical state of alkali feldspars: Implications for the behavior of the crust during impacts[END_REF] [START_REF] Kobsch | The critical point and the supercritical state of alkali feldspars: Implications for the behavior of the crust during impacts[END_REF]; SiO 2 [START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF][START_REF] Green | Bulk properties and near-critical behaviour of SiO 2 fluid[END_REF]; forsterite, MgSiO 3 [START_REF] Xiao | Critical vaporization of MgSiO 3[END_REF]; and enstatite, Mg 2 SiO 4 [START_REF] Townsend | Liquid-vapor coexistence and critical point of Mg 2 SiO 4 from ab initio simulations[END_REF]. critical characteristics, where their critical temperature approach that of SiO 2 . The critical point of SiO 2 as calculated in this work shows a good agreement with previous reported values that placed it at 5400 K, 0.55 g/cm 3 and 0.19 GPa [START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF] and, more recently, between 5000-6000 K, 0.5 -1.00 g/cm 3 and 0.2 -0.5 GPa [START_REF] Green | Bulk properties and near-critical behaviour of SiO 2 fluid[END_REF].

The vast majority of works presented in Figure 3.2 estimate critical points by applying ab initio molecular dynamic techniques [START_REF] Kobsch | Behavior of feldspars during the Giant Impact[END_REF][START_REF] Kobsch | The critical point and the supercritical state of alkali feldspars: Implications for the behavior of the crust during impacts[END_REF][START_REF] Xiao | Critical vaporization of MgSiO 3[END_REF][START_REF] Townsend | Liquid-vapor coexistence and critical point of Mg 2 SiO 4 from ab initio simulations[END_REF][START_REF] Green | Bulk properties and near-critical behaviour of SiO 2 fluid[END_REF], with differences in the use of such techniques by varying parameters and statistical ensembles. For instance, the energy cutoff was set at 500 eV by [START_REF] Xiao | Critical vaporization of MgSiO 3[END_REF] and 800 eV by [START_REF] Townsend | Liquid-vapor coexistence and critical point of Mg 2 SiO 4 from ab initio simulations[END_REF]. Higher accuracy is generally expected when using higher values for energy cutoff, however, simulations take oftentimes longer to calculate. Observing that the system has properly converged at a lower cutoff, results will be precise and calculated in a shorter period of time. The critical point of silica published by [START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF] is estimated from an analytical equation of state, where its construction is discussed in Chapter 1. The use of fundamentally different techniques towards the same goal is mutually beneficial, since they can complement and corroborate each other. Recent work has also shown benefits in adding an experimental bias onto numerical simulations [START_REF] White | Communication: Improved ab initio molecular dynamics by minimally biasing with experimental data[END_REF].

Structure of the fluid

Pair distribution functions

We describe the structure of the fluid starting from a pair distribution function analysis. We compute the PDF for Si-O and O-O pairs at all densities (Figure 3.3) and temperatures (Figure 3.4). We present a detailed table containing positions of maxima and minima on the complete span of pair distribution functions curves in the appendix C. The nearest neighbors are represented in the first coordination sphere, where the first peak represents the average bond length of the represented pair. The second coordination sphere accounts for the general polymerization of the fluid. As mentioned in section 2.6, the structure of the PDF reflects directly the structure of the fluid. The beginning of the PDF reflects the repulsive forces acting between the atoms that prevent them to overlap, which can differ up to 15% as a function of variations in temperature-density of the specified pair.

Figure 3.3 shows a comparison of pair distribution functions in the fluid below (4000 K) and above (7000 K) critical conditions for Si-O and O-O pairs spanning all calculated densities. Panels on the left represent the function at 4000 K and panels on the right the function at 7000 K. The upper part stands for Si-O, whereas the lower part stands for O-O. At all conditions, we observe a general broadening of the peaks as the density decreases. This width increase suggests a looser packing of the structure, which is also consistent with dissociation of molecules related with increasing vaporization of the fluid. With the decrease in density, we also generally observe the placement of the first minimum well above zero, indicating an increasing dynamic behavior of the particles [START_REF] Green | Bulk properties and near-critical behaviour of SiO 2 fluid[END_REF]. Such behavior becomes remarkably significant for O-O pairs at 4000 K in densities equal and below 0.27 g/cm 3 , where the fluid would be well inside the gas phase region.

At 4000 K, the fluid is below critical and mostly inside the liquid-vapor coexistence dome depending on the density (see Figure 3.1). The position of the first peak of the Si-O pair varies from 1.64 Å at 2.33 g/cm 3 to 1.63 Å at 0.20 g/cm 3 , with a mean value of 1.63 Å. First minima values vary from 2.41 Å at 2.33 g/cm 3 to 2.51 Å at 0.20 g/cm 3 , with a mean value of 2.49 Å. However, at 7000 K where the fluid is at supercritical conditions, we observe the position of the first peak varying from 1.64 Å at 2.33 g/cm 3 to 1.57 Å at 0.27 g/cm 3 , with a mean value of 1.61 Å, and first minima values varying from 2.56 Å at 2.33 g/cm 3 to 2.84 Å at 0.27 g/cm 3 , with an average of 2.69 Å. We observe that the second coordination sphere is more prominent at 4000 K, suggesting a higher level of organization and degree of polymerization of the fluid in comparison with 7000 K.

Considering the O-O pair, at 4000 K the main peak is located between 2 and 3 Å and is related to the positioning of oxygens within the tetrahedron. It varies from 2.68 Å at 2.33 g/cm 3 to 2.72 Å at 0.20 g/cm 3 , averaging over 2.72 Å. We observe an incipient peak between 1 a 2 Å at densities equal to 2.00 g/cm 3 or lower. It becomes more significant at lower densities and it represents the formation of diatomic oxygen molecules. At 7000 K, tetrahedra oxygen peak shifts slightly to the right, varying from 2.67 Å at 2.33 g/cm 3 to 2.86 Å at 0.27 g/cm 3 , with an average of 2.80 Å. The molecular oxygen peak is present at all densities and becomes dominant at densities equal to 0.78 g/cm 3 or lower.

By fixing the value of the density, we are able to compare and evaluate the effects of temperature onto the organization of the structure. We study the structural behavior variation at 2.33 g/cm 3 and 0.27 g/cm 3 for every isotherm, as shows in Figure 3.4. With increasing temperature, the peaks also become wider and the second coordination sphere decreases. Our mean values of bond length of the Si-O pair vary from 1.63 Å at 4000 K to 1.62 Å at 4500 and 5000 K to 1.61 Å at 5500 K and above temperatures. We observe a dependency of molecular oxygen formation with the temperature, represented by the emerging of the characteristic peak that increases in height with increasing temperature in both compared densities. This peak becomes prevalent at temperatures equal or superior to 5500 K, meaning that they become dominant when the fluid is in the supercritical domain.

For both densities, the first trough is placed above zero and the second coordination sphere decreases with increasing temperature. At 2.33 g/cm 3 , the mean distance of the Si-O has a small variation with temperature, approximated at 1.64 Å at all temperatures. For the O-O tetrahedron pair, the distance varies between 2.67 and 2.68 Å with an unspecific relation with temperature. We also observe the effect of thermal expansion in pair distribution functions at 2.33 g/cm 3 , marked by the widening of the peaks. However, this effect is less evident at low density, where we mostly observe a general asymmetry and shift of the main peaks. At 0.27 g/cm 3 , Si-O bonds vary from 1.63 Å at 4000 K to 1.57 Å at 7000 K, with an average value of 1.59 Å.

Our results are consistent with reported data in the literature. Typical values of Si-O bond Left panels show the variation of the distribution at 2.33 g/cm 3 , whereas right panels show the variation of the distribution at 0.27 g/cm 3 . The total radial distance is giving in half of the size of the simulation box. Therefore, the horizontal axis for the 0.27 g/cm 3 diagrams go up to 15 Å given that the lattice for SiO 2 at such conditions is 30 Å, whereas for the diagram with density of 2.33 g/cm 3 the radial distance goes up to 7 since the lattice for those simulations were of 14.554 Å.

length can vary from 1.564 up to 1.720 Å depending on materials, thermodynamic conditions and methods applied, as reported by numerous researchers [START_REF] Baur | The prediction of bond length variations in Silicon-Oxygen bonds[END_REF][START_REF] Gibbs | Correlations between Si-O bond length, Si-O-Si angle and bond overlap populations calculated using extended Hükel molecular orbital theory[END_REF][START_REF] Baur | Silicon-Oxygen bond lengths, bridging angles Si-O-Si and synthetic low tridymite[END_REF][START_REF] Dove | Direct measurement of the si-o bond length and orientational disorder in the high-temperature phase of cristobalite[END_REF][START_REF] Green | Bulk properties and near-critical behaviour of SiO 2 fluid[END_REF] and by the Computational Chemistry Comparison and Benchmark Database of the National Institute of Standards and Technology of the United States. The average value, however, is placed around 1.62 Å. The tetrahedral oxygen distancing is typically close to 2.63 Å [START_REF] Dove | Direct measurement of the si-o bond length and orientational disorder in the high-temperature phase of cristobalite[END_REF], whereas diatomic oxygen bond length is reported varying from 1.21 to 1.36 Å [START_REF] Zaichenko | Pathways to triplet or singlet oxygen during the dissociation of alkali metal superoxides: Insights by multireference calculations of molecular model systems[END_REF]. Peak shift and asymmetry have been reported in pair distribution functions and could be related to an increase in anharmonicity of the potential with increasing temperatures and the shape of typical inter-atomic interaction potentials [START_REF] Ding | Temperature effects on atomic pair distribution functions of melts[END_REF].

The decrease in the average expected bond length of the Si-O pair with respect to the decrease in density is evidenced by the shift to left of the main peaks. This behavior can be explained by the increasing dissociation of the fluid at lower densities, meaning that not all of the Si-O pair is inside the tetrahedron and is increasingly in the molecular form. The same mechanism can explain the variation in bond values observed as a function of temperature at Ab initio study of silica and hydrated silica during giant impacts 0.27 g/cm 3 . The density-temperature Si-O bond length variation is shown as a function of energy in Figure 3.5.

-7 -6.5 -6 -5. 2.33 g/cm 3 2.00 g/cm 3

1.00 g/cm 3 0.52 g/cm 3 0.27 g/cm 3 The shift to the right observed in the oxygen tetrahedral distancing peaks can be related to the thermal expansion [START_REF] Ding | Temperature effects on atomic pair distribution functions of melts[END_REF]. Moreover, molecular oxygen bond length was reported to decrease in a nearly direct relation with the increasing of stretching frequency in some compounds [START_REF] Cramer | Variable character of O-O and M-O bonding in side-on (η 2 ) 1:1 metal complexes of O 2[END_REF], however it can also be related to the formation of differently charged diatomic oxygen species.

Coordination numbers

Crystallographic principles can be applied to melts [START_REF] Mysen | The structure of silicate melts[END_REF]. It is known that the coordination of silica in solids can vary directly with pressure, notably from 4-folded at ambient pressures to 6-folded Si-O at high pressures. In the fluid domain, the organization of SiO 2 can also adapt to the conditions presented. We calculate the fraction of all Si-O species in the fluid starting from the pair distribution function and present it for each density throughout the simulated time, as shown in Figure 3.6. We once more compare the results at conditions below and above critical, at 4000 and 7000 K, respectively.

At 4000 K, there is a major dominance of 4-folded silica at all densities and the fluid is still well polymerized. Nevertheless, as we move towards the gas domain at lower densities, coordination numbers gradually decrease with the formation of smaller Si-O clusters. As mentioned previously, at 4000 K the fluid is inside the liquid-vapor coexistence dome. At 0.20 g/cm 3 , the proportion of 3-folded silica becomes important with a significant increase in coordinations 1 and 2 as well. This behavior is consistent with previous observation on the shift of the Si-O peak in the pair distribution functions described above.

At 7000 K, the dominance of 4-folded silica disappears and the coordination number becomes highly dependent on the density. Lower densities have an increase in lower coordinations and decrease in proportions of coordinations 5 and 6. The average CN transitions from 4 at 2.33 g/cm 3 , to 3 at 0.78 g/cm 3 and, finally to 2-1 at 0.27 g/cm 3 . The 3-folded silica gradually Ab initio study of silica and hydrated silica during giant impacts increases at densities from 2.33 to 1.23 g/cm 3 and starts to decrease from 1.00 to 0.27 g/cm 3 , when the abundance of CN 1 and CN 2 is more significant. This behavior suggests that the degree of dissociation in the supercritical regime relies on the density of the fluid.

We further analyze behavioral variations of the system with respect to temperature at fixed densities (Figure 3.7). On the left panel, we describe CN proportions at 2.33 g/cm 3 and on the right panel, we describe those at 0.27 g/cm 3 . At normal density (2.33 g/cm 3 ), silica is primarily in tetrahedral configuration, evidenced by the dominance of peaks at CN 4. There is a direct relation between temperature and the proportion of species in the fluid. The dissociation is favored with the increase in temperature, where 2 and 3-folded silica is formed. However, the mean pressure varies from 2.71 GPa at 4000 K to 8.47 GPa at 7000 K, which also account for the formation of 5 and 6-folded silica. Conversely, CN 4 decreases in proportion with the escalation of temperature, but remains significant.

At 0.27 g/cm 3 , the main coordination number varies directly with temperature. It shifts from CN 4 at 4000 K to CN 1 at 7000 K in a much more even distribution. Higher packing of silica gradually decreases in favor of the lesser packed. The silica is no longer dominantly in tetrahedra and is much more depolymerized. Coordination number 4000 K 7000 K Figure 3.6: Proportion of coordination numbers per density at fixed temperature. On the left, the panel shows proportions from 0.20 to 2.33 g/cm 3 at 4000 K. On the right, the panel shows proportions from 0.27 g/cm 3 to 2.33 g/cm 3 at 7000 K.

Polymerization

Silicate materials are known for the long polymerization chains. The extent of the polymerization alters the viscosity, which ultimately can affect the distribution of volatiles in the protolunar disk [START_REF] Charnoz | Evolution of the protolunar disk: Dynamics, cooling timescale and implantation of volatiles onto the earth[END_REF]. We establish the polymerization degree as a relation of the total weighted sum ratio (V total ) of the coordination of silicon with respect to oxygen (represented as Si-O x ) and that of oxygen with respect to silicon (represented as O-Si y ), as shown in Figure 3.8 and better detailed in section 2.6.2. Pair distribution functions are restricted to the positioning of atoms in space rather than an analysis of bonding electrons between atoms. Therefore, we have to interpret carefully the results given by such method. 2.33 g/cm 3 0.27 g/cm 3 Coordination Number Left panel shows the variation in proportion of coordination numbers at a density of 2.33 g/cm 3 from 4000 to 7000 K, whereas the right panel show the variation at a density of 0.27 g/cm 3 at the equal temperatures.

However, we observe a strong correlation between thermodynamic properties and the structure of the fluid. As we decrease density, we exit the liquid domain and enter the domain of the gas phase. It becomes increasingly dominated by smaller molecular clusters, which in comparison would have expectedly a lower viscosity and diffuse more rapidly. In Figure 3.8, we present the polymerization ratio as a function of density for every isotherm. The dashed line sets the critical limit, where isotherms below are in supercritical state. Above this line, we explicitly separate the domains of the main phases where the points are located with respect to the positioning inside the phase diagram built and described in section 3.1.

The V total ratio overall varies from 1.498 (Figure 3.9b) to 2.000 (Figure 3.9a). A table containing the weighted sum of x values in Si-Ox, of y values in O-Si y and the final total ratio V total per temperature and density is given in Appendix D. At 2.33 g/cm 3 the fluid is well polymerized at 4000 K, where V is equal to 2.000. With the increase in the temperature, the network gradually starts to break and the value of V decreases to 1.965. This deviation pattern between isotherms concerning V values is repeated throughout the full extent of simulated densities, due to the steeper variation observed at higher temperatures. The deviation is significantly large at 0.27 g/cm 3 , where the gap between polymerization degrees at 4000 and 7000 K is of more than 25%. This behavior indicates that the temperature has a stronger effect at low densities and the break of the network happens more rapidly at such conditions.

Inside the liquid domain, we observe a steady linear behavior directly related to the densitytemperature variations. In the gas domain, the behavior is more erratic and less predictable in the sense that no linearity is observed. However, the general descending trend of polymerization is maintained. It is worth to note that metastable phases are not defined in this diagram, but they should be considered inside the individual gas and liquid domains. In supercritical state, the structure of the fluid also depends on the thermodynamic properties and dissociation rate continues to vary accordingly. Gray areas show the gas domain region at low densities and the liquid domain region at higher densities. In the middle of these two regions is the coexistence area. The dashed line marks the limit of the supercritical regime, where all points below this line are above critical conditions.

In a different analysis, we compare the V total ratio as a function of pressure (Figure 3.10). Here, black lines are showed for reference. They mark the maximum value of V total at 2.000 and the minimum value of pressure at zero. It becomes evident the relation of the network stability with pressure, specially at temperatures above critical (i.e. 5500, 6000 and 7000 K). Isotherms below critical conditions show steady V values until they reach zero pressure, when they vertically drop in a stronger dependency with density. Nevertheless, supercritical systems show a smoother variation of the polymerization degree as a function of pressure, indicating a stronger effect of pressure in the structure of the fluid at the specified temperatures.

Ab initio study of silica and hydrated silica during giant impacts

Dynamic behavior of atoms

Diffusion is a key property and highly dependent on the degree of polymerization of the fluid and viscosity. The self-diffusion coefficients of silicon and oxygen are calculated based on the transport rate obtained from their mean square displacement slope. The mean square displacements are given in Appendix E. We compare self-diffusion coefficients at different temperatures and densities, on the full scope of simulated conditions of the pure SiO 2 system. Displacements of oxygen are generally slightly larger than that of silicon at similar conditions Ab initio study of silica and hydrated silica during giant impacts due to their difference in atomic mass.

Figure 3.11 shows the self-diffusion coefficients as a function of density of Si and O. At fixed temperature, we observe a decrease in the diffusion with increasing density. Generally, the decrease is steeper at lower densities and less expressive at higher densities. The rate of displacement in a direct relation with the temperature, incrementally escalating as temperature increases. In another perspective, we evaluate the variation of the self-diffusion coefficients as a function of temperature at fixed densities (Figure 3.12). We do linear regression of the distribution of simulated points to describe the behavior of the fluid with respect to the thermal variations. Silicon and oxygen atoms diffuse more rapidly at lower densities in a direct linear relation with temperature. Left side panel shows coefficients for silicon at 2.33 g/cm 3 (brown), 1.00 g/cm 3 (pink) and 0.27 g/cm 3 (green). Right side panel shows coefficients for oxygens at same densities.

We describe the characteristic diffusion distance of silicon and oxygen atoms from the respective self-diffusion coefficients at each pressure-density-temperature condition by applying the equation below:

x = √ 2D * t (3.1)
where x is the distance given in meters, D is the self-diffusion coefficient and t is time. We evaluate transport of atoms through time at four different frames: 100yrs, 100 kyr, 100 Myr and 1Gyr. We present results obtained from solving the characteristic diffusion distance equation as a function of self-diffusion coefficient, as shown in Figure 3.13. The comparison of four timeintervals concerns the potential magnitude of time scales implicated in the formation of the Moon and further rocky bodies of the Solar System for comparison purposes. Different lunar formation models consider different duration of accretion, retraction of the disk and posterior equilibration, where these periods could vary from a day [START_REF] Canup | Forming a Moon with an Earth-like composition via a giant impact[END_REF], tens of years [START_REF] Lock | The origin of the Moon within a terrestrial synestia[END_REF] to thousands of years [START_REF] Pahlevan | Equilibration in the aftermath of the lunar-forming giant impact[END_REF]. [START_REF] Charnoz | Evolution of the protolunar disk: Dynamics, cooling timescale and implantation of volatiles onto the earth[END_REF] proposed that the complete accretion of the Moon happened in a period between 10 to 10 4 years, where the temperatures varied between 5000 and 10000 K in the first 10 years. Duration of lunar magma ocean (LMO) with the persistence of elevated temperatures could also vary from 1 to several millions of years [START_REF] Lock | Geochemical constraints on the origin of the moon and preservation of ancient terrestrial heterogeneities[END_REF] and magma oceans throughout the Solar System are reported to persist up to 100 Myr (Elkins-Tanton, 2012). Accretion of planets, core formation of planetary bodies and formation of primary crusts in the Solar system also would follow the latter timescale [START_REF] Chambers | Making more terrestrial planets[END_REF][START_REF] Solomatov | Magma oceans and primordial mantle differentiation[END_REF][START_REF] Taylor | Planetary crusts: their composition, origin and evolution[END_REF]. However, the condensation of first solids would be in the order of a few million years and posterior processes involved in the geological evolution of the Earth could reach up to 1 Gyr.

The distance travelled by atoms depend on their self-diffusion coefficient, which is directly associated with the density of the system, where particles within a denser state tend to travel slower than particles on a more decompress system. As a matter of comparison, we point out the variation of distances with densities of 2.33 g/cm 3 and 0.27 g/cm 3 and we present summarized results in Table 3.1. In a 100 yr interval, when accretion and retraction of the disk in most lunar forming scenarios would be completed, silicon particles could reach between 4.62 and 35.8 m at 4000 K and between 20.6 and 69.5 m at 7000 K. However, transport of oxygen atoms is faster than silicon, and these particles could travel distances from 5.36 to 39.9 m at 4000 K and from 22.6 to 74.3 m at 7000 K.

Introduction

We studied four hydrated systems where the concentration of water increases from 3.61 to 23.07 wt%. The general chemical formulation of each of these systems can be written as: 9H 2 O.72SiO 2 , 18H 2 O.72SiO 2 , 36H 2 O.72SiO 2 and 72H 2 O.72SiO 2 . The supercells are built from an initial pure SiO 2 supercell, where H 2 O molecules were incrementally added, as detailed on Chapter 2. Following the same method applied for silica simulations, we varied temperature and cell volume in order to construct the respective equations of state of each of the aforementioned systems. The choice of the temperature for each composition was based on assumptions made with respect to the effects of the water concentration onto the general thermodynamic behavior. Therefore, we privileged higher temperatures for 9H 2 O.72SiO 2 and 18H 2 O.72SiO 2 , and lower temperatures for 36H 2 O.72SiO 2 and 72H 2 O.72SiO 2 . Each individual composition was simulated at three different nominal values of temperature and at ten or more density points. With respect to the density, we started with an intermediate value and then explored further lower and higher density states in order to sample liquid and gas domains.

The construction of the equation of state and location of the critical point follow the same principles applied for pure silica calculations. However, due to a limited number of simulation points on the gas phase region, the location of the gas spinodal is undefined. The low constraint of the third-order polynomial fit in such region does not allow a reliable equivalence between first-order derivative and the spinodal point, which prevents us from establishing accurate critical points for hydrated system. Nevertheless, the presence of liquid spinodals and the general assessment of fluid behavior allowed us to estimate a critical boundary based on temperature where, below such boundary, the system is in liquid-gas phase and, above this boundary, it is in supercritical state.

Structural properties were studied following the same methodology described on Chapter 2. We present results on pair distribution functions and coordination number analysis. The study on the polymerization of hydrated system is also presented, but demonstrated values are normalized due to the addition of oxygen atoms present in water molecules. Transport properties are calculated from mean square displacements and shown as diffusion coefficient values. Results of static and dynamic properties are frequently compared with pure SiO 2 data and the values of thermodynamic variables, polymerization, mean square displacement and diffusion coefficients can be found on Appendix B, D, E and F.2, respectively.

Equations of state

With the average pressure obtained from our calculations, we construct equations of state of four individual hydrated systems:

9H 2 O.72SiO 2 , 18H 2 O.72SiO 2 , 36H 2 O.72SiO 2 and 72H 2 O.72SiO 2 .
The results are presented in a pressure-density-temperature diagram where each composition comprises three different temperatures. Figure 4.1 shows the equations of state of each system, with individual panels representing results of a single composition at their varying temperatures and densities, as indicated on the top of each panel. Circles represent simulated conditions and dashed lines indicate third-order polynomial fits. As mentioned on the previous section, the positioning of gas spinodals is undefined, limiting the retrieval of information concerning the critical point location.

On the top left panel we show isotherms of 9H 2 O.72SiO 2 at 3000, 4000 and 5000 K, where densities varied from 0.34 up to 1.82 g/cm 3 and pressure varies from -0.58 to 1.74 GPa. Fitted curves at 3000 and 4000 K show a local minimum value of pressure around 1.2 to 1.3 g/cm 3 , which indicates the preferential positioning for the liquid spinodal on both isotherms. However, at 5000 K, we observe a monotonical descent of the curve, suggesting that liquid and gas phases are undistinguishable. Such behavior indicates, therefore, that the system is in supercritical state. On the top right panel we can observe a similar behavior for 18H 2 O.72SiO 2 simulations. The polynomial fits of isotherms at 3000 and 4000 K show a local minimum value around 1.6 to 1.8 g/cm 3 , which indicates the probable location of the liquid spinodal. At 5000 K, the fluid presents the same behavior as in the system previously described, where the absence of spinodals indicate the supercritical character of the system.

Bottom panels show the two compositions of hydrated silica with the largest concentrations of water. On the left side, we display results obtained from 36H 2 O.72SiO 2 calculations. Pressure-density values are indicated for temperatures at 2000, 3000 and 4000 K. Here, fitted curves at 2000 and 3000 K show a local minimum around 1.00 g/cm 3 and, therefore, liquid and gas are present at such temperatures. However, at 4000 K the continuous descent of the curve suggests that these phases cannot be separated anymore, meaning that the system is in supercritical state. The bottom right panel shows results calculated from 72H 2 O.72SiO 2 simulations.

On this system, we have the highest concentration of water of our dataset and it is where we can observe the largest variation in the thermodynamic behavior with respect to pure SiO 2 and to the previous hydrated systems. We simulated three different temperatures: 2000, 3000 and 4000 K. This composition shows separation of liquid and gas phases only at 2000 K, where it was possible to establish both spinodal points. In this scenario, the gas spinodal is placed at 0.65 g/cm 3 and the liquid spinodal is at 0.79 g/cm 3 . Given that there is no temperature simulated below 2000 K that could assist on the construction of a reliable phase diagram for this system, the critical point was not placed. However, the system shows supercritical behavior at temperatures of 3000 and 4000 K. In order to assess the sole effect of the augmentation in water concentration onto thermodynamic properties of the different systems, we compare isotherms of different compositions at equal temperatures, as shown in Figure 4.2. We present the results in four panels with fixed temperatures: 2000, 3000, 4000 and 5000 K. The upper left panel shows the equation of state of 36H 2 O.72SiO 2 (green) and 72H 2 O.72SiO 2 (pink) at 2000 K. On the upper right panel, we show isotherms of every hydrated composition at 3000 K, where blue is 9H 2 O.72SiO 2 and dark green is 18H 2 O.72SiO 2 . The lower left panel shows equations of state of anhydrous and hydrated systems at 4000 K, where dark blue represents pure SiO 2 simulation points. On the lower right panel we present isotherms at 5000 K of pure SiO 2 in green, 9H 2 O.72SiO 2 in blue and 18H 2 O.72SiO 2 in light blue.

It is possible to observe that systems with similar density at same temperature display an increase in pressure values with increasing water concentration. We highlight the diagram displaying 4000 K isotherms in the bottom left panel of Figure 4.2, where all five compositions simulated in this present study are compared. Here, it becomes evident the effect of water in the phase state of the various systems by the gradual displacement of spinodal points and increase in pressure. The change in the limit of stability of the phases suggests that systems are in different states at similar temperature-density conditions. For instance, if we consider densities around 1.5 g/cm 3 at 4000 K, pure SiO 2 is inside gas domain, whereas 18H 2 O.72SiO 2 is inside the liquid-gas coexistence dome. At this same conditions, an increase in water contents can drive the system to supercritical conditions, as demonstrate by the behavior of 36H 2 O.72SiO 2 and 72H 2 O.72SiO 2 .

0 0.5 1 1.5 -1 -0.5 0 0.5 1 Pressure (GPa) 72H 2 O.72SiO 2 36H 2 O.72SiO 2 0 0.5 1 1.5 -1 -0.5 0 0.5 1 72H 2 O.72SiO 2 36H 2 O.72SiO 2 18H 2 O.72SiO 2 9H 2 O.72SiO 2 0 0.5 1 1.5 -0.5 0 0.5 1 Density (g/cm 3 ) Pressure (GPa) 72H 2 O.72SiO 2 36H 2 O.72SiO 2 18H 2 O.72SiO 2 9H 2 O.

Influence of water on the critical point of SiO 2

On the previous section we were able to observe the variation in the physicochemical behavior of the various systems by the change in the aspect of equations of state. Although we cannot establish a reliable phase diagram for hydrated compositions, we outline a critical boundary limit based on the thermodynamic data gathered in this work in a diagram relating temperature as a function of water concentration in wt%, as shown in Figure 4.3.

We start by placing the critical point of water, as reported in the literature [START_REF] Wagner | The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use[END_REF], at 647 K and 100 wt% H 2 O. We, then, place the critical point of pure SiO 2 at 5250 K and 0 wt% H 2 O, as calculated in this work. The choice of the temperature was based on a mean value between upper (5500 K) and lower (5000 K) constraints. These two points mark the supercritical limit for pure H 2 O and pure SiO 2 systems and are represented as purple diamonds on Figure 4.3. We continue by adding the information concerning silica simulations phases at 5000, 4500 and 4500 K. At such conditions, SiO 2 is in liquid-gas regime and is, therefore, represented by black crosses.

We proceed with the analysis of the phase state of each one of the SiO 2 -H 2 O intermediary compositions. From the analysis of the respective equations of state, we estimate the phase of each system at the calculated temperatures. As previously discussed, the presence of spinodal points in a equation of state indicates that the material is in the liquid-gas regime, whereas the absence of such points leads to the assumption that the same material is now in the supercritical regime. We observe that simulations on 9H 2 O.72SiO 2 at 3000 and 4000 K present the local minima and maxima indicative of the spinodal points. The same pattern is exhibit by 18H 2 O.72SiO 2 at 3000 and 4000K, 36H 2 O.72SiO 2 at 2000 and 3000 K, and 72H 2 O.72SiO 2 at 2000 K. At these conditions, the aforementioned setups are in the liquid-gas region and are, therefore, represented as black crosses on Figure 4.3. Simulations at 5000 K on 9H 2 O.72SiO 2 and 18H 2 O.72SiO 2 , 4000 K on 36H 2 O.72SiO 2 and 3000 K on72H 2 O.72SiO 2 showed a distinct behavior, where the monotonical descent of equations of state curves indicate that these materials are in the supercritical state at such conditions. We represent these fluids in supercritical regime as green squares.

By knowing the physical states of our systems, we are now able to estimate a boundary where we are likely to find the critical limit of all silicate systems containing from 0 to 100 wt% of H 2 O. The purple dashed curve on Figure 4.3 connects the two known critical points and sets such limit of critical behavior. The gray area below this curve illustrates subcritical states, where systems are inside liquid-gas domain. Above the trendline, all fluids are in the supercritical state.

72SiO 2 9H 2 O.72SiO 2 18H 2 O.72SiO 2 36H 2 O.72SiO 2 72H 2 O.72SiO 2

Structural properties variations

From the same standpoint as applied for anhydrous simulations, we study the structural properties of the hydrated systems. We describe the behavior of each individual composition in terms of pair distribution functions, coordination numbers and polymerization analysis. Because we want to assess the influence of the presence of water into the system, we compare systems within similar quantities. For instance, in the following sections of this thesis, we are going to compare the different responses of each composition at the same simulated temperature (i.e. 4000 K) and/or at similar densities.

Pair distribution functions

The fluid structure of 9H 2 O.72SiO 2 , 18H 2 O.72SiO 2 , 36H 2 O.72SiO 2 and 72H 2 O.72SiO 2 is studied in terms of interatomic distances of oxygen-oxygen pairs at a nominal temperature of 4000 K. According to our calculations, at this temperature, systems containing 9 and 18 molecules of water are in the liquid-gas regime region, whereas systems with 36 and 72 water units are in the supercritical state. We individually represent and analyze pair distribution functions of each of the aforementioned compositions and then compare the observed behaviors in order to establish potential composition-property relational patterns. We compute PDF values at all densities and present the results as a function of radial distance in angstroms, as presented in In a global analysis of the four systems, the nearest neighbors represented in the first coordination sphere vary between 1.825 and 3.78 Å and indicate the distancing of oxygen atoms present in silica tetrahedra. Similar to the behavior observed in dry compositions, as the packing of the fluid structure becomes looser with the decrease in density, PDF peaks become wider, consistent with dissociation and vaporization of the fluid. Furthermore, distancing of atoms on tetrahedra can vary according to the density, where dense packing can force atoms to be closer and lead to a shift to the left of the peak, as clearly observed and exemplified on 72H 2 O.72SiO 2 2.77 g/cm 3 (Figure 4.4d). We also observe the placement of the first minimum above zero, evidencing an increased dynamism of the particles [START_REF] Green | Bulk properties and near-critical behaviour of SiO 2 fluid[END_REF], reaching its maximum at lowest densities.

As mentioned on Section 3.2.1, the initial portion of the pair distribution function plot reflects the repulsive forces acting between atoms. As we reach radial distances of 1 Å, we start to observe the typical peaks of diatomic oxygen, shown in detail on inserted plots of Figure 4.4. Overall, the presence of O2 is stronger at low densities, when lower pressures favor dissociation, and less present at higher densities, where it tends to stay in tetrahedra. Diatomic oxygen peak on 9H 2 O.72SiO 2 vary from 1.085 to 1.475 Å at 0.34 g/cm 3 , which represents the widest range of bond distance for the dioxygen pair in this system. Peaks present at 0.42 and 0.54 g/cm 3 are less wider and expressive. The 18H 2 O.72SiO 2 simulations follow the same trend as the aforementioned composition, where peaks can vary between 1.065 to 1.535 Å. Considering 36H 2 O.72SiO 2 and 72H 2 O.72SiO 2 , variations on peak width is less significant, where they can vary from 1.075 to 1.475 at 0.78 g/cm 3 for a system with 36 water molecules and from 1.075 to 1.575 Å at 0.43 g/cm 3 on a system with 72 water molecules. With respect to the prominence of O 2 distributions studied in this analysis, we observe an inverse relation between the presence of molecular oxygen and the concentration of water in the system, meaning that dissociation of O 2 molecules into the vapor is less frequent when quantities of water are sufficiently high. Given that we observe a steady decrease of O 2 following the trend of water concentration, we cannot associate this variation to any phase transition (i.e. liquid-vapor-supercritical), but only to the effect of H 2 O.

Coordination numbers

Here in the following section we assess the fluid organization by evaluating the structure of the atomic pair. We calculate the fraction of coordinated Si-O at each density on every system, dry and wet, at 4000 K and show the results on Figure 4 sitions hardly overlap, we chose to compare results within two distinct sets of values: a lower density that varies between 0.35 to 0.43 g/cm 3 and represented by squared symbols, and a higher density varying from 0.93 to 1.03 g/cm 3 , represented by circles. The comparison of coordination number showed little to no deviations as a direct response to variations in composition. It is important to remark that SiO 2 , 9H 2 O.72SiO 2 and 18H 2 O.72SiO 2 are in liquid-vapor domain, whereas 36H 2 O.72SiO 2 and 72H 2 O.72SiO 2 are in supercritical regime at 4000 K. The previous analysis of anhydrous silica discussed on Section 3.2.2 showed a correlation between thermodynamic properties and the structure of the fluid. However, within the conditions presented on Figure 4.5, it is not possible to observe such correlation, where the increase in dissociation is marked by an augmentation of lower coordination numbers towards vapor and supercritical states.

Nevertheless, hydrogen atoms can react with silica to form orthosilicic acid molecules, where the general formula is Si(OH) 4 . The presence of such molecule was observed in a 36H 2 O.72SiO 2 simulation at 4000 K and 1.03 g/cm 3 , as shown in Figure 4.6. Its presence in other systems as well as its quantification need to be further investigated. Considering the structure of the Si-O pair in a system rich in Si(OH) 4 , the prevailing coordination number is still 4, given that four oxygen atoms are coordinated around one silicon. However, the Si-O-H organization has an impact over the polymerization of the fluid due to a network break of silica. This fact is corroborated by the behavior observed regarding O-Si coordinations. On Figure 4.7 we show five panels with the calculated fraction of O-Si coordination number distributions at every density of anhydrous and hydrated simulations at 4000 K. Here we can observe a gradual increase in the dissociation of the fluid, marked by an increase in coordination 1 in a direct relation with the decrease in density and the increase in water concentration.

Polymerization

The presence of water can affect polymerization structure by liberation of orthosilicic acid, as demonstrated on the previous section. The mechanisms involved in this process are not going to be detailed in this work, instead our focus is in quantifying the extent of depolymerization on various hydrated fluids.

As detailed on Section 3.8, we calculate the degree of polymerization by estimation of the variable V total . It gives the relation between coordination fractions of silicon with respect to oxygen atoms (Si-O x ) over coordination fractions of oxygen with respect to silicon (O-Si y ). However, in the case of hydrated compositions, the additional oxygen atoms introduced in different proportions into the systems require that a factor of correction is added into the equation. Thus, we normalize V total for the number of oxygens contained in each individual system, defining Vn total . Here, polymerization values vary between 0.011 and 0.014, where 0.011 indicates a depolymerized system and 0.014 indicates a fully polymerized one. A table with discrete values of Si-Ox, O-Siy and Vn total is given on Appendix D.2. Similar to the patterns observed in pure silica simulations, hydrated compositions show a decrease in calculated polymerization values with density decrease. In comparison, the variations observed here are less homogeneous than those observed in anhydrous systems. Such behavior cannot be linked to any of the studied properties so far and needs further investigation.

We also observe that Vn total values are overall smaller as water concentration increases. Furthermore, high density states containing water already show a considerable degree of depolymerization. Such fact is presumably due to the formation of orthosilicic acid by hydrolysis of Si-O-Si network bonds. Considering only values for densities above 2 g/cm 3 , it is possible to observe a dramatic decrease in polymerization from SiO 2 to 72H 2 O.72SiO 2 values.

Above the supercritical line three systems are present: SiO 2 , 9H 2 O.72SiO 2 and 18H 2 O.72SiO 2 . Inside the liquid domain, these systems present an almost linear behavior. Towards the coexistence and the gas domain, such behavior becomes more erratic, where the dissociation of 9H 2 O.72SiO 2 can be found to be locally more intense than SiO 2 . Below the supercritical boundary line, we observe an important gap between 36H 2 O.72SiO 2 and 72H 2 O.72SiO 2 , meaning that the system is highly dissociated in the latter. Density (g/cm 3 )

Diffusion coefficient (m 2 /s) 72H 2 O.72SiO 2 36H 2 O.72SiO 2 18H 2 O.72SiO 2 9H 2 O.72SiO 2 72SiO 2 0 0.5 1 1.5 2 2.5 3 10 -8 10 -7
Density (g/cm 3 )

72H 2 O.72SiO 2 36H 2 O.72SiO 2 18H 2 O.72SiO 2 9H 2 O.72SiO 2 72SiO 2
Silicon Oxygen Figure 4.9: Comparison of diffusion coefficients of all simulated systems at 4000 K. Circles represent hydrated system and diamond symbols represent anhydrous silica.

Transport properties variations

The presence of water influences not only the thermodynamics, but also transport properties. Diffusion coefficients were calculated for hydrated compositions following the same principles applied for anhydrous systems. Discrete values of the bulk hydrated dataset at all temperatures are presented in Appendix F.2 and F.3. We compare self-diffusion coefficients of silicon and oxygen at 4000 K for hydrated and anhydrous compositions. The mobility of the referred atoms increase together with the concentration of H 2 O, as shown in Figure 4.9. This effect becomes more prominent at higher densities, where a vertical spread of the coefficients is observed, indicating a stronger influence of water at such conditions. At lower densities, this different becomes less significant and coefficients can overlap.

Analyzing the displacement of particles over time, we can also establish a correlation between the maximum distance travelled and the presence of water. In the presence of water, silicon and oxygen particles tend to travel longer distances with respect to dry environments. We see in Figure 4.10 the variation in the position of atoms as a function of the self-diffusion coefficient for silicon and oxygen at 4000 K of the dry silica system (blue) and 72H 2 O.72SiO 2 (pink) at different scales of time. We mark specific densities as a matter of comparison between systems in order to establish reference points which explicitly evidence the leap in travelled distances.

The lower limit of density simulated in this study for 72H 2 O.72SiO 2 is 0.43 g/cm 3 , whereas for pure SiO 2 this limit would be 0.20 g/cm 3 . Considering that low density states would implicate particles diffusing at a higher rate, we observe still a longer distance travelled by silicon in a dry system with respect to a wet system due to this density difference. Notwithstanding, the effects of water with respect to oxygen atoms at lower densities are more significant, evidenced by a larger gap between 1.00 g/cm 3 SiO 2 and 1.01 g/cm 3 72H 2 O.72SiO 2 , but also causing atoms at 0.43 g/cm 3 to travel similar distances as 0.20 g/cm 3 SiO 2 .

Up to 100 years, silicon particles in a pure SiO 2 environment would be able to travel between 4.62 and 35.8 meters, whereas in a wet environment the distances could reach between 8.74 and 25.7 meters, depending on density conditions. Considering oxygen atoms, those distances would vary from 5.36 and 39.9 meters, and from 9.96 to 37.1 meters in a dry and wet system, 2.77 g/cm 3 2.33 g/cm 3 2.77 g/cm 3

1.00 g/cm 3

1.01 g/cm 3

1.00 g/cm 3

1.01 g/cm 3 respectively. If we scale up to 100 Myr, silicon in pure SiO 2 could achieve distances between 4621.4 and 35823.7 meters, whereas oxygen could travel from 5363.2 to 39926.9 meters, also depending on the density of the system. Now regarding larger time scales, in 1 Gyr silicon particles could travel up to 30597.6 meters at 1.00 g/cm 3 in pure SiO 2 system and 53016.1 meters at 1.008 g/cm 3 in a 72H 2 O.72SiO 2 system, as for oxygen it could reach 32712.9 to 66896.1 meters, with respect to the same conditions. Discrete values of characteristic diffusion distances are presented in Table 4.1.

Ab initio study of silica and hydrated silica during giant impacts Chapter 5

Implications on lunar formation

Oxygen contents

Contents of oxygen in a particular environment are usually given in terms of its fugacity, which represents the effective pressure of the gas [START_REF] Atkins | Atkins' Physical Chemistry[END_REF]. It has been reported that variation in oxygen fugacity can affect rock melting, speciation of volatile species during outgassing, composition of planetary atmosphere, valence ratio of transition metals in minerals and more [START_REF] Lin | Oxygen controls on magmatism in rocky exoplanets[END_REF]. On terrestrial rocks, a variation between a Bowen trend and a Fenner trend can be observed as a result of oxygen abundance variation, where tholeiitic series are favored in reducing conditions and calc-alkaline series are a result of more oxidizing environments [START_REF] France | A new method to estimate the oxidation state of basaltic series from microprobe analyses[END_REF]. Experiments reproducing solar system redox conditions revealed a changing in the positioning of the liquidus curve of basalt depending on oxygen abundance [START_REF] Lin | Oxygen controls on magmatism in rocky exoplanets[END_REF], suggesting that an increase in oxygen contents can decrease the temperature of rock melting.

Measurements of paleo redox environmental conditions can be assessed by the mineralogical assembly, stablished in a series of mineral buffers described in the literature (Lindsley, 1991). The oxygen fugacity of each buffer trend is given in a gradual variation as a function of temperature, where the magnetite-heamatite (MH) assembly indicates settings with the highest oxygen partial pressures and the quartz-iron-fayalite (QIF) assembly represents the lowest value, which is associated to reducing conditions. Although calculations aiming to study the chemistry of the protolunar disk showed more oxidizing conditions during lunar formation [START_REF] Visscher | Chemistry of the Protolunar Disk, 44th Annual Lunar and Planetary Science Conference, Lunar and Planetary Science Conference[END_REF], mineral buffers and recent ab initio simulations favor a reducing environment for the Moon. Lunar mineral controls have been studied in terms of basalt and mantle compositions that bracket the redox state between IW (iron-wüstite) and IW-1 [START_REF] Stolper | Experimental petrology of eucritic meteorites[END_REF][START_REF] Shearer | Oxygen fugacity of mare basalts and the lunar mantle application of a new microscale oxybarometer based on the valence state of vanadium[END_REF][START_REF] Palme | Cosmochemical estimates of mantle composition[END_REF][START_REF] Deng | A magma ocean origin to divergent redox evolutions of rocky planetary bodies and early atmospheres[END_REF].

The concentration of molecular O 2 can be a proxy for temperature-pressure-density conditions under which the Moon condensed. In a reduced environment, the gas phase is dominated by species oxygen-poor and it is mostly formed by H 2 and CO molecules. Considering pair distribution function analysis on dry fluids presented on Section 3.2.1 of this work, such conditions could not be in place if the fluid is in supercritical state at temperatures above 5500 K, given the abundance of dioxygen in the system. The same principle can be used to assess possible density states, where the threshold is at 0.78 g/cm 3 . Therefore, the lunar condensation must have happened in a disk region where temperatures were lower than 5500 K and densities were higher than 0.78 g/cm 3 .

Assessment of the partial pressure of oxygen remain important to understand the composition and equilibration of phases in the protolunar disk, since it can alter the fractionation of species into the gas and ultimately impact mixing processes. Furthermore, assuming that concentration of molecular oxygen increases, so does the redox potential. In these conditions, O 2 will eventually be converted to H 2 O through a couple of reactions and the energy release can raise disk entropy and favor mixing.

Equilibration in the protolunar disk

Equilibration is a key aspect on lunar forming scenarios. The impact between two objects with distinct compositions calls for a postimpact process responsible for mixing components and resulting in two isotopically homogeneous bodies. Mechanisms involved in this process can dictate the degree of mixing among phases and ultimately affect possible outcomes.

The amount of mass ejected into the disk after the impact is an important factor for disk equilibration and lunar accretion. It relates to the amount of energy transferred on impact and the thermodynamic properties of materials involved. Models applied for lunar forming scenarios simulations rely on equations of state to describe the behavior of such materials under varying conditions, where M-ANEOS is the most commonly used [START_REF] Melosh | A New and Improved Equation of State for Impact Computations[END_REF]. This equation of state improved the previous monoatomic representation of the gas phase into as gas composed of diatomic molecules. We demonstrate by our structural analysis that such assumption leads to a misrepresentation of the gas, since the dissociation of the fluid often happens in clusters containing more than two atoms. Because of that, the polyatomic character of the gas requires less energy for vaporization and affects its energetic barrier of formation, which can potentially increase the amount of gas in the disk [START_REF] Canup | Dynamics of lunar formation[END_REF] and relax constraints on size of impactor and energy involved.

Furthermore, silicate materials modeled in lunar forming scenarios often are described in terms of the equation of state of forsterite, which has a critical point estimated above 6000 K [START_REF] Townsend | Liquid-vapor coexistence and critical point of Mg 2 SiO 4 from ab initio simulations[END_REF]. However, if we consider the nature of rock forming silicates with their respective critical points, as presented in Figure 3.2, together with the estimations of redox states during lunar formation and the O 2 availability at different temperatures, we limit temperatures protolunar disk temperatures in a dry scenario up to 5000 K. Considering hydrated scenarios, the lowering of critical temperatures allows systems to remain inside the supercritical regime, where mixing happens more easily, delaying the liquid-vapor fractionation and ultimately boosting chemical equilibration.

Chemical equilibration is also strongly impacted by the transfer of species among phases. Diffusion is the main mechanism involved in the transport away from the vapor-liquid boundary. Variation in the coefficients of diffusion can affect the rate on which particles move in their respective environments and, hence, favor or not equilibration. From our simulations, we were able to demonstrate how the increase in the temperature can improve diffusion on pure silicate systems. The decrease in density has a significant effect as well, where lesser dense states tend to diffuse faster. Such behavior is boosted by the presence of water. In a hydrated postimpact lunar debris disk, particles travel more rapidly with respect to anhydrous systems and, therefore, mixing and equilibration are again favored. It is worth to note that the duration of the protolunar disk can be critical for elemental distribution and equilibration, where the timescale on which diffusion happens will directly affect the atomic travelled length, both on anhydrous or hydrated phases.

Timescales of the protolunar disk have been associated with the viscosity of the fluid [START_REF] Charnoz | Evolution of the protolunar disk: Dynamics, cooling timescale and implantation of volatiles onto the earth[END_REF]. Abnormal viscous vapor can escape outwards or efficiently flow towards the surface of the Earth, leading to a significant disk mass loss. Because of that, cooling timescales at these conditions can be ten times shorter than those of non-viscous disks. Variations in the polymerization are known to affect the viscosity of fluids, where highly polymerized systems tend to present higher viscosity and will, therefore, contribute to a faster cooling. In opposition, depolymerized compositions show lower viscosity and will favor long lasting disks. The hydrolysis of silica with the formation of silicic acid break polymer networks and leads towards less viscous states.

Chapter 6

Conclusions

We simulated SiO 2 applying ab initio molecular dynamics methods over six different temperatures, from 4000 to 7000 K, in a density range from 0.2 g/cm 3 up to 2.33 g/cm 3 , corresponding to pressures up to 8.5 GPa. This region corresponds to the liquid-vapor transition and critical point region. From our simulations, we successfully built a SiO 2 phase diagram delimitating the coexistence dome within the modeled conditions. We predicted the critical point at temperatures between 5000 and 5500 K, in a density range from 0.6 to 0.85 g/cm 3 and at a pressure interval of 0.15 to 0.25 GPa. We demonstrated that above 5500 K, SiO 2 is in supercritical state where no liquid-vapor separation is possible for this material. Given the protolunar disk temperature conditions, silica condense close to the photosphere of the disk and can rapidly re-vaporize as it descends to towards the midplane, depending on disk temperature profiles.

We also demonstrated that the dissociation of SiO 2 is more complex than established previously in current equations of state [START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF]. Molecular clusters from SiO 2 dissociation form different species with varying number of atoms, where O 2 is also present. The number of atoms in the molecules present in the gas depends largely on the temperature and it is dominantly diatomic only in the supercritical regime at 7000 K. The gas phase shows clusters with coordination numbers varying from 2 to 5, depending on pressure-temperature-density conditions, which suggests that the mass of the protolunar disk could be underestimate in current lunar forming scenarios simulated that employ smoothed-particles hydrodynamics methods.

In order to characterize the polymerization of the fluid, we developed and applied a new metric that allows us to assess the degree of dissociation in terms of the relation between Si-O and O-Si coordination numbers, which indicates the potential amount of overall bonding in the fluid. This method needs to be applied carefully and is still in need of further analysis. From the analysis of the diffusion coefficients of particles in the pure SiO 2 system, we established characteristic diffusion distances for silicon and oxygen between 4000 and 7000 K for both particle types. In the estimate timeframe of complete lunar accretion (100yr), we demonstrate that silicon particles can reach a maximum distance of 0.0051 m at 4000 K and 0.0124 m at 7000 K, and oxygen particles can travel distances up to 0.0071 m at 4000K and 0.0132 m at 7000 K.

We introduced water molecules in the SiO 2 system creating four hydrated systems, in a silica-water proportion of 8:1, 4:1, 2:1 and 1:1. We simulate each system at three different temperatures, varying between 2000 and 5000 K and construct phase diagrams of each individual hydrated system. We demonstrated a variation of the distribution of phases with respect to the content of water, ultimately decreasing the critical point temperature with increasing water concentration. Due to a limited number of simulation points in the gas region, the precise location of the respective critical points was not possible of being established for these systems. However, we have estimated a limit of stability of the supercritical fluid in relation to the temperature as a function of the concentration of water. Water affects the distribution of phases and can delay the vapor-liquid separation and, consequently, the partitioning of chemical species. The presence of water also influences the characteristic diffusion distances of elements and, together with variations in phase behaviors, it can increase the mobilization of materials.

Models of lunar formation consider a giant impact as the main genetic mechanism, where two main theories prevail: the canonical model and the fast-spinning Earth model. In both cases, the elemental distribution is not completely resolved since current methods cannot precisely describe the behavior of volatiles after impact, which leads to speculation concerning their final distribution. The volatile budget is not only a matter of partitioning into liquid and vapor phases, but also of how volatiles affect disk properties as liquid fractionation occurs. Furthermore, the development of more accurate equations of state is necessary to properly describe the dissociation of molecules that can alter the total final mass. Ultimately, the understanding of the behavior of the materials in post-impact conditions play a necessary role giving clues that lead to more accurate models. Our results show that current techniques employed in the search of an accurate model for lunar formation do not satisfy the complexity of the disk dynamics and still need improvements. Nevertheless, it is clear that the complexification of models with accurate variables seem to favor the canonical impact model proposed by [START_REF] Canup | Forming a Moon with an Earth-like composition via a giant impact[END_REF].

List of Figures

A.2 Number of atoms

The convergence of the energy with respect to the number of atoms was tested with four different configurations: 48, 96, 144 and 216 atoms. The energy was normalized by the number of atoms and compared with respect with the energy value of the system with the highest number of atoms. Here, we obtained a difference of more than 0.02 eV/atom. Considering this, we chose to proceed with simulations containing 216 atoms. The simulation of larger systems is unfeasible in terms of time constraints.

Appendix B

Thermodynamic summary

Thermodynamic variables calculated by applying the methods described in Section 2. Values are given for the bulk dataset of pure SiO 2 setups (Table B.1) and hydrated systems (Table B.2 and B.3).

B.1 Pure SiO 2

1

  Earth-Moon system . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 Isotopic composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Impact Theories . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3.1.1 Canonical Impact model . . . . . . . . . . . . . . . . . . . . 8 1.3.1.2 Fast-spinning Earth impact model . . . . . . . . . . . . .

Figure 1

 1 photon-dominated complex molecules dust-gas interaction

Figure 1 . 2 :

 12 Figure 1.2: The Earth and the Moon schematic representation. Please note that the thickness of terrestrial inner layers is shown in depth with respect to the surface, whereas the lunar inner structure is indicated as radius from the center. a) Representation of the internal structure of the Earth. On the left side, values of depth are shown for the core-mantle boundary (CMB) at 2891 km, for inner core-outer core boundary (ICB) at 5150 km and for Earth's total radius (6371 km). Respective pressure, mass and volume values are also shown for each boundaries.Figure from Li and Fei (2014). b) Estimate internal structure of the Moon. The solid inner core is represented in orange with radius of 240 km, predicted fluid outer core is indicated in yellow with a radius of 330 km, partial melt enveloping the core is in dark yellow with a radius of 480 km, mantle is indicated in green and crust in grey. Figure adapted from Weber et al. (2011).

Figure 1 . 3 :

 13 Figure 1.3: Bulk Silicate Earth composition normalized to CI chondrites as a function of the half-mass condensation temperature in K. Colors represent different geochemical behaviors. Black circles are highly siderophile and chalcophile elements; open circles are halogens; red circles are moderately siderophile elements; green circles are slightly siderophile elements; yellow circles are litophile elements; blue circles are volatile litophile elements. Figure from Halliday (2014).

Figure 1 . 4 :

 14 Figure 1.4: Pressure of water atmosphere as a funcition of planet mass and initial water content produced by degassing of magma ocean. Dashed lines represent initial bulk water contents in the magma ocean. Bold lines represent the water content in the final 1 vol% liquid in magma ocean, close to full crystallization of the planet. Shaded area represents the range of water content in achondritic meteorites. Figure from Elkins-Tanton (2011).

Figure 2 . 1 :

 21 Figure 2.1: Multiscale modeling scheme on a time-length diagram. Ab initio methods are given in quantum scale, described from physical laws. At a coarser scale, molecular dynamics (MD) methods apply hard-sphere models at larger time scales. Modelisation complexity increases towards finer scales. Figure modified from Andrew J. Whittle's Research Group at MIT, Multiscale modeling of clay-water systems (accessed: June 16th 2022, https: // ajw-group.mit.edu/ multiscale-modeling-clays).

Figure 2 . 2 :

 22 Figure 2.2: Schematic representation of the DFT-MD cycle as implemented in VASP ® . Illustration from Lorenzen (2012)

Figure 2 . 3 :

 23 Figure 2.3: Initialization of the SiO 2 system structure. On the right upper side, the unit cell is represented along the a axis and on the right lower side it is represented along the c axis for a complete view of the cell. On the left side, is the final supercell with the inital atomic position adopted as starting configuration of the pure SiO 2 simulations. Red spheres represent oxygen atoms and yellow spheres represent silicon atoms.

Figure 2 . 4 :

 24 Figure 2.4: Unit cells of SiO 2 with interatomic water. On the left, the unit cell contains one molecule of water and on the right side, it contains 4 molecules of water. Hydrogen atoms are presented as white spheres.

Figure 2 . 5 :Figure 2 . 6 :

 2526 Figure 2.5: Overview of anydrous SiO 2 simulations. Each individual point represents one simulated system in its fixed density and temperature. Different colors represent different temperatures.

Figure 2

 2 Figure2.7: Phase diagram schematics in pressure-density-temperature space. Isotherms are indicated as solid lines varying with pressure (vertical axis) and density (horizontal axis). Local minima and maxima values of isotherms are spinodal points and define the limit of the gasliquid coexistence dome (dark-blue area). Dashed lines represent the Maxwell construction of equal areas where the pressure is the equilibrium pressure and sets the metastability limit of homogeneity of the gas (A) and the liquid (B). Strictly gas or liquid regions are indicated as gray areas. Critical point (CP) conditions of pressure (P c ), density (ρ c ) and temperature (T c ) are indicated in red. The supercritical fluid regime is represented as the red area and is bounded by T c .

Figure 2

 2 Figure 2.8: Representation of the pair distribution function analysis. a) Central atom in brown and neighboring atoms in green. The first coordination sphere is at radial distance r0 and the second coordination sphere is at radial distance r1. b) Schematic diagram showing the pair distribution function of an atomic pair. Dashed line shows the integral of the function and the red point represents the derivative of the curve located over the first minimum through, indicating the coordination number of the system.

Figure 3 .

 3 Figure 3.1: Equations of state of pure SiO 2 and phase diagram. The diagram shows the variation of pressure as a function of density along different isotherms. Diamond symbols represent simulation points, where the different colors represent different temperatures. Dashed colored lines along the isotherms are third-order polynomials fitted for each individual temperature. Crosses mark liquid and gas spinodal given by local minima and maxima points. The points indicated with (x) are second-order derivatives of the given isotherms and the light-red dashed line is a linear fit of such derivatives. The spinodal points along with the linear fit, bracket the critical point region, indicated as the hatched area in the diagram. Gray filled area indicates the liquid-vapor coexistence dome.

Figure 3 . 2 :

 32 Figure 3.2: Critical points of major rock-forming silicates. Diagram of temperature as a function of density with the critical point of SiO 2 as calculated in this work. Silicates described in the literature are also shown for reference: anorthite, CaAl 2 Si 2 O 8 (Kobsch, 2020); orthoclase, KAlSi 3 O 8 , and albite, NaAlSi 3 O 8[START_REF] Kobsch | The critical point and the supercritical state of alkali feldspars: Implications for the behavior of the crust during impacts[END_REF]; SiO 2[START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF][START_REF] Green | Bulk properties and near-critical behaviour of SiO 2 fluid[END_REF]; forsterite, MgSiO 3[START_REF] Xiao | Critical vaporization of MgSiO 3[END_REF]; and enstatite, Mg 2 SiO 4[START_REF] Townsend | Liquid-vapor coexistence and critical point of Mg 2 SiO 4 from ab initio simulations[END_REF].

Figure 3 . 3 :

 33 Figure 3.3: Pair distribution function analysis for the full interval of simulated densities at selected temperatures. Left side panels show the distribution at 4000 K and panels on the right side shows distributions at 7000 K. Upper panels show Si-O pairing and lower panels shows O-O pairing. The span of densities at 4000 K goes from 2.33 to 0.20 g/cm 3 . The span of densities at 7000 K goes from 2.33 to 0.27 g/cm 3 .

Figure 3 . 4 :

 34 Figure 3.4: Pair distribution function of the Si-O and O-O pair at different temperatures and fixed densities. Upper panels show Si-O pairing and lower panels show O-O pairing for all temperatures in both cases.Left panels show the variation of the distribution at 2.33 g/cm 3 , whereas right panels show the variation of the distribution at 0.27 g/cm 3 . The total radial distance is giving in half of the size of the simulation box. Therefore, the horizontal axis for the 0.27 g/cm 3 diagrams go up to 15 Å given that the lattice for SiO 2 at such conditions is 30 Å, whereas for the diagram with density of 2.33 g/cm 3 the radial distance goes up to 7 since the lattice for those simulations were of 14.554 Å.

Figure 3 . 5 :

 35 Figure 3.5: Position of the first peak of the Si-O pair distribution function of the full extent of simulated temperatures as a function of energy. Dashed lines represent isochoric variation of the position-energy relation.

Figure 3 . 7 :

 37 Figure 3.7: Proportion of coordination numbers per temperature at fixed density. Left panel shows the variation in proportion of coordination numbers at a density of 2.33 g/cm 3 from 4000 to 7000 K, whereas the right panel show the variation at a density of 0.27 g/cm 3 at the equal temperatures.

Figure 3 . 8 :

 38 Figure 3.8: Polymerization of the silicate fluid in different domains. The total ratio of Si-O x over O-Si y is represented as a function of density.Gray areas show the gas domain region at low densities and the liquid domain region at higher densities. In the middle of these two regions is the coexistence area. The dashed line marks the limit of the supercritical regime, where all points below this line are above critical conditions.

Figure 3 Figure 3 .

 33 Figure 3.9: Visualization of different degrees of polymerization. (a) Pure SiO 2 at 4000 K and 2.33 g/cm 3 , where V total = 2. (b) Pure SiO 2 at 7000 K and 0.27 g/cm 3 , where V total = 1.498.

Figure 3 .

 3 Figure 3.11: Diffusion coefficients of Si and O. Diagrams are presented as a function of density for all temperatures. Left side panel: silicon; right side panel: oxygen.

Figure 3 .

 3 Figure 3.12: Diffusion coefficients of Si and O as a function of temperature at fixed densities. Left side panel shows coefficients for silicon at 2.33 g/cm 3 (brown), 1.00 g/cm 3 (pink) and 0.27 g/cm 3 (green). Right side panel shows coefficients for oxygens at same densities.

Figure 3 . 13 : 2

 3132 Figure 3.13: Displacement of atoms over time. Characteristic diffusion distance as a function of self-diffusion coefficients. Silicon distances are shown on the left-side and oxygen distances are shown on the right-side. Comparison between displacements at 4000 (blue) and 7000 K (red) in different time intervals.

Figure 4 . 2 :

 42 Figure 4.2: Comparison between pressure-density trends at same nominal temperatures with respect to the concentration of H 2 O in the system. Circles represent simulation points and dashed-lines represent third-order polynomial fits.

Figure 4 . 3 :

 43 Figure 4.3: The critical point of SiO 2 and the critical point of H 2 O (Wagner and Pruß, 2002): trend of supercritical state between pure SiO 2 and pure H 2 O with respect to the concentration of H 2 O in the silicatic system. The critical points are represented with purple diamond symbols; green squares represent simulations manifesting supercritical behavior; black crosses represent simulations where spinodals are present. The green filled squares indicates compositions where the fluid would be found below critical conditions. Colored arrows indicate points related to each individual composition.

Figure 4 . 4 .

 44 Each panel illustrate a single composition, where upper left panel (a) shows 9H 2 O.72SiO 2 , upper right panel (b) shows 18H 2 O.72SiO 2 , lower left panel (c) shows 36H 2 O.72SiO 2 and lower right panel (d) shows 72H 2 O.72SiO 2 .

Figure 4 . 4 :

 44 Figure 4.4: Pair distribution functions of oxygen-oxygen pairs at 4000 K for all hydrated systems. a) 9H 2 O.72SiO 2 ; b) 18H 2 O.72SiO 2 ; c) 36H 2 O.72SiO 2 ; d) 72H 2 O.72SiO 2 . Inserted panels show in detail the distribution between 0.8 and 2 Å indicating the presence of O 2 molecules.

Figure 4 . 5 :

 45 Figure 4.5: Comparison of coordination number proportions of Si-O in anhydrous and hydrated systems at 4000 K at similar densities. Circles represent values for densities between 0.93 and 1.03 g/cm 3 and squares represent values for densities between 0.35 and 0.43 g/cm3 . Each color indicates one individual system, where dark blue is pure SiO 2 , light blue is 9H 2 O.72SiO 2 , dark green is 18H 2 O.72SiO 2 , light green is 36H 2 O.72SiO 2 and pink is 72H 2 O.72SiO 2 .

Figure 4 .

 4 Figure 4.6: Orthosilicic acid molecule spotted in a simulation of 36H 2 O.72SiO 2 at 1.03 g/cm 3 and 4000 K. Spheres represent individual atoms, where yellow is silicon, red is oxygen and white is hydrogen.

Figure 4 . 7 :

 47 Figure 4.7: Distribution of O-Si coordinations in anhydrous and hydrated systems at 4000 K at all simulated densities. Each panel shows an individual composition: SiO 2 on the top left, 9H 2 O.72SiO 2 on the top middle, 18H 2 O.72SiO 2 on the top right, 36H 2 O.72SiO 2 on the bottom middle and 72H 2 O.72SiO 2 on the bottom right of the figure.

Figure 4 . 8 :

 48 Figure 4.8: Polymerization analysis of pure SiO 2 and hydrated systems at 4000 K. The Vn total value is the V total value normalized by the number of oxygens in each system.

Figure 4 .

 4 Figure 4.10: Displacement of silicon and oxygen in dry and wet systems atoms in four different intervals of time. The distance is presented as a function of the self-diffusion coefficients of each atom at given pressure-density conditions. Pure SiO 2 results are presented in blue and 72H 2 O.72SiO 2 results are presented in pink. Dashed lines mark the characteristic diffusion distance of specified densities within the same color code. Considered time scales are 100 yr, 100 kyr, 100 Myr and 1 Gyr.

Figure A. 2 :

 2 Figure A.2: Convergence test of energy with respect to the number of atoms. The energy is given per atom for a given simulated system. We test systems containing 48, 96, 144 and 216 atoms.

Figure E. 2 :

 2 Figure E.2: Mean square displacements of hydrous fluids. Curves represent the displacement of individual densities at given temperatures, for silicon (left panels), oxygen (middle panels) and hydrogen (right panels). Colors indicate temperature, where red is 5000 K, orange is 4000 K, green is 3000 K and blue is 2000 K.

  Figure F.1: Complete dataset of self-diffusion coefficients of Si, O and H for every hydrated silicate system at all temperatures.

2.1 Multiscale modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Ab initio molecular dynamics simulations . . . . . . . . . . . . . . . 2

  .2.1 Schrödinger equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Born-Oppenheimer approximation . . . . . . . . . . . . . . . . . . . . 2.2.3 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3.1 Hohenberg and Kohn theorem . . . . . . . . . . . . . . . . .

	2.2.3.2 Kohn-Sham equations and exchange-correlation functional .
	2.2.4 Molecular dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3 Initial configurations . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Computational details . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Equations of state and phase diagrams . . . . . . . . . . . . . . . . 2.6 Structural properties . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Polymerization analysis . . . . . . . . . . . . . . . . . . . . . . . . . .

	2.6.1 Atomic pair distribution function and coordination numbers . . . . . .
	2.6.2

2.7 Transport properties . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 2 .

 2 1: System size scaling of various methods and their maximum number of atoms. Table from[START_REF] Lee | Computational materials science: an introduction[END_REF].

	Methods	Scaling	N max (atoms)
	Schrödinger	O(e N )	1
	Hartree-Fock	O(N 4 )	∼50
	DFT	O(N)-O(N 3 )	∼200
	MD	O(N)	10 (7-11)

Table 2 . 2 :

 22 Summary of calculated systems at their fixed conditions.

	SiO 2	4000 -7000	2.33 -0.20	216	1152
	9H 2 O.72SiO 2	3000 -5000	2.72 -0.34	243	1224
	18H 2 O.72SiO 2	3000 -5000	2.29 -0.35	270	1296
	36H 2 O.72SiO 2	2000 -4000	2.45 -0.38	324	1440
	72H 2 O.72SiO 2	2000 -4000	2.77 -0.43	432	1728

System

Temperatures (K) Densities (g/cm 3 ) N • atoms N • electrons

Chapter 3 Pure SiO 2 system Contents 3.1 Phase diagram and the critical point

. . . . . . . . . . . . . . . . . . 26 3.2 Structure of the fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2.1 Pair distribution functions . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2.2 Coordination numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2.3 Polymerization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Dynamic behavior of atoms . . . . . . . . . . . . . . . . . . . . . . . 36

Table 3 .

 3 1: Characteristic diffusion distances travelled of silicon and oxygen atoms in different time intervals at 4000 and 7000 K. Values referring to 2.33 g/cm 3 simulations are indicated by x min and values referring to 0.27 g/cm 3 simulations are given as x max . Distances are given in meters.

			4000 K				7000 K
	Time	Silicon		Oxygen		Silicon		Oxygen
	(years)	x min	x max	x min	x max	x min	x max	x min	x max
	100	4.62E+00 3.58E+01 5.36E+00 3.99E+01 2.06E+01 6.95E+01 2.26E+01 7.43E+01
	100k	1.46E+02 1.13E+03 1.70E+02 1.26E+03 6.51E+02 2.20E+03 7.15E+02 2.35E+03
	100M	4.62E+03 3.58E+04 5.36E+03 3.99E+04 2.06E+04 6.95E+04 2.26E+04 7.43E+04
	1G	1.46E+04 1.13E+05 1.70E+04 1.26E+05 6.51E+04 2.20E+05 7.15E+04 2.35E+05

  .5. Different colors represent different compositions, where dark blue is SiO 2 , light blue is 9H 2 O.72SiO 2 , dark green is 18H 2 O.72SiO 2 , light green is 36H 2 O.72SiO 2 and pink is 72H 2 O.72SiO 2 . Given that densities among compo-

	g(r) O-O	5 10 15	a	0.8 0 0.3 0.6 0.9 9H 2 O.72SiO 2 1.2	1.6	2	2.72 2.21 1.82 1.52 1.28 1.09 0.93 0.80 g/cm 3	4 8 12	b	0.8 0 0.3 0.6 0.9 18H 2 O.72SiO 2 1.2	1.6	2	0.83 2.81 1.32 1.13 0.97 1.57 1.89 2.29 g/cm 3
										0.70							0.73
										0.54							0.56
		0								0.42 0.34	0						0.44 0.35
	g(r) O-O	10 15	c	36H 2 O.72SiO 2 0.9 0.6 0.8 1.2 0 0.3	1.6	2	2.45 2.02 1.68 1.42 1.20 1.03 0.89	4 6 8 10	d	72H 2 O.72SiO 2 0.9 0.6 0.8 1.2 0 0.3	1.6	2	2.77 2.28 1.90 1.60 1.36 1.17 1.01
		5								0.78							0.88
										0.60 0.47	2						0.68 0.53
		0								0.38	0						0.43
		0	1	2	3	4	5	6	7		0	1	2	3	4	5	6	7
				Radial distance (Å)						Radial distance (Å)

  3 . Each color indicates one individual system, where dark blue is pure SiO 2 , light blue is 9H 2 O.72SiO 2 , dark green is 18H 2 O.72SiO 2 , light green is 36H 2 O.72SiO 2 and pink is 72H 2 O.72SiO 2 .

Table 4 .

 4 1: Characteristic diffusion distances of silicon and oxygen of 72H 2 O.72SiO 2 at 4000 K.Values from simulations at 0.43 g/cm 3 are represented as x min , whereas values from simulations at 2.77 g/cm 3 are indicated as x max . Values are given in meters.

	Time	Silicon		Oxygen
	(years)	x min	x max	x min	x max
	100	8.74E+00 2.57E+01 9.96E+00 3.71E+01
	100k 2.76E+02 8.12E+02 3.15E+02 1.17E+03
	100M 8.74E+03 2.57E+04 9.96E+03 3.71E+04
	1G	2.76E+04 8.12E+04 3.15E+04 1.17E+05

  1.1 Protoplanetary disk schematic structure . . . . . . . . . . . . . . . . . . . . . . 1.2 The Earth and the Moon internal structure . . . . . . . . . . . . . . . . . . . . . 1.3 Bulk Silicate Earth composition . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Residual water content in magma oceans. . . . . . . . . . . . . . . . . . . . . . . 2.1 Multiscale modeling scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Schematic representation of the DFT-MD cycle as implemented in the VASP ® code. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Initial SiO 2 supercell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Example of initial unit cells including water . . . . . . . . . . . . . . . . . . . . 2.5 Overview of anydrous SiO 2 simulations. . . . . . . . . . . . . . . . . . . . . . . . 2.6 Overview of hydrated system simulations. . . . . . . . . . . . . . . . . . . . . . 2.7 Phase diagram schematics in pressure-density-temperature space. . . . . . . . . 2.8 Representation of the pair distribution function analysis. . . . . . . . . . . . . . 3.1 Equations of state of anhydrous SiO 2 and phase diagram. . . . . . . . . . . . . . 3.2 Critical points of major rock-forming silicates. . . . . . . . . . . . . . . . . . . . 3.3 Pair distribution function analysis at fixed temperatures. . . . . . . . . . . . . . 3.4 Pair distribution function analysis at fixed densities. . . . . . . . . . . . . . . . . 3.5 Position of the first peak of the Si-O PDF as a function of energy. . . . . . . . . 3.6 Fraction of coordination numbers per density at fixed temperature. . . . . . . . 3.7 Fraction of coordination numbers per temperature at fixed densities. . . . . . . . 3.8 Polymerization of the silicate fluid in different domains. . . . . . . . . . . . . . . 3.9 Visualization of different degrees of polymerization. . . . . . . . . . . . . . . . . 3.10 Polymerization as a function of pressure at different temperatures. . . . . . . . . 3.11 Diffusion coefficients of Si and O at various temperatures. . . . . . . . . . . . . . 3.12 Diffusion coefficients of Si and O at fixed densities. . . . . . . . . . . . . . . . . 3.13 Characteristic diffusion distance of Si and O as a function of self-diffusion coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Equations of state of four hydrated silicate systems. . . . . . . . . . . . . . . . .

4.2 Comparison of equations of state at fixed temperatures. . . . . . . . . . . . . . . 4.3 The trend of the supercritical state between SiO 2 and H 2 O critical points. . . . . 4.4 Pair distribution functions of oxygen-oxygen pairs at 4000 K for all hydrated systems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table B .

 B 1: Summary of thermodynamic variables at all densities of pure SiO 2 simulations.

	Temperature (K)	Density (g/cm 3 )	Pressure (GPa)	Internal Energy (eV/atom)	Temperature (K)	Density (g/cm 3 )	Pressure (GPa)	Internal Energy (eV/atom)
	4000	2.33	2.7077	-7.0110	5000	0.20	0.0866	-5.5974
	4000	2.00	0.3693	-6.9642	5500	2.33	5.0527	-6.3046
	4000	1.75	-0.3578	-6.9067	5500	2.00	2.4481	-6.2090
	4000	1.46	-0.3980	-6.8459	5500	1.75	1.4891	-6.1054
	4000	1.23	-0.3462	-6.7659	5500	1.46	0.8267	-6.0008
	4000	1.00	-0.2608	-6.6950	5500	1.23	0.5727	-5.9049
	4000	0.90	-0.1846	-6.6506	5500	1.00	0.3688	-5.8022
	4000	0.68	-0.0950	-6.6027	5500	0.90	0.3453	-5.7727
	4000	0.52	-0.0640	-6.5642	5500	0.78	0.2877	-5.7213
	4000	0.46	-0.0495	-6.5373	5500	0.68	0.2708	-5.6503
	4000	0.41	-0.0288	-6.4677	5500	0.59	0.2594	-5.6031
	4000	0.37	-0.0255	-6.4570	5500	0.52	0.2301	-5.5733
	4000	0.33	-0.0064	-6.3724	5500	0.46	0.2131	-5.5119
	4000	0.27	-0.0027	-6.4273	5500	0.41	0.1851	-5.4856
	4000	0.20	0.0016	-6.3376	5500	0.37	0.1797	-5.4898
	4500	2.33	3.2169	-6.7764	5500	0.33	0.1650	-5.4278
	4500	2.00	0.9215	-6.6943	5500	0.27	0.1579	-5.3912
	4500	1.75	0.2119	-6.6326	5500	0.20	0.1272	-5.2682
	4500	1.46	-0.0605	-6.5122	6000	2.33	6.1940	-6.0965
	4500	1.23	-0.0611	-6.4224	6000	2.00	3.3847	-5.9952

Table B .1 continued from previous page

 B

	Temperature (K)	Density (g/cm 3 )	Pressure (GPa)	Internal Energy (eV/atom)	Temperature (K)	Density (g/cm 3 )	Pressure (GPa)	Internal Energy (eV/atom)
	4500	1.00	-0.0605	-6.3715	6000	1.46	1.2663	-5.7939
	4500	0.90	-0.0351	-6.3110	6000	1.23	0.8416	-5.7219
	4500	0.68	0.0073	-6.2273	6000	1.00	0.6035	-5.6162
	4500	0.52	0.0277	-6.1482	6000	0.78	0.4744	-5.4752
	4500	0.46	0.0572	-6.0556	6000	0.68	0.4089	-5.4208
	4500	0.41	0.0276	-6.1224	6000	0.59	0.3818	-5.3719
	4500	0.33	0.0404	-6.0298	6000	0.52	0.3082	-5.4002
	4500	0.27	0.0314	-6.0412	6000	0.41	0.2704	-5.3064
	4500	0.20	0.0357	-5.9858	6000	0.33	0.2320	-5.2599
	5000	2.33	4.0394	-6.5386	6000	0.27	0.2036	-5.1920
	5000	2.00	1.5888	-6.4377	7000	2.33	8.4697	-5.7463
	5000	1.75	0.7823	-6.3510	7000	2.00	5.1673	-5.6333
	5000	1.46	0.3857	-6.2393	7000	1.46	2.2210	-5.4431
	5000	1.23	0.2134	-6.1467	7000	1.23	1.5643	-5.3401
	5000	1.00	0.1772	-6.0483	7000	1.00	1.0994	-5.2424
	5000	0.90	0.1489	-6.0195	7000	0.78	0.7933	-5.1398
	5000	0.68	0.1296	-5.9359	7000	0.68	0.6554	-5.1175
	5000	0.52	0.1295	-5.7999	7000	0.59	0.5969	-5.0580
	5000	0.46	0.1370	-5.7715	7000	0.52	0.5046	-5.0265
	5000	0.41	0.1243	-5.7081	7000	0.41	0.4140	-4.9467
	5000	0.33	0.1163	-5.6388	7000	0.33	0.3371	-4.8846
	5000	0.27	0.0909	-5.7488	7000	0.27	0.2939	-4.8171
	B.2 Hydrated systems

Table B .

 B 2: Summary of thermodynamic variables at all densities of 9H 2 O.72SiO 2 and 18H 2 O.72SiO 2 simulations. SiO x ) V (OSi y ) V total T (K) Density (g/cm 3 ) V (SiO x ) V (OSi y ) V total used in the calculation of the Vn total variable and estimation of ation of hydrated fluids.

	9H 2 O.72SiO 2 Density (g/cm 3 ) Pressure (GPa) 2.72 9.2315 2.21 2.1894 1.82 0.4683 1.52 -0.2960 1.28 -0.5316 1.09 -0.5828 0.93 -0.3806 0.81 -0.3005 0.70 -0.2401 0.54 -0.1574 0.42 -0.1165 0.34 -0.0426 2.72 9.5219 2.21 2.7279 1.82 0.2806 1.52 -0.1073 1.28 -0.1487 1.09 -0.2214 0.93 -0.1470 0.81 -0.1074 0.70 -0.0395 0.54 -0.0367 0.42 -0.0204 0.34 0.0022 2.72 11.2396 2.21 4.1241 1.82 1.7441 1.52 0.7954 1.28 0.4046 1.09 0.3115 0.93 0.2741 0.81 0.2045 0.70 0.1793 0.54 0.1542 0.42 0.1366 0.34 0.1252 Table B.3: Summary of thermodynamic variables at all densities of 36H 2 O.72SiO 2 and 18H 2 O.72SiO 2 Temperature (K) Internal Energy (eV/atom) Temperature (K) Density (g/cm 3 ) Pressure (GPa) Internal Energy (eV/atom) 3000 -6.9996 3000 2.81 12.6817 -6.6934 3000 -7.0510 3000 2.29 4.1784 -6.7701 3000 -7.0617 3000 1.89 1.0657 -6.8031 3000 -7.0167 3000 1.57 -0.0608 -6.8113 3000 -7.0109 3000 1.32 -0.3871 -6.7583 3000 -6.9721 3000 1.13 -0.3488 -6.7505 3000 -6.9975 3000 0.97 -0.4503 -6.7447 3000 -6.9210 3000 0.83 -0.3336 -6.7096 3000 -6.9197 3000 0.73 -0.1655 -6.7202 3000 -6.9209 3000 0.56 -0.1251 -6.7130 3000 -6.8771 3000 0.44 -0.1322 -6.6926 3000 -6.9219 3000 0.35 -0.0697 -6.6084 4000 -6.6679 4000 2.81 13.2465 -6.4068 4000 -6.7208 4000 2.29 4.6785 -6.4576 4000 -6.6789 4000 1.89 1.3782 -6.4479 4000 -6.5984 4000 1.57 0.2894 -6.4083 4000 -6.5377 4000 1.32 -0.0404 -6.3550 4000 -6.5138 4000 1.13 -0.0322 -6.2966 4000 -6.4722 4000 0.97 -0.0363 -6.2807 4000 -6.4121 4000 0.83 -0.0179 -6.2414 4000 -6.4370 4000 0.73 -0.0398 -6.2396 4000 -6.3524 4000 0.56 -0.0070 -6.1778 4000 -6.3196 4000 0.44 -0.0100 -6.1470 4000 -6.2126 4000 0.35 0.0138 -6.0733 5000 -6.3032 5000 2.81 15.4476 -6.0795 5000 -6.2694 5000 2.29 6.1715 -6.0711 5000 -6.1586 5000 1.89 2.6336 -6.0077 5000 -6.0582 5000 1.57 1.3080 -5.9222 5000 -5.9951 5000 1.32 0.7585 -5.8535 5000 -5.9350 5000 1.13 0.4863 -5.8360 5000 -5.8408 5000 0.97 0.3528 -5.7223 5000 -5.8120 5000 0.83 0.2883 -5.6816 5000 -5.7372 5000 0.73 0.2432 -5.6332 5000 -5.7172 5000 0.56 0.1981 -5.5746 5000 -5.5804 5000 0.44 0.1791 -5.4868 5000 -5.5078 5000 0.35 0.1365 -5.4735 72H 2 O.72SiO 2 simulations. 36H 2 O.72SiO 2 72H 2 O.72SiO 2 Temperature (K) Density (g/cm 3 ) Pressure (GPa) Internal Energy (eV/atom) Temperature (K) Density (g/cm 3 ) Pressure (GPa) Internal Energy (eV/atom) 2000 2.45 8.1903 -6.5350 2000 2.28 8.4823 -6.0375 2000 2.02 2.8131 -6.5598 2000 1.90 3.5643 -6.0820 2000 1.68 0.9752 -6.5600 2000 1.60 1.3130 -6.0868 1.42 0.1136 -6.5815 2000 1.36 0.5190 0.41 1.615 2.576 2.765 3.795 0.41 1.588 2.730 2.810 4.108 Table D.1 continued from previous page 1.20 3.785 1.580 2.396 1.36 3.860 1.450 2.662 -6.0903 1.20 -0.3530 -6.5687 2000 1.17 0.1334 -6.0749 1.03 -0.7476 -6.5427 2000 1.01 0.0235 -6.0758 0.89 -0.5618 -6.5391 2000 0.88 -0.0730 -6.0694 0.78 -0.4151 -6.5335 2000 0.68 -0.1851 -6.0688 0.60 -0.2942 -6.5286 2000 0.53 -0.0448 -6.0741 0.47 -0.1380 -6.5231 2000 0.43 -0.0838 -6.0627 0.38 -0.0884 -6.5186 3000 2.77 20.2077 -5.7673 2.45 8.7030 -6.3476 3000 2.28 9.4887 -5.8795 2.02 3.1716 -6.4179 3000 1.90 4.2268 -5.9192 1.68 1.0607 -6.4399 3000 1.60 1.9811 -5.9495 1.42 0.3925 -6.4180 3000 1.36 0.8725 -5.9541 1.20 -0.2613 -6.4074 3000 1.17 0.5333 -5.9391 1.03 -0.3431 -6.4063 3000 1.01 0.2623 -5.9338 0.89 -0.1956 -6.3966 3000 0.88 0.1530 -5.9327 0.78 -0.1355 -6.3846 3000 0.68 0.0936 -5.9422 0.60 -0.0502 -6.3888 3000 0.53 0.0358 -5.9403 0.47 -0.0528 -6.3289 3000 0.43 0.0336 -5.9015 0.38 -0.0375 -6.3017 4000 2.77 22.7022 -5.5590 2.45 9.3557 -6.0777 4000 2.28 10.7820 -5.6331 2.02 3.7986 -6.1069 4000 1.90 5.1142 -5.6682 1.68 1.4695 -6.0876 4000 1.60 2.6280 -5.6740 1.42 0.6898 -6.0463 4000 1.36 1.4392 -5.6623 1.20 0.3585 -6.0243 4000 1.17 0.9040 -5.6284 1.03 0.1467 -5.9965 4000 1.01 0.5772 -5.6174 0.89 0.1075 -5.9762 4000 0.88 0.4216 -5.6123 0.78 0.0910 -5.9346 4000 0.68 0.2505 -5.5702 0.60 0.0534 -5.9390 4000 0.53 0.1589 -5.5550 0.47 0.0509 -5.8216 4000 0.43 0.1278 -5.5287 0.38 0.0387 -5.8487 ----0.33 1.612 2.562 2.743 3.884 0.27 1.611 2.561 2.749 3.895 0.20 1.604 2.573 2.753 3.931 2.33 1.638 2.475 2.672 3.755 2.00 1.637 2.515 2.695 3.769 1.75 1.635 2.540 2.712 3.786 1.46 1.631 2.542 2.722 3.826 1.23 1.629 2.587 2.730 3.843 1.00 1.621 2.583 2.738 3.877 0.90 1.622 2.574 2.742 3.867 0.68 1.614 2.618 2.775 3.985 0.52 1.607 2.627 2.760 3.965 0.46 1.604 2.636 2.758 3.950 0.41 1.602 2.617 2.776 3.987 0.33 1.594 2.640 2.783 3.998 0.27 1.600 2.612 2.785 0.20 1.587 2.676 2.782 2.33 1.639 2.497 2.673 2.00 1.636 2.535 2.711 1.75 1.633 2.556 2.726 3.840 1.46 1.629 2.580 2.740 3.871 1.23 1.624 2.605 2.746 3.897 1.00 1.618 2.621 2.766 3.924 0.90 1.614 2.612 2.766 3.933 0.78 1.613 2.642 2.770 3.962 0.68 1.606 2.665 2.775 3.978 0.59 1.604 2.658 2.775 3.895 0.52 1.602 2.658 2.779 3.993 0.46 1.597 2.685 2.791 4.029 0.41 1.593 2.712 2.798 4.032 0.37 1.589 2.702 2.792 4.022 0.33 1.589 2.741 2.799 4.083 0.27 1.582 2.705 2.802 4.008 0.20 1.573 2.746 2.824 4.064 2.33 1.640 2.511 2.675 3.813 2.00 1.638 2.558 2.707 3.861 1.46 1.626 2.600 2.757 3.902 1.23 1.625 2.581 2.775 1.00 1.614 2.637 2.774 3.971 0.78 1.605 2.661 2.835 3.865 0.68 1.603 2.681 2.825 3.885 0.59 1.596 2.693 2.806 4.010 0.52 1.597 2.672 2.793 4.073 4500 0.68 3.343 1.733 1.929 6000 1.00 3.024 1.620 4500 0.52 3.240 1.693 1.914 6000 0.78 2.764 1.551 4500 0.46 3.074 1.641 1.873 6000 0.68 2.690 1.512 4500 0.41 3.194 1.670 1.913 6000 0.59 2.561 1.489 4500 0.33 2.997 1.629 1.840 6000 0.52 2.567 1.476 1.42 3.835 1.603 2.393 1.60 3.911 1.451 2.696 1.739 1.68 3.901 1.624 2.402 1.90 4.044 1.496 2.703 1.719 2.02 4.011 1.659 2.417 2.28 4.213 1.545 2.726 1.779 2.45 4.182 1.720 2.432 2.77 4.607 1.642 2.806 1.782 36H 2 O.72SiO 2 72H 2 O.72SiO 2 1.867 3.895 4000 1.46 3.826 1.917 1.996 5500 1.75 3.619 1.845 4000 1.23 3.774 1.893 1.994 5500 1.46 3.448 1.779 1.938 4000 1.00 3.704 1.869 1.982 5500 1.23 3.316 1.732 1.915 4000 0.90 3.664 1.857 1.973 5500 1.00 3.129 1.667 1.877 4000 0.68 3.600 1.823 1.974 5500 0.90 3.052 1.633 1.869 4000 0.52 3.595 1.806 1.991 5500 0.78 2.981 1.616 4000 0.46 3.517 1.813 1.939 5500 0.68 2.870 1.579 4000 0.41 3.449 1.770 1.949 5500 0.59 2.754 1.545 4000 0.37 3.466 1.766 1.963 5500 0.52 2.684 1.529 4000 0.33 3.342 1.746 1.914 5500 0.46 2.609 1.483 4000 0.27 3.379 1.781 1.897 5500 0.41 2.546 1.469 4000 0.20 3.259 1.709 1.907 5500 0.37 2.519 1.458 4500 2.33 4.029 2.015 1.999 5500 0.33 2.449 1.447 4500 2.00 3.885 1.950 1.993 5500 0.27 2.281 1.380 4500 1.75 3.809 1.913 1.992 5500 0.20 2.028 1.295 4500 1.46 3.687 1.855 1.988 6000 2.33 3.927 1.985 4500 1.23 3.586 1.826 1.964 6000 2.00 3.737 1.904 4500 1.00 3.519 1.797 1.958 6000 1.46 3.346 1.745 4500 0.90 3.455 1.765 1.957 6000 1.23 3.166 1.669 1.897 0.34 3.382 1.639 2.063 0.35 3.416 1.584 2.157 1.918 0.42 3.568 1.719 2.076 0.44 3.553 1.633 2.175 1.963 0.54 3.598 1.724 2.087 0.56 3.575 1.634 2.188 1.978 0.70 3.656 1.737 2.105 0.73 3.669 1.651 2.222 1.566 0.81 3.662 1.728 2.119 0.83 3.666 1.665 2.202 1.654 0.93 3.707 1.757 2.110 0.97 3.723 1.677 2.221 1.693 1.09 3.764 1.780 2.114 1.13 3.742 1.686 2.219 1.728 1.28 3.788 1.798 2.106 1.32 3.811 1.719 2.217 1.733 1.52 3.849 1.815 2.120 1.57 3.872 1.741 2.223 1.759 1.82 3.923 1.854 2.116 1.89 3.945 1.778 2.218 1.756 2.21 4.027 1.899 2.120 2.29 4.049 1.817 2.229 1.782 2.72 4.204 1.982 2.121 2.81 4.324 1.937 2.233 1.818 9H 2 O.72SiO 2 18H 2 O.72SiO 2 1.845 (g/cm 3 ) V (SiO x ) V (OSi y ) V total Density (g/cm 3 ) V (SiO x ) V (OSi y ) V total mmary of values Density 1.962 3.817 4000 1.75 3.896 1.949 1.999 5500 2.00 3.768 1.906 1.977 3.786 4000 2.00 3.945 1.973 1.999 5500 2.33 3.969 1.995 1.989 5000 0.27 2.693 1.510 1.783 7000 0.27 1.890 1.262 1.498 4.024 4000 2.33 4.023 2.011 2.000 5000 0.20 2.448 1.463 1.674 5000 0.33 2.574 1.470 1.751 7000 0.33 2.033 1.309 1.553 3.925 0.33 1.586 2.686 2.810 4.075 0.27 1.574 2.736 2.829 3.995 2.33 1.642 2.556 2.674 3.851 2.00 1.634 2.562 2.730 3.892 1.46 1.628 2.609 2.778 3.955 1.23 1.621 2.622 2.789 4.006 1.00 1.613 2.655 2.805 0.78 1.603 2.687 2.818 4.083 0.68 1.605 2.681 2.818 0.59 1.597 2.735 2.822 0.52 1.594 2.740 2.832 0.41 1.590 2.777 2.829 0.33 1.577 2.828 2.852 0.27 1.570 2.838 2.856 4.655 T (K) Density (g/cm 3 ) V (SiO x ) V (OSi y ) V total T (K) 5000 0.41 2.724 1.533 1.778 7000 0.41 2.175 1.354 1.606 (g/cm 3 ) V (SiO x ) V (OSi y ) V total 5000 0.46 2.818 1.572 1.793 7000 0.52 2.333 1.399 1.668 Density 5000 0.52 2.902 1.590 1.825 7000 0.59 2.403 1.430 1.680 4.108 Table D.1: Summary of the polymerization data for every SiO 2 system simulated. 5000 0.68 3.120 1.647 1.894 7000 0.68 2.467 1.450 1.702 4.204 5000 0.90 3.216 1.699 1.894 7000 0.78 2.559 1.484 1.725 4.218 5000 1.00 3.294 1.719 1.917 7000 1.00 2.775 1.557 1.782 4.054 5000 1.23 3.452 1.781 1.938 7000 1.23 2.949 1.619 1.822 4.112 5000 1.46 3.548 1.809 1.961 7000 1.46 3.167 1.692 1.872 Polymerization data 5000 1.75 3.714 1.878 1.977 7000 2.00 3.578 1.852 1.932 4.035 Appendix D T (K) Density (g/cm 3 ) 1.03 3.757 1.566 2.398 1.17 3.774 1.418 2.662 0.89 3.725 1.568 2.375 1.01 3.732 1.401 2.663 V (4500 0.27 3.011 1.602 1.880 6000 0.41 2.405 1.413 0.78 3.631 1.560 2.327 0.88 3.724 1.426 2.611 1.702 4500 0.20 2.935 1.562 1.879 6000 0.33 2.226 1.356 0.60 3.635 1.571 2.313 0.68 3.603 1.409 2.558 1.641 5000 2.33 3.993 2.002 1.994 6000 0.27 2.096 1.334 0.47 3.465 1.516 2.285 0.53 3.605 1.417 2.545 1.571 5000 2.00 3.831 1.928 1.987 7000 2.33 3.912 1.991 1.965 0.38 3.467 1.518 2.284 0.43 3.521 1.398 2.518

Table F .2 continued from previous page

 F Table F.1: Self-diffusion coefficients values of silicon and oxygen atoms on pure SiO 2 systems in the ensemble of simulated conditions. Table F.2: Self-diffusion coefficients values of silicon, oxygen and hydrogen atoms on 9H 2 O.72SiO 2 and 18H 2 O.72SiO 2 systems in the ensemble of simulated conditions.

			Table F.1 continued from previous page
		T	Density 9H 2 O.72SiO 2 D Si	D O	T	Density	D Si 18H 2 O.72SiO 2 D O
	T	(K) Density	(g/cm 3 ) D Si	(m 2 /s) D O	(m 2 /s) D H	(K)	(g/cm 3 ) T Density (m 2 /s) D Si	(m 2 /s) D O	D H
	(K)	(g/cm 3 )	0.46 (m 2 /s)	1.78E-07 1.75E-07 6000 (m 2 /s) (m 2 /s)	(K)	0.67 (g/cm 3 ) 2.00E-07 2.24E-07 (m 2 /s) (m 2 /s)	(m 2 /s)
	T (K) 4000 4000 4000 0.34 2.72 2.21 1.82 1.52 1.28 1.09 0.93 0.81 0.70 0.54 0.42 0.34 2.72 2.21 1.82 1.52 4000 1.28 4000 4000 1.09 F.2 Hydrated systems Density (g/cm 3 ) D Si (m 2 /s) (m 2 /s) D O 2.33 3.39E-09 4.56E-09 5000 T (K) 2.00 5.03E-09 7.07E-09 5500 1.75 5.61E-09 7.03E-09 5500 0.41 1.64E-07 1.84E-07 6000 0.33 2.28E-07 2.33E-07 6000 3000 4.76E-09 5.26E-09 9.24E-09 3000 Density (g/cm 3 ) 0.20 2.33 2.00 0.59 0.52 0.35 5.16E-07 7.00E-07 D Si (m 2 /s) D O (m 2 /s) 3.19E-08 3.98E-08 2.50E-07 2.85E-07 8.56E-09 2.66E-08 1.71E-07 2.50E-07 3.06E-07 0.27 2.96E-07 3.39E-07 6000 0.41 4000 5.71E-09 7.60E-09 3.80E-08 4000 2.81 8.03E-09 1.07E-08 5.30E-08 3.63E-07 4.02E-07 0.20 2.78E-07 3.44E-07 6000 0.33 4000 4.78E-09 6.19E-09 3.24E-08 4000 2.29 5.66E-09 9.50E-09 4.30E-08 4.23E-07 5.59E-07 2.33 2.31E-08 2.61E-08 6000 0.27 4000 5.39E-09 7.35E-09 2.92E-08 4000 1.89 7.18E-09 1.02E-08 4.32E-08 5.78E-07 7.28E-07 2.00 2.23E-08 2.99E-08 7000 2.33 4000 7.69E-09 9.67E-09 4.18E-08 4000 1.57 8.83E-09 1.17E-08 5.13E-08 6.72E-08 8.10E-08 1.75 2.97E-08 4.15E-08 7000 2.00 4000 1.16E-08 1.61E-08 7.31E-08 4000 1.32 1.77E-08 2.17E-08 5.93E-08 7.82E-08 1.07E-07 1.46 4.24E-08 5.26E-08 7000 1.46 4000 1.51E-08 1.87E-08 4.04E-08 4000 1.13 2.22E-08 2.75E-08 6.97E-08 1.15E-07 1.37E-07 1.23 6.12E-08 7.07E-08 7000 1.23 4000 2.12E-08 2.34E-08 6.84E-08 4000 0.97 2.30E-08 3.20E-08 8.54E-08 1.48E-07 1.76E-07 1.00 9.04E-08 8.94E-08 7000 1.00 4000 3.77E-08 4.20E-08 7.22E-08 4000 0.83 3.09E-08 3.76E-08 7.64E-08 1.80E-07 2.03E-07 0.90 9.41E-08 9.97E-08 7000 0.78 4000 3.48E-08 4.18E-08 1.06E-07 4000 0.73 4.61E-08 4.91E-08 1.08E-07 2.76E-07 2.67E-07 0.67 1.71E-07 1.62E-07 7000 0.67 4000 4.75E-08 5.85E-08 1.78E-07 4000 0.56 6.42E-08 9.13E-08 2.83E-07 2.71E-07 3.16E-07 0.52 1.83E-07 2.03E-07 7000 0.59 4000 5.34E-08 8.71E-08 4.34E-07 4000 0.44 7.11E-08 1.18E-07 4.29E-07 3.19E-07 3.67E-07 0.46 2.33E-07 2.55E-07 7000 0.52 4000 1.28E-07 1.45E-07 2.62E-07 4000 0.35 1.20E-07 1.63E-07 4.13E-07 3.55E-07 3.88E-07 0.41 2.19E-07 2.88E-07 7000 0.41 5000 1.95E-08 2.17E-08 1.12E-07 5000 2.81 2.01E-08 2.25E-08 1.18E-07 4.69E-07 4.71E-07 0.33 3.45E-07 4.04E-07 7000 0.33 5000 2.20E-08 2.67E-08 1.06E-07 5000 2.29 2.55E-08 3.24E-08 1.25E-07 6.46E-07 6.93E-07 0.27 4.07E-07 4.58E-07 7000 0.27 5000 3.60E-08 4.10E-08 1.11E-07 5000 1.89 3.04E-08 4.03E-08 1.05E-07 7.66E-07 8.76E-07 5000 4.71E-08 5.32E-08 1.35E-07 5000 1.57 4.58E-08 5.12E-08 1.43E-07 3.57E-08 4.38E-08 1.46 6.74E-09 8.72E-09 5500 1.75 5000 5.92E-08 6.97E-08 1.46E-07 5000 1.32 5.43E-08 6.95E-08 1.63E-07 4.53E-08 5.80E-08 1.23 1.25E-08 1.25E-08 5500 1.46 7.12E-08 7.55E-08 1.00 1.48E-08 1.70E-08 5500 1.23 5000 6.96E-08 9.18E-08 1.58E-07 5000 1.13 6.36E-08 7.22E-08 1.81E-07 5000 0.93 8.62E-08 1.01E-07 2.43E-07 5000 0.97 8.82E-08 1.14E-07 2.17E-07 8.73E-08 9.38E-08 4000 0.90 1.89E-08 2.55E-08 5500 1.00 5000 0.81 9.52E-08 1.17E-07 2.70E-07 5000 0.83 9.46E-08 1.17E-07 2.64E-07 1.09E-07 1.25E-07 4000 0.68 3.32E-08 3.51E-08 5500 0.90 5000 0.70 1.34E-07 1.51E-07 2.88E-07 5000 0.73 1.26E-07 1.37E-07 2.25E-07 1.46E-07 1.54E-07 4000 0.52 3.78E-08 5.81E-08 5500 0.78 5000 0.54 1.70E-07 1.93E-07 3.29E-07 5000 0.56 1.73E-07 1.92E-07 3.58E-07 1.33E-07 1.67E-07 4000 0.46 7.52E-08 1.16E-07 5500 0.67 1.87E-07 2.00E-07 4000 0.41 7.20E-08 8.03E-08 5500 0.59 T Density D Si D O D H T Density D Si D O D H 1.97E-07 2.49E-07 9H 2 O.72SiO 2 5000 0.42 2.20E-07 2.82E-07 3.59E-07 5000 0.44 2.41E-07 3.10E-07 4.94E-07 18H 2 O.72SiO 2 5000 0.34 2.60E-07 3.30E-07 7.19E-07 5000 0.35 2.68E-07 3.39E-07 6.67E-07
	(K)	4000 (g/cm )	0.37 (m 2 /s)	8.43E-08 7.73E-08 5500 (m 2 /s) (m 2 /s)	(K)	0.52 (g/cm 3 ) 2.34E-07 2.59E-07 (m 2 /s) (m 2 /s)	(m 2 /s)
	3000 3000 3000 3000 3000 3000 3000 3000 3000 3000 3000	4000 4000 2.72 4000 2.21 4500 1.82 4500 1.52 4500 1.28 4500 1.09 4500 0.93 4500 0.81 4500 0.70 4500 0.54 4500 0.42	0.33 0.27 5.25E-10 8.52E-10 7.40E-09 1.04E-07 1.54E-07 5500 1.28E-07 1.76E-07 5500 3000 0.46 0.41 0.20 2.03E-07 2.53E-07 5500 0.36 3.37E-10 5.34E-10 3.46E-09 3000 2.33 1.05E-08 1.16E-08 5500 0.33 1.47E-10 3.94E-10 3.66E-09 3000 2.00 1.37E-08 1.79E-08 5500 0.27 2.83E-10 4.36E-10 2.79E-09 3000 1.75 1.59E-08 2.03E-08 5500 0.20 1.95E-10 3.25E-10 1.68E-09 3000 1.46 2.39E-08 2.67E-08 6000 2.33 4.21E-10 6.48E-10 2.94E-09 3000 1.23 3.72E-08 4.30E-08 6000 2.00 3.41E-10 3.22E-09 4.62E-08 3000 1.00 4.48E-08 5.37E-08 6000 1.46 9.90E-10 1.47E-09 4.85E-09 3000 0.90 5.65E-08 6.47E-08 6000 1.23 1.80E-09 1.57E-09 4.04E-09 3000 0.67 8.68E-08 9.29E-08 6000 1.00 1.29E-09 1.59E-09 3.69E-09 3000 0.52 8.95E-08 1.42E-07 6000 0.78 5.20E-09 7.10E-09 3.08E-08 3000	2.74E-07 3.22E-07 2.81 1.17E-09 1.77E-09 1.20E-08 2.99E-07 3.36E-07 2.29 7.21E-10 1.24E-09 8.68E-09 3.07E-07 3.80E-07 1.89 3.60E-10 7.57E-10 5.11E-09 3.56E-07 4.73E-07 1.57 2.54E-11 7.31E-11 1.37E-08 5.55E-07 6.94E-07 1.32 2.14E-10 1.46E-09 1.01E-08 6.09E-07 7.17E-07 1.13 3.65E-10 1.01E-09 5.58E-09 4.18E-08 5.45E-08 0.97 5.85E-10 3.73E-09 2.80E-08 5.18E-08 6.12E-08 0.83 1.98E-09 5.40E-09 3.31E-08 7.45E-08 9.65E-08 0.73 7.27E-10 2.31E-09 1.39E-08 1.12E-07 1.19E-07 0.56 7.68E-10 5.40E-09 4.22E-08 1.33E-07 1.74E-07 1.73E-07 2.21E-07 0.44 3.23E-09 6.34E-09 3.67E-08

Table F .

 F 3: Self-diffusion coefficients values of silicon, oxygen and hydrogen atoms on 36H 2 O.72SiO 2 and 72H 2 O.72SiO 2 systems in the ensemble of simulated conditions.

			36H 2 O.72SiO 2				72H 2 O.72SiO 2
	T	Density	D Si	D O	D H	T	Density	D Si	D O	D H
	(K)	(g/cm 3 )	(m 2 /s)	(m 2 /s)	(m 2 /s)	(K)	(g/cm 3 )	(m 2 /s)	(m 2 /s)	(m 2 /s)
	2000	2.45	1.97E-10 4.63E-10 2.16E-09	2000	2.28	7.35E-10 2.31E-09 8.81E-09
	2000	2.02	8.03E-11 4.06E-10 2.94E-09	2000	1.90	7.96E-10 2.37E-09 7.03E-09
	2000	1.68	1.08E-10 1.56E-09 1.18E-08	2000	1.60	6.63E-10 4.00E-09 1.37E-08
	2000	1.42	3.28E-11 1.45E-09 9.41E-09	2000	1.36	1.04E-09 5.97E-09 2.24E-08
	2000	1.20	9.91E-11 2.35E-09 1.30E-08	2000	1.17	1.91E-09 1.20E-08 3.59E-08

Table F .3 continued from previous page

 F

			36H 2 O.72SiO 2				72H 2 O.72SiO 2
	T	Density	D Si	D O	D H	T	Density	D Si	D O	D H
	(K)	(g/cm 3 )	(m 2 /s)	(m 2 /s)	(m 2 /s)	(K)	(g/cm 3 )	(m 2 /s)	(m 2 /s)	(m 2 /s)
		1.03	2.37E-10 6.62E-09 4.17E-08	2000	1.01	1.84E-09 1.33E-08 4.18E-08
		0.89	2.52E-10 8.96E-09 5.10E-08	2000	0.88	1.35E-09 1.49E-08 5.69E-08
		0.78	6.36E-10 9.12E-09 1.05E-07	2000	0.68	1.49E-09 2.88E-08 9.10E-08
		0.60	4.32E-10 8.75E-09 5.35E-08	2000	0.53	3.41E-09 4.50E-08 1.39E-07
		0.47	5.10E-09 1.86E-08 8.69E-08	2000	0.43	9.62E-09 6.69E-08 2.14E-07
		0.38	6.63E-09 4.09E-08 2.50E-07	3000	2.77	4.73E-09 6.30E-09 3.49E-08
		2.45	2.01E-09 3.60E-09 2.07E-08	3000	2.28	3.40E-09 6.15E-09 2.18E-08
		2.02	9.40E-10 2.32E-09 1.18E-08	3000	1.90	3.15E-09 8.17E-09 2.97E-08
		1.68	6.63E-10 1.90E-09 8.46E-09	3000	1.60	4.69E-09 1.06E-08 2.62E-08
		1.42	1.52E-09 2.83E-09 9.29E-09	3000	1.36	2.34E-09 1.21E-08 3.45E-08
		1.20	1.21E-09 4.24E-09 3.07E-08	3000	1.17	3.89E-09 1.60E-08 6.78E-08
		1.03	1.51E-09 5.38E-09 3.55E-08	3000	1.01	5.61E-09 2.16E-08 8.08E-08
		0.89	1.72E-09 1.28E-08 7.76E-08	3000	0.88	2.76E-09 2.61E-08 8.54E-08
		0.78	1.79E-09 9.56E-09 7.82E-08	3000	0.68	1.62E-08 5.69E-08 1.49E-07
		0.60	3.14E-09 1.40E-08 1.21E-07	3000	0.53	1.73E-08 8.23E-08 2.53E-07
		0.47	2.08E-08 5.08E-08 2.07E-07	3000	0.43	5.06E-08 1.11E-07 3.10E-07
		0.38	2.17E-08 5.71E-08 2.96E-07	4000	2.77	1.21E-08 1.57E-08 6.52E-08
		2.45	8.10E-09 1.14E-08 4.79E-08	4000	2.28	1.42E-08 2.01E-08 5.93E-08
		2.02	9.61E-09 1.40E-08 4.48E-08	4000	1.90	1.35E-08 2.05E-08 6.44E-08
		1.68	1.16E-08 1.80E-08 5.14E-08	4000	1.60	2.53E-08 3.29E-08 7.98E-08
		1.42	1.55E-08 2.51E-08 6.98E-08	4000	1.36	2.41E-08 3.69E-08 8.60E-08
		1.20	2.11E-08 3.09E-08 8.70E-08	4000	1.17	3.40E-08 5.19E-08 1.11E-07
		1.03	2.97E-08 4.01E-08 1.16E-07	4000	1.01	4.46E-08 7.09E-08 1.30E-07
		0.89	3.15E-08 5.75E-08 1.76E-07	4000	0.88	4.96E-08 9.06E-08 1.70E-07
		0.78	5.36E-08 7.75E-08 1.73E-07	4000	0.68	7.97E-08 1.29E-07 2.58E-07
		0.60	7.45E-08 1.15E-07 2.68E-07	4000	0.53	8.98E-08 1.65E-07 3.41E-07
		0.47	1.06E-07 1.60E-07 3.79E-07	4000	0.43	1.05E-07 2.19E-07 4.21E-07
		0.38	1.27E-07 1.77E-07 4.21E-07

Acknowledgements

I'm grateful to the École Normale Supérieure de

List of Tables Appendix A

Convergence tests

A.1 ENCUT

The cutoff energy was tested at the lowest and highest temperature of interest, meaning at 4000 K and 7000 K. The choice of the density was made based on the same assumption, where the conditions tested were the highest and lowest density value at their specific temperature, meaning (2.33, 0.27 and 0.20 g/cm 3 ). The energy value was normalized by the number of atoms (216) and by the energy of the highest cutoff energy value. The variation observed in the energy is of less than 1 meV/atom, which lead to the choice of a small cutoff energy (550 eV) in order to increase the speed of the calculations. Circles represent the simulation at the chosen cutoff energy. Lattice parameters correspond to 2.33 g/cm 3 (14.554 Å), 0.20 g/cm 3 (a33.000 Å) and 0.27 g/cm 3 (a30.000 Å). The solid line represents the approximate ENCUT value chosen as a parameter in the scope of simulations.

Appendix C Pair distribution function analysis

From pair distribution function analysis of dry silica simulations, we extract mean values of first peaks and troughs for Si-O and O-O pairs, as detqiled in Table C below. Figure E.1: Mean square displacement of (a) silicon and (b) oxygen on pure silica simulations at all temperatures and all densities. Different colors represent simulations at different temperatures, where dark blue indicate calculations at 4000 K, light blue is 4500 K, green is 5000 K, yellow is 5500 K, orange is 6000 K and red is 7000 K results.

E.2 Hydrated systems