Iron as a building block material of the Earth naturally received significant attention. Considerable efforts have been made to determine its thermodynamic and thermophysical properties up to the Earth's inner core's conditions. However, its physical properties in the low-density regime are less explored, and notably the position of the liquid-vapor equilibrium line and of the critical point are lacking. The missing information inhibits developing a complete equation of state that covers the released state after shock waves, and thus hinders the characterization of large planetary impacts.

The present study aims at closing the knowledge gap on the liquid-vapor equilibrium dome of iron. For this we exploit molecular dynamics and Monte Carlo methods where the energy and the forces are estimated by the density functional theory. We then employ statistical and thermodynamics methods to construct the position of the critical point, build the liquid-vapor dome, and characterize the physical properties of the fluid iron.

First we determine the position of the critical point from ab initio molecular dynamics simulations along several isotherms. The simulation results provide the position of the liquid spinodal above 3000 K, and the gas spinodal close to the critical point. We bracket the position of the critical point in the 9000-9350 K temperature range, and 1.85-2.40 g/cm 3 density range, corresponding to 4-7 kbars pressure range. Additionally, we characterize the structure and the transport properties of the fluid iron over a wide density and temperature range, with a particular focus on the supercritical state.

Then we compute two Hugoniot lines starting with two realistic initial conditions. By comparing the entropy values calculated along these Hugoniot lines to that at the boiling point, we find that the pressure required to reach the onset vaporization is significantly lower than previous estimates. It suggests that previous hydrodynamic simulations underestimate the iron vapor production, and that the core of Theia underwent partial vaporization during the giant impact. Similarly, we find that a large fraction of the planetesimals falling on Earth during the late veneer must have had their cores undergoing partial vaporization. The readily achieved partial core vaporization would enhance the iron-silicates equilibration, which helps explain geochemical observations. At last, we determine the liquid-vapor equilibrium line of iron. For this, we have extended and implemented the Gibbs ensemble Monte Carlo method coupled with the finite-temperature density functional theory. The first benchmark test to sodium shows a good agreement with available experimental results. We then apply this technique to iron and calculate its liquid density in equilibrium with the vapor phase. We also show the importance of magnetism diminishes as approaching the critical point.

Résumé

Le fer, en tant que brique élémentaire de la Terre, a reçu beaucoup d'attentions. Des efforts considérables ont été mis en oeuvre pour déterminer ses propriétés thermodynamiques et thermophysiques à des conditions atteignant celles du noyau terrestre. Cependant, ses propriétés physiques dans le domaine des faibles densités sont moins explorées, et il manque en particulier la position de la courbe d'équilibre liquidegaz et du point critique. Les informations manquantes entravent le développement d'une équation d'états complète qui couvrirait l'état de détente post onde de choc, et donc empêchent la caractérisation des grands impacts planétaires.

Cette étude vise à réduire le fossé de connaissances sur l'équilibre liquide-gaz du fer. Pour cela nous utilisons la Dynamique Moléculaire et la méthode Monte Carlo dans lesquelles les énergies et les forces sont estimées à partir de la théorie de la fonctionnelle densité. Nous utilisons ensuite des méthodes statistiques et thermodynamiques pour construire la position du point critique, le dôme liquide-gaz, et caractériser les propriétés physiques du fer à l'état de fluide.

Tout d'abord nous avons déterminé la position du point critique à partir de simulations de dynamique moléculaire ab initio selon plusieurs isothermes. Les résultats des simulations nous ont donné la position du spinodal liquide au-dessus de 3000 K, et du spinodal gazeux à proximité du point critique. La position du point critique est estimée entre 9000-9350 K et 1.85-2.4 g/cm 3 , ce qui correspond à 4-7 kbars. Nous avons également caractérisé la structure et les propriétés de transport du fer fluide pour une large gamme de densités et températures, avec une attention particulière sur l'état supercritique.

Ensuite nous avons calculé deux courbes Hugoniot à partir de deux conditions initiales réalistes. En comparant les valeurs d'entropie calculées le long de ces courbes à celle du point d'ébullition, nous avons trouvé que la pression requise pour atteindre le seuil de vaporisation est significativement plus basse que précédemment estimée. Cela suggère que les simulations hydrodynamiques précédentes sous-estiment la production de vapeur de fer, et que le noyau de Théïa aurait subi une vaporisation partielle lors de l'impact géant. De même nous avons trouvé qu'une grande fraction des planétésimaux ayant frappé la Terre lors du vernis tardif ont dû voir leur noyau vaporisé partiellement. La facilité avec laquelle les noyaux se vaporisent devrait améliorer l'équilibration fer-silicate, ce qui permettrait d'expliquer les observations géochimiques.

Enfin, nous avons déterminé l'équilibre liquide-gaz du fer. Pour cela nous avons amélioré et implémenté la méthode Monte Carlo dans l'ensemble de Gibbs couplée avec la théorie de la fonctionnelle densité en températures finies. Le premier test de référence avec le sodium nous as donné un bon accord avec les résultats expérimentaux. Nous avons donc appliqué cette technique au fer et calculé sa densité liquide à l'équilibre avec la phase vapeur. Nous avons également montré que l'importance du magnétisme diminue à l'approche du point critique.

Introduction

The Moon is the most well-studied satellite in our solar system, whose investigation can be dated back to the epoch of Hipparchus and Ptolemy. In the ancient time, the primary attempt was to observe the lunar orbital motion, document its moving path and estimate its speed relative to a fixed star. These tedious and time-consuming astrophysical observations leads to the discovery of several exciting phenomena like the lunar precession. After Isaac Newton formulated the laws of motion, which laid the foundation of celestial mechanics, the focus has been shifted to explain the observed lunar motion from the perspective of mutual gravitational attractions among the Earth-Moon-Sun system. The modern research into the Moon was mostly motivated by the Apollo missions. These missions returned tons of the lunar rocks. From measurements of their chemical composition, our understanding of the lunar interior structure has deepened considerably. These new findings, in return, stimulate scientists from different disciplines to research into the possible origin of the Moon. This chapter will focus on the giant impact theory, which is considered as the most plausible mechanism to explain the Moon's formation. We will first review the major geochemical observations that the giant impact theory must reconcile. We then point out the total angular momentum at the time of the giant impact is still under debate. Possible values are anywhere between once and twice the present value of the Earth-Moon system. Based on this, we can divide the giant impact theory into two categories: the canonical impact model and the high-angular momentum impact model. Next, we give a brief synopsis of the geochemical consequence for these two models. With these background information in mind, we narrow down to the topic of iron, which is the main component of planets' core, and explore its behavior during giant impacts. By examining the phase diagram of iron, we identify a knowledge gap on the thermodynamic properties of iron in the low-density region, which would inhibit developing an accurate equation of state. As an essential ingredient in the giant impact modeling, the equation of state affects the amount of vapor produced as well as the energy distribution across the Earth during the impact. This chapter will end up with research problems on which the present study will focus.

The possible origin scenarios for the Moon

Many theories have been devised to explain the Moon's origin, including the fission theory, the co-formation theory, the capture theory, and the giant impact hypothesis. In the fission theory, the proto-Earth rotates so rapidly that it becomes dynamically unstable. A lunar-sized mass from the Earth's mantle is separated and injected into the Earth's orbit to be the satellite. Therefore, it is expected that the Moon has a similar isotopic composition and is depleted in iron relative to the Earth. In the co-formation theory, the Moon forms at the same time as its parent Earth. Under such an explanation, in the early solar system, gravity would pull materials together to form the Moon simultaneously as gravity bonds particles together to form the Earth. Such a moon would be expected to have a very similar isotope composition to the Earth. In the capture theory, a rocky body formed elsewhere in the solar system and was later captured by the Earth's gravitational field. Since the body captured has a wide range of possibilities in its isotope composition and iron content, various chemical consequences are possible. However, each theory stated above fails to explain one of the major characteristics of the Earth-Moon system (see [START_REF] Stevenson | Origin of the moon-the collision hypothesis[END_REF], for a review). The fission theory needs three times more angular momentum than the present Earth-Moon system. But there is no mechanism available to supply such large angular momentum during the planetary accumulation process. The co-formation theory is also challenging due to the difficulties in explaining the lunar iron depletion relative to the Earth (see Section 1.3). It is more reasonable to expect a very similar iron content if both bodies form together. In the capture theory, the likelihood of capturing a lunar-sized body to the present lunar orbit is nearly zero.

Influenced by Safronov's work on the planet formation through collisions between planetary embryos, both [START_REF] Hartmann | Satellite-sized planetesimals and lunar origin[END_REF] and [START_REF] Cameron | The origin of the moon[END_REF] proposed the giant impact hypothesis as a plausible explanation for the formation of the Moon. In this model, appropriate impact conditions satisfying the Earth-Moon system's angular momentum constraints would generate a proto-lunar disk with most of the materials originating from the silicate mantle of the impactor or of the proto-Earth, then a Moon is accreted from such a disk. Therefore, the hypothesis can easily explain the lunar iron depletion relative to the Earth. In addition, the giant impact would happen naturally since the collision between planets is thought to be common in the late accretion stage. These characters have rendered the giant impact hypothesis the leading theory to explain the Moon's origin.

The geochemical constraints

Measuring the lunar rocks' chemical composition returned by Apollo missions and their comparisons to the Earth provide important constraints on which the giant impact hypothesis must reconcile. Both the Moon and the Earth have a very similar interior structure, with an iron core in the center and a silicate mantle lying above. But the average density is different with 3.34 g/cm 3 for the Moon compared to 5.51 g/cm 3 for the Earth, indicating a relative lunar iron depletion. There are two potential iron reservoirs in the Moon: the silicate mantle/crust and the metallic core. Although it is debated on the exact size of the Moon's core, it is generally accepted that its size is at most 400 km in radius (see Canup, 2004b, for a review). The combination of seismic analysis and geological modeling sets an upper bound of the total amount of iron, which is around 8-10 wt% in the Moon compared with 30 wt% in the Earth.

Despite the significant difference in the total amount of iron, a remarkable similarity in the oxygen isotope composition between the Moon and the Earth's silicate mantle was revealed by measuring the oxygen isotope composition of the lunar rock samples. Initially, it was considered that the oxygen isotopes were homogeneous in the solar system, which can explain the observed similarity in the oxygen isotope between the Earth and Moon. However, [START_REF] Franchi | The oxygen-isotopic composition of earth and mars[END_REF] found an offset in fractionation line between the Earth and Mars, indicating oxygen isotopes were distributed heterogeneously in the early solar system. Subsequent N-body simulations also suggest the possibility of a similar oxygen isotope composition between the impactor and the proto-Earth is very low [START_REF] Pahlevan | Equilibration in the aftermath of the lunarforming giant impact[END_REF]. As the oxygen isotope remains unchanged during subsequent geological processes after the Moon's formation, the isotopic indistinguishability provides a strong argument that the Moon's origin is closely linked with the Earth's mantle, and thus provide a strong constraint on developing the giant impact theory.

The Earth and the Moon also show indifference in titanium isotope composition ( 50 Ti/ 47 Ti) [START_REF] Zhang | The proto-earth as a significant source of lunar material[END_REF]. Since titanium was distributed heterogeneously in the solar system, the impactor should differ its titanium isotope composition to the proto-Earth. Then we would expect the Moon to develop an offset in 50 Ti/ 47 Ti compared to the Earth, but they do not. The titanium-isotope homogeneity in the Earth-Moon system further confirms that the Moon is accreted from the proto-lunar disk, where most materials are from the Earth's mantle.

Two studies have recently measured the tungsten isotopic composition and found the Moon has a slightly higher 182 W/ 184 W ratio than the Earth [START_REF] Kruijer | Lunar tungsten isotopic evidence for the late veneer[END_REF][START_REF] Touboul | Tungsten isotopic evidence for disproportional late accretion to the earth and moon[END_REF]. Since the tungsten isotopes vary throughout the solar system, it is reasonable to expect that the giant impact event would have involved two objects with different isotopic compositions of tungsten. In order to explain the slight difference in 182 W/ 184 W, both studies suggest that the Earth and Moon had the same tungsten isotopic composition after the giant impact, and the difference in the 182 W/ 184 W ratio between the Earth and Moon is due to the addition of the late veneer, which is thought to have a chondritic composition. The model proposed also requires siderophile elements stripping from the mantle by metal from the impactor's core. Otherwise, the tungsten data would be too abundant due to the addition of the late veneer, which would not be self-consistent with the highly siderophile element concentration in the terrestrial mantle.

Several volatile elements like K and Zn are found to differ their isotopic composition in the Moon and Earth. [START_REF] Wang | Potassium isotopic evidence for a high-energy giant impact origin of the moon[END_REF] found the lunar rocks are slightly enriched in the heavy isotope of K compared to the Earth. They suggest it may result from the incomplete condensation from a bulk silicate Earth vapor. During the condensation, the lighter isotope would prefer to stay in the vapor phase, causing the enrichment of heaver isotope in the solid or liquid phases. [START_REF] Wang | Potassium isotopic evidence for a high-energy giant impact origin of the moon[END_REF] also indicate the potassium isotope result is inconsistent with the canonical impact model but supports the high-angular momentum giant impact model for the origin of the Moon. However, as discussed by [START_REF] Canup | Lunar volatile depletion due to incomplete accretion within an impact-generated disk[END_REF], the preferential accretion of volatile-rich melt in the inner disk to the Earth, rather than to the Moon, could also contribute to the volatile elements depletion and their isotope fractionation in the Moon. Therefore, the enrichment of the heavy K isotope in the Moon may only support the presence of a proto-lunar disk from which the Moon is accreted, but cannot be used to distinguish different impact models. [START_REF] Paniello | Zinc isotopic evidence for the origin of the moon[END_REF] found lunar rocks are enriched in the heavy isotope of Zn and have a lower Zn concentration than the Earth. Since the equilibrium condensation would not produce such a large degree of fractionation, they suggest these variations are due to the lunar magma ocean's large-scale evaporation, which not only decreases the concentration of the volatile elements but also causes the fractionation.

In summary, the Moon is strongly depleted in iron but has a similar refractory element isotope composition to the Earth. Besides, the Moon is enriched in the heavier isotope of the volatile elements, which is caused by the evaporation process after the Moon's formation and/or the equilibrium condensation accompanying the Moon's formation.

1. 4 The dynamical constraint from the lunar orbital evolution

To find a set of impact parameters that can produce a massive and iron-depleted proto-lunar disk, we also need the dynamic constraint of the total angular momentum, which includes the rotational angular momentum of the Earth and the Moon, and the orbital angular momentum of the lunar motion. We know that the Moon is almost in a circular cycle revolving around the Earth with an eccentricity of 0.0549 and a slight inclination of 5.14 • to the ecliptic plane (Fig. 1.1). From these real measurements, we can calculate the current total angular momentum of the Earth-Moon system. However, this value may not be the total angular momentum at the time when the Moon formed. Therefore, we need to find a way to recover its evolving history and obtain the total angular momentum at the time when it completed the accretion and detached from the proto-lunar disk. This is not an easy task since the Moon is continuously affected by the tidal interaction exerted by the Earth and Sun to recess its orbit.

The lunar orbital evolution

We can partially solve this problem by celestial mechanics and an essential contribution has been made by [START_REF] Touma | Evolution of the earth-moon system[END_REF] and [START_REF] Touma | Resonances in the early evolution of the earth-moon system[END_REF]. As these calculations are extremely complex due to the employment of generalized coordinates like action-angle variables, only a general frame is outlined. Starting with the derivation of an analytic Hamiltonian for the Earth-Moon-Sun system, [START_REF] Touma | Evolution of the earth-moon system[END_REF] made several assumptions like the lunar orbit is a circle to reduce the problem into a manageable size. Together with averaging over the short orbital time scale, the remaining degrees of freedom in the Hamiltonian is the obliquity of the Earth to the ecliptic plane, the mutual inclination of the lunar orbit to the Earth's equator, the inclination of lunar orbit to the ecliptic plane, and the precession rate of the lunar orbit. These variables are constantly changing due to the Earth's oblateness and tidal interactions. Although the physical origin of the tidal torque can easily be understood as the presence of a phase lag due to the mismatch between the orbital and rotational motion, for practical implementations, we need to employ different tidal models to approximate it as a function of the orbital elements of the Earth and the Moon. Fortunately, the detailed evolving history of the lunar orbit shows indifference to a particular model. The most striking finding is that the mutual inclination has a value of 12 • when approaching the Earth [START_REF] Goldreich | History of the lunar orbit[END_REF][START_REF] Touma | Evolution of the earth-moon system[END_REF]. However, it disagrees with the giant impact theory which predicts the formation of the Moon is near to the Earth's equator and thus the initial value is close to zero.

1.4.2 The high mutual inclination and the initial total angular momentum

The high mutual inclination problem has been re-investigated by [START_REF] Touma | Resonances in the early evolution of the earth-moon system[END_REF]. With direct numerical simulations, they discover two new phenomena called evection and eviction resonance that have a dramatic impact on the dynamic evolution of the Earth-Moon system. Resonance is a phenomenon of describing increased oscillation amplitude when the frequency of a periodically applied force matches the natural frequency of the system. We can exemplify this interesting phenomenon by the vibration of a harmonic oscillator. If the driving frequency is equal to the natural frequency of the oscillator, the resonance is achieved and the amplitude of the oscillation increases dramatically. This concept can be extended into the planetary science where resonance means a commensurability amongst the frequencies of precession or orbital motions between two planets or stars. If it happens, a periodic gravitational influence will excite the orbital eccentricity or inclination to high values just like the vibration amplitude of an oscillator. The continuous action from the Earth's tidal torque will make the lunar orbital expand and the precession period of the perigee increase (see Figure 13 in [START_REF] Touma | Resonances in the early evolution of the earth-moon system[END_REF]. When this period is close to the orbital period (about one year) of the Earth, the evection resonance happens, and pumps up the eccentricity of the lunar orbit until the tides induced by the Moon and by the Earth cancel out so that the Moon stops moving outward. Then the lunar tide becomes slightly stronger than the tide in the Earth causing the contraction of the lunar orbit. During this stage, the spin rate of the Earth declines in a steady manner until escaping from the evection resonance. After that, the Moon encounters a mixed evection-inclination resonance twice which excites the mutual inclination to 12 • . [START_REF] Touma | Resonances in the early evolution of the earth-moon system[END_REF] found the Moon escapes the evection resonance very fast and only a tiny amount of angular momentum is removed from the Earth-Moon system. However, [START_REF] Ćuk | Making the moon from a fast-spinning earth: a giant impact followed by resonant despinning[END_REF] suggest the contraction period is prolonged causing a significant de-spinning of the Earth and efficient removal of the angular momentum to the heliocentric orbit of the Earth. They suggest the initial total angular momentum at the time of the giant impact can be twice as large as the current value. By employing a different tidal model, [START_REF] Wisdom | Early evolution of the earth-moon system with a fastspinning earth[END_REF] and [START_REF] Ward | Analytical model for the tidal evolution of the evection resonance and the timing of resonance escape[END_REF] found too much or too little angular momentum was removed from the Earth-Moon system if the initial angular momentum is high, which is inconsistent with the present Earth-Moon system. Different results obtained from applying different tidal models suggest more work is needed to investigate whether the evection resonance is a viable mechanism to remove the excess angular momentum. Ćuk et al. (2016) proposed a new origin scenario for the Moon. In this model, the Moon is accreted from the disk generated by a high-angular momentum impact with an initial high-obliquity (70 The evolution of the mutual inclination as a function of radius. The mutual inclination must be at least 12 • when the Moon forms. However, the giant impact theory predicts the moon formed at the Earth's equator predicting an initial value of 0 • . This discrepancy can be reconciled with the Moon-Sun resonance [START_REF] Touma | Resonances in the early evolution of the earth-moon system[END_REF], the disk-satellite resonance [START_REF] Ward | Origin of the moon's orbital inclination from resonant disk interactions[END_REF] or the Laplace plane transition (Ćuk et al., 2016), which can excite the mutual inclination to a high value (see the black dashed line in b). This figure is reproduced with permission from [START_REF] Touma | Evolution of the earth-moon system[END_REF]. (c) [START_REF] Ćuk | Making the moon from a fast-spinning earth: a giant impact followed by resonant despinning[END_REF] have shown the Earth's rotational rate slows down due to the Moon-Sun resonance which can be used to drain away the excess angular momentum. They suggest an initial high angular momentum at the time of the giant impact is dynamically feasible. However, subsequent studies using different tidal models found too much [START_REF] Wisdom | Early evolution of the earth-moon system with a fastspinning earth[END_REF] or too little [START_REF] Ward | Analytical model for the tidal evolution of the evection resonance and the timing of resonance escape[END_REF] angular momentum be removed from the Earth-Moon system, raising the question whether the Moon-Sun resonance is still a viable mechanism. This figure is reproduced with permission from [START_REF] Ćuk | Making the moon from a fast-spinning earth: a giant impact followed by resonant despinning[END_REF].

transition will excite the lunar inclination to 30 • and reduce the obliquity of the Earth to 20 • . This process will also remove the excess angular momentum from the Earth-Moon system and transfer to the heliocentric orbit of the Earth. Once the Moon has passed through the Laplace phase transition, it will undergo the Cassini state transition to achieve the spin-orbit resonance causing a large lunar obliquity (over 30 • ). The high-obliquity tides in the Moon, in return, will strongly damp the lunar inclination to reach the present value of 5 • . However, [START_REF] Tian | Vertical angular momentum constraint on lunar formation and orbital history[END_REF] found the vertical component of the total angular momentum in the Earth-Moon system, which is perpendicular to the Earth's orbital plane, is almost conserved.

They suggest an initial high-obliquity (70 • ) Earth with a high angular momentum proposed in Ćuk et al. (2016) would result in too large vertical angular momentum that is inconsistent with the present observations.

Summary

The total angular momentum at the time the Moon formed is an essential part of constraining the impact conditions. If there is no mechanism to remove or add angular momentum to the Earth-Moon system, its initial angular momentum is the same as the present Earth-Moon system. [START_REF] Ćuk | Making the moon from a fast-spinning earth: a giant impact followed by resonant despinning[END_REF] proposed the Moon-Sun resonance would remove a large amount of angular momentum from the Earth-Moon system. Therefore, the initial angular momentum could be up to two times higher than the present observation. However, [START_REF] Wisdom | Early evolution of the earth-moon system with a fastspinning earth[END_REF] and [START_REF] Ward | Analytical model for the tidal evolution of the evection resonance and the timing of resonance escape[END_REF] found the resonance appears to remove too much to too little angular momentum if the initial angular momentum is high, causing inconsistency with the present Earth-Moon system. Ćuk et al. (2016) suggest that the instability associated with the Laplace plane transition would remove considerable angular momentum from the Earth-Moon system. However, [START_REF] Tian | Vertical angular momentum constraint on lunar formation and orbital history[END_REF] suggest an initial high obliquity Earth with a high angular momentum proposed in Ćuk et al. (2016) cannot produce the present Earth-Moon system.

1.5 The giant impact theory

Material injection mechanism

Before we are going to discuss in detail the impact process and its geochemical implications, a concise introduction of the material injection mechanism is given. The energetic collision between the impactor and the proto-Earth creates powerful shock waves, which will compress constitutive materials very rapidly to a high pressure and temperature condition. When the shock wave reaches the free space assumed to be a vacuum, to sustain the zero pressure interface a rarefaction wave or a relief wave must be produced. Then it reflects from the free space and travels into the compressed materials to make them expand [START_REF] Forbes | Shock Wave Compression of Condensed Matter[END_REF]. As the pressure drops and volume rises during the expansion, the particle velocity of shocked materials would increase, so is the total energy which is defined as the sum of the kinetic In the canonical model (a), the graze collision of a Mars-sized projectile with proto-Earth at a velocity of 10 km/s can generate a disk (b) and leave the total angular momentum of post-impact structure close to the present Earth-Moon system. The disk material is made of silicates liquids with 20 wt% vapor and primarily derived from the impactor, from which the Moon is accreted. Thus the Moon should have a distinct isotope composition compared to the Earth [START_REF] Canup | Origin of the moon in a giant impact near the end of the earth's formation[END_REF]. [START_REF] Canup | Accretion of the moon from an impact-generated disk[END_REF] suggest a large amount of mass must be injected directly beyond the Roche limit to form a single lunar-mass satellite. In (b), all materials outside of the Roche limit are represented by a single Moon as the accretion process is very fast on a timescale of years. Due to the Jeans instability the high density magma disk inside the Roche limit will clump and then be sheared apart by the tidal force, which results in an effective viscosity and is much larger than the molecular viscosity [START_REF] Stevenson | Origin of the moon-the collision hypothesis[END_REF]. The entropy file in (e) represented by the blue dashed line is typical for the canonical impact, which indicates the entire structure is thermally stratified. In the high-angular momentum impact model (d), a smaller impactor with a high velocity of 30 km/s hitting a rapid spinning proto-Earth can produce a disk with enough material mixing between colliding bodies to account for the measured isotopic ratios [START_REF] Ćuk | Making the moon from a fast-spinning earth: a giant impact followed by resonant despinning[END_REF]. The post-impact structure is highly thermally stratified (red line in (e)) and the silicate transitions smoothly from liquid to the supercritical fluid to vapor [START_REF] Lock | The origin of the moon within a terrestrial synestia[END_REF], which was named as synestia. Besides, there is a smooth change in the angular velocity from the corotating inner region to the sub-Keplerian disk-like region (f). energy and gravitational potential energy. For the high-angular momentum impact model, we need to include the rotational energy from the pre-impact spin to the total energy as well, which will make the injection more easily.

If the total energy of the injected materials is negative, they will follow a Keplerian orbit with the periapse on the Earth if there are no other materials block their pathway to inject. After one orbital period they will fall back and re-impact with the Earth. In contrast, if the injected materials have a positive total energy, they will escape the Earth's gravitational field. In order to form a proto-lunar disk, we need other mechanisms to lift the periapse of some parts of shocked materials with a negative total energy above the Earth.

One possible mechanism is the gravitational torque [START_REF] Stevenson | Origin of the moon-the collision hypothesis[END_REF], which pumps the angular momentum into the injected materials to make their orbits lift and avoid to re-impact with the Earth. This mechanism is very similar to the current recession of the lunar orbit due to the tidal torque exerted by the Earth. Another mechanism is the pressure gradient, which becomes important if there is significant vaporization happening during the giant impact. As compression is an irreversible process, the shocked materials will gain entropy in the course to reach the peak pressure and temperature conditions. Then an expansion is followed due to relief waves. As the decompression process happens very fast, convection, viscous dissipation and radiation play a very limited role. We can safely assume this process is isentropic. If the decompressed materials intersect with the liquid-vapor dome, it will turn into the liquid-vapor mixture. The significant volume change from the condensed phase (liquid or solid) to the vapor phase upon vaporization causes an abrupt increase in the particle velocity. Then the outflow materials is subjected to the pressure gradients which is able to lift the periapse of injected materials (see Figure 6 in [START_REF] Stevenson | Origin of the moon-the collision hypothesis[END_REF]. Due to the complex impact geometry, how much materials are injected into the proto-lunar disk and which mechanism is dominant can only be determined from hydrodynamic simulations.

The canonical impact model

As the laboratory-scale experiments are not able to simulate such planetary-scale impacts, our understanding of the giant impact mostly comes from hydrodynamic simulations (Fig. 1.2). In this section, I summarize the outcome for the canonical impact modeling. The high angular momentum impact scenario will be analysed in the next section. For simulation results, we are more interested in what kind of impact parameters will produce the present Earth-Moon system, which includes the total mass (M T ) as a sum of the impactor and the proto-Earth, the impact-tototal-mass ratio (γ), impact angle (b) and the total angular momentum (L imp ). The simulation details are well beyond my interests.

Benz and Cameron pioneered the application of the SPH method to the giant impact simulations. Interestingly, a series of papers were published in Icarus with the same title 'The origin of the Moon and the Single Impact Hypothesis' [START_REF] Benz | The origin of the moon and the singleimpact hypothesis I[END_REF][START_REF] Benz | The origin of the moon and the singleimpact hypothesis II[END_REF][START_REF] Benz | The origin of the moon and the singleimpact hypothesis III[END_REF][START_REF] Cameron | The origin of the moon and the single impact hypothesis V[END_REF][START_REF] Cameron | The origin of the moon and the single impact hypothesis IV[END_REF]. Due to the limit on the computational speed, all simulations except for in [START_REF] Cameron | The origin of the moon and the single impact hypothesis V[END_REF] have a low resolution. For instance, the total number of particles in the generated proto-lunar disk is only 30 [START_REF] Benz | The origin of the moon and the singleimpact hypothesis II[END_REF]. From these simulations, they conclude that 1. the low-mass impactors with the impact-to-total-mass ratio less than 0.12 produces iron-rich disk;

2. the gravitational torque is more important than the pressure gradient for materials emplacement, which still holds even in the present high-resolution hydrodynamic simulations. [START_REF] Cameron | The origin of the moon and the single impact hypothesis V[END_REF] first performed the high-resolution simulations with a total number of particles around N ∼ 10 4 . The most successful impact to produce the current Earth-Moon system is with γ=0.3 and M T =0.65M E , where M E is the mass of the Earth. This case is also called the early-Earth impact scenario and the Earth needs to acquire extra 0.35M E by the late accretion. Since the Moon will receive a proportional amount of material as well and there is no mechanism available to filter out iron in these materials, the initial iron-depleted Moon will become ironrich again causing an inconsistency with the geochemical observations. Therefore, the early-Earth impact scenario is not favorable. Canup et al. (2001) re-examined the simulation results in [START_REF] Cameron | The origin of the moon and the single impact hypothesis V[END_REF] and proposed a scaling law to describe the results of satellite-forming impact simulations. They found the disk mass tends to increase and iron content decreases with increasing b for 0.4 < b < 0.8, which is independent of M T . The maximum yielding of the massive and iron-depleted disk is at b ∼ 0.8. [START_REF] Canup | Origin of the moon in a giant impact near the end of the earth's formation[END_REF] were guided by this trend and revisited the small impactor case with γ < 0.12. A total of 36 impact simulations were run with γ = 0.108 -0.115, b = 0.70 -1.0, L imp = L EM and M T = M E , where L EM is the total angular momentum of the present Earth-Moon system. They found a massive proto-lunar disk can be generated to allow the accretion of a single moon. They also confirm most of the materials in the lunar disk is from the impactor with a mass fraction of 0.6-0.74, indicating that the Moon will inevitably have a different oxygen isotope ratio compared to the Earth's mantle.

1.5.3 The high-angular momentum impact model [START_REF] Ćuk | Making the moon from a fast-spinning earth: a giant impact followed by resonant despinning[END_REF] have proposed a new giant impact scenario which can produce a Moon being isotopically similar to the Earth. In this model, a small impactor with a mass around 0.026-0.1 M E but with a high velocity (30 km/s) hits a rapidly spinning Earth. It results in a vapor-dominated disk with enough material mixing between colliding bodies to account for the measured isotopic ratios. The resulted disk has a much high angular momentum and more massive than the one generated by the canonical impact model. However, the successful Moon-forming impact leaves the Earth-Moon system with an excess angular momentum. [START_REF] Ćuk | Making the moon from a fast-spinning earth: a giant impact followed by resonant despinning[END_REF] propose the Earth-Moon system can lose angular momentum by the orbital resonance between the Sun and Moon, which is still under debate. [START_REF] Lock | The origin of the moon within a terrestrial synestia[END_REF] found the post-impact structure is highly thermally stratified and the silicate transitions smoothly from the liquid phase to the supercritical phase, then to the vapor phase. This special structure was named as synestia. There is also a smooth change in the angular velocity from the corotating inner region to the sub-Keplerian disk-like outer region. The outer part of the disk-like region is likely to be well mixed due to the falling condensates and vertical fluid motion [START_REF] Lock | The origin of the moon within a terrestrial synestia[END_REF]. However, a whole mixing in the synestia may be difficult. On the one hand, synestia is thermally stratified, meaning the outer part is hot and has a lower density. Lifting denser materials from the inner part to the outer part needs to overcome the gravitational force. On the other hand, the monotonic increase of the specific angular momentum from the inside-out creates a barrier. The exchange of a large amount of materials between the inner part and outer part in the synestia will decrease the angular momentum of the outer part, result in the collapse of the disk-like region and leave too little materials in the disk [START_REF] Melosh | New approaches to the moon's isotopic crisis[END_REF].

Summary

The giant impact model is the leading theory to explain the formation of the Moon.

The classic canonical impact model fails to explain the Moon's isotopic similarity to the Earth, while the high-angular momentum impact model has problems with the removal of the excess angular momentum to match the present Earth-Moon system.

The behaviour of iron during giant impacts

Both the canonical and high-angular momentum impact models predict the impactor's core merge rapidly into the proto-Earth's core at a timescale of hours [START_REF] Kraus | Impact vaporization of planetesimal cores in the late stages of planet formation[END_REF]. Since the chemical exchange between the sinking iron from the impactor and the silicates in the magma ocean requires a much longer time, there is a very limited chemical equilibration between the impactor's core and the Earth's silicate mantle. However, the hafnium-tungsten isotope studies [START_REF] Kruijer | Lunar tungsten isotopic evidence for the late veneer[END_REF][START_REF] Touboul | Tungsten isotopic evidence for disproportional late accretion to the earth and moon[END_REF] have suggested a substantial level of metal-silicate equilibration is needed to remove siderophile elements almost entirely from the rocky mantle to explain the small excess of 182 W of the Earth relative to the Moon. There are two possible mechanisms to enhance equilibration. The first one is to mechanically break the impactor's core into small pieces since the smaller metal pieces would have longer falling time and shorter chemical equilibration time. [START_REF] Dahl | Turbulent mixing of metal and silicate during planet accretion-and interpretation of the Hf-W chronometer[END_REF] show only iron fragments less than 10 km in diameter would equilibrate with silicates. [START_REF] Kendall | Differentiated planetesimal impacts into a terrestrial magma ocean: Fate of the iron core[END_REF] have futher revised this value by performing hydrodynamic simulations and found iron blobs with a radius of 100 km is in full equilibrium with the magma ocean. It should be noted that all Moon-forming impact simulations do not include the strength of materials, causing the overestimation of the size of fragmented iron core [START_REF] Barr | On the origin of earth's moon[END_REF]. The second mechanism is vaporisation. After cooling down, the vaporised materials would condensate into a distribution of small droplets on a centimetre level, and thus enhance the metal-silicate equilibration. [START_REF] Kraus | Impact vaporization of planetesimal cores in the late stages of planet formation[END_REF] developed an experimental technique to determine the shock pressure required for vaporization of iron along the principal Hugoniot line. They found the starting vaporization pressure is around 415 GPa compared to 817 GPa from the previous theoretical estimate using ANEOS (Pierazzo et al., 1997). Therefore, previous hydrodynamic simulations with ANEOS may underestimate the production of iron vapor.

Before the giant impact, the proto-Earth's core may grow with stable compositional stratification because higher abundances of light elements would be incorporated into the liquid metal as a result of the increasing metal-silicate equilibration pressure and temperature during accretion [START_REF] Jacobson | Formation, stratification, and mixing of the cores of earth and venus[END_REF]. The stable stratification would inhibit the outer-core convection and prevent from generating a geodynamo. [START_REF] Jacobson | Formation, stratification, and mixing of the cores of earth and venus[END_REF] proposed that an energetic giant impact may pump enough energy to homogenise the core as long as at least 4% of the total energy is deposited in the core. However, it remains as an open question on the exact amount of energy distributed into the Earth's core during the giant impact, which requires significant numerical simulations to clarify.

The role of equations of state

During the giant impact, the colliding interface and silicates in the disk can be heated up to 10 4 K and 7000-8000 K, respectively, even in the canonical impact model (Canup, 2004a). More surprisingly, the core of the impactor can reach as high as 45000 K (Canup, 2004a). It remains unclear whether these reported temperatures are physically reasonable and not caused by the equations of state of iron and silicates used in these simulations. For any giant impact simulation, the equations of state (EOS) are needed to describe the thermodynamic properties over a wide range of temperature, pressures and density. It will inevitably affect the energy distribution [START_REF] Nakajima | Melting and mixing states of the earth's mantle after the moon-forming impact[END_REF] and alter the after-impact dynamic evolution of the Earth such as the core-mantle equilibration and homogenising the compositional strafication of the Earth's core.

There are two sets of equations of state that have been widely used in the hydrodynamic simulations. The first one is called Tillotson EOS, which has a relatively simple analytic formula that makes computations very fast. The shortage is that Tillotson EOS does not provide any information on the physical state of materials. Therefore, simulation results are lack of vapor which may affect the efficiency of the pressure gradient to inject materials into the proto-lunar disk [START_REF] Stevenson | Origin of the moon-the collision hypothesis[END_REF].

The second one is ANEOS, which is more sophisticated and provides a thermodynamically consistent description of the phase diagram of iron and silicates. However, ANEOS treats the vapor phase as a monatomic mixture of atoms rather than the molecular cluster, resulting in an unrealistic large amount of entropy needed to vaporise solid or liquid into the gas. This problem has been solved in the improved version of ANEOS or M-ANEOS [START_REF] Melosh | A hydrocode equation of state for SiO 2[END_REF]. Canup (2004a) has compared the hydrodynamic simulation results using the Tillotson EOS with that of M-ANEOS, and found no difference in the origin of materials of the proto-lunar disk, where more than 60 wt% is from the impactor. Only more vapor is produced using M-ANEOS as expected, which is still not enough to make the pressure gradient as the major injection mechanism.

The phase diagram of iron

As a building-block of the terrestrial planet, iron naturally receives a lot of attention on its thermodynamic properties under various conditions. Many experimental techniques have been developed to determine its phase diagram and equation of state.

The laser-heated or resistance-heated diamond anvil cell is capable of compressing iron up to 400 GPa and 6000 K [START_REF] Sinmyo | Melting curve of iron to 290 gpa determined in a resistance-heated diamond-anvil cell[END_REF]. To reach pressures into the terapascal (TPa) range and temperatures up to 10,000 K or more, the dynamic compression is the only method available. It has been used to measure the physical properties of iron up to 1.4 TPa [START_REF] Smith | Equation of state of iron under core conditions of large rocky exoplanets[END_REF]. The rapid advance in the computer capability in the last two decades has made ab initio simulations feasible.

It is capable of accessing pressure and temperature conditions that are difficult to measure experimentally. In addition, theoretical simulations can directly compute the entropy (e.g. [START_REF] Alfe | The melting curve of iron at the pressures of the earth's core from ab initio calculations[END_REF], which is vital to build up an accurate equation of state and is challenging to obtain from experiments.

Figure 1.3 displays the phase diagram for iron in a pressure range of 0-360 GPa and a temperature range of 0 -8000 K. At ambient conditions, the most stable phase of iron adopts a body-centred cubic (α-bcc) structure and is in a ferromagnetic state. The Curie temperature, which marks a transition from the ferromagnetic state to the paramagnetic state, is about 1043 K. A further phase transition to the face-centred cubic (γ-fcc) structure occurs at 1200 K. In the temperature range of 1670 K to 1815 K, iron adopts the δ-bcc structure again. Above 1815 K, iron starts to melt and becomes liquid. At ambient temperature and with increasing pressure, iron changes into the hexagonal close-packed structure (ε-hcp) at around 16 GPa. This structure is stable up to at least 400 GPa and 6000 K, although recent studies suggest there is a phase transition from ε-hcp to bcc structure [START_REF] Belonoshko | Stabilization of body-centred cubic iron under inner-core conditions[END_REF] at the inner core conditions (360 GPa and 6000 K). The triple point between α-bcc, ε-hcp and γ-fcc phase of iron is around 10.5 GP and 753 K. Along the melting curve, there are two triple points. The first one is at 5 GPa and 2000 K coexisted with δ-bcc and γ-fcc phase of iron. The second one is at 90 GPa and 3000 K coexisted with ε-hcp and γ-fcc phase of iron.

A thermodynamically consistent EOS describing the Helmholtz free energy as a function of volume and temperature (F (V, T )) requires the entropy information. It is too expansive to perform a series of simulations over a wide range of pressure and temperature conditions to directly determine the entropy. The common practice is to build up a analytic EOS with several adjustable parameters that can be fitted [START_REF] Sinmyo | Melting curve of iron to 290 gpa determined in a resistance-heated diamond-anvil cell[END_REF]. The pressure at the core-mantle boundary (CMB) is around 130 GPa. The pressure at the inner-core boundary (ICB) is about 330 GPa. The empty red square is the melting temperature at ICB determined by the ab initio method [START_REF] Alfe | Temperature of the inner-core boundary of the earth: Melting of iron at high pressure from first-principles coexistence simulations[END_REF]. (b) The liquid-vapor equilibrium line in the temperature-density space. The black solid line and red dotted line represents the liquid-vapor equilibrium line estimated by [START_REF] Fortov | Shock waves and equations of state of matter[END_REF] and [START_REF] Grosse | The densities of liquid iron and nickel and an estimate of their critical temperature[END_REF], repectively. The blue point denotes the density of liquid iron at the boiling point of 1 bar. The magenta shared region is the experimental conditions done by [START_REF] Hixson | Sound speed and thermophysical properties of liquid iron and nickel[END_REF]. (c) The liquid-vapor equilibrium line in the temperature-pressure space. The supercritical state is represented by the shaded region. from the measured or calculated pressure-density-temperature relations and phase boundary of iron (e.g. [START_REF] Dorogokupets | Thermodynamics and equations of state of iron to 350 gpa and 6000 k[END_REF]. The derived EOS works very well to predict pressure as a function of temperature or density. However, the obtained entropy from F (V, T ) are sometimes questionable since it is related to the partial derivative of Helmholtz free energy to the temperature at a fixed volume and thus depends on the particular EOS model. For iron, the magnetic entropy in the paramagnetic phase further increases the complexity. [START_REF] Ruban | Temperature-induced longitudinal spin fluctuations in fe and ni[END_REF] employed a microscopic phenomenological model to study the effect of the longitudinal spin fluctuation and found hcp-Fe may acquire a significant magnetic moment at high pressure and temperature. Therefore, future studies are required to clarify whether we need to consider magnetic entropy for all phases of iron.

Above studies are mainly on the condensed phases (>7.8 g/cm 3 ) of iron. The research into the low-density region near the liquid-vapor equilibrium line is very scarce. [START_REF] Grosse | The densities of liquid iron and nickel and an estimate of their critical temperature[END_REF] estimated the critical point of iron based on the rectilinear law and the theorem of corresponding states. The rectilinear law states the density at the boiling point (5.828 g/cm 3 at around 3160 K) is 4.35 times more than that at the critical point. With this in mind, we can determine the critical density which is 1.34 g/cm 3 . Then the critical temperature (7000 K) can be extrapolated from the relation between the liquid density and temperature spanning from the melting point to the boiling point. The corresponding law states for a metal, the entropy change along the liquid-vapor equilibrium line is equal for all materials at the same reduced temperature which is defined as the ratio between the physical temperature and the critical temperature. [START_REF] Grosse | The densities of liquid iron and nickel and an estimate of their critical temperature[END_REF] have measured the entropy change at the boiling point of 3160 K for iron. This entropy change is equal to that for mercury at a reduced temperature of 0.31. Therefore, the critical temperature of iron is estimated to be about 10000 K (3160/0.31). Considering the significant difference of 3000 K given by these two methods, more work is needed to determine the liquid-vapor equilibrium and the critical point of iron. [START_REF] Hixson | Sound speed and thermophysical properties of liquid iron and nickel[END_REF] measured the density of liquid iron up to 4000 K in an argon gas atmosphere at a pressure of 0.2 GPa, where the high temperature was generated by the electrical-pulse-heating technique. With an argon-ion laser backlighting the sample, the density was calculated by measuring its sample diameter. [START_REF] Fortov | Shock waves and equations of state of matter[END_REF] has fitted a semi-empirical EOS by using available experimental data from [START_REF] Grosse | The densities of liquid iron and nickel and an estimate of their critical temperature[END_REF] and [START_REF] Hixson | Sound speed and thermophysical properties of liquid iron and nickel[END_REF], and predicted the critical point of iron at 8786 K and 1.638 g/cm 3 .

Proposed research

The lack of reliable information on the thermodynamic properties of iron in the low-density region prevents us from developing an accurate equation of state to describe its thermodynamic response to shock waves. Consequently, it would affect the impact energy distribution and thus cause different geochemical and geodynamic consequences such as the amount of iron vaporized and the extent of iron-silicates equilibration. The aim of the present work is to close the knowledge gap of the liquid-vapor equilibrium of iron and to improve our understanding of the behavior of iron during giant impacts. In particular, we 1. assess whether the planets' core undergoes partial vaporization during giant impacts (Chapter 3), 2. implement a new technique to compute the liquid-vapor equilibrium and the critical point for a metallic system and validate it with sodium for which experimental data are available (Chapter 4), 3. apply this technique to iron (Chapter 5).

Before that, a brief introduction of first-principles simulation is provided, which is fundamental for our research (Chapter 2). At last, a summary of the results is given and future work is outlined (Chapter 6).

Introduction

In the present study, we employ ab initio molecular dynamics (AIMD) and Monte Carlo method to determine the physical properties of iron in the low-density region since it is experimentally difficult to access the low-density and high-temperature conditions. This chapter's focus is on the molecular dynamics method, and the introduction of the Monte Carlo method will be left to Chapter 4. The roadmap of this chapter is as following:

1. We start with the time-dependent Schrodinger equation that governs electron and nucleus dynamics (Eq. 2.1). By applying the Born-Oppenheimer approximation, we show the fast electrons would adiabatically follow the slow nuclei, and stay at the ground state for a large band gap system or the thermal ground state for the small band gap system. If we treat nuclei as classic particles and neglect their quantum effects, the nuclei follow Newton's second law of motion (Eq. 2.10 and 2.9). The force exerted on the ions is obtained from the position derivative of the electronic potential energy. For a small band gap system, we show the electronic potential energy should include the contribution of the excited states (Eq. 2.11). For a large band gap system, the effect of the excited states is not pronounced, and thus the electronic ground state energy acts as the electronic potential energy.

2. The remaining problem is how to obtain the electronic potential energy, where the electron excited states must be included since we focus on the metallic iron at high temperature. The direct method is to find all eigenfunctions and eigenvalues for the electronic Hamiltonian. But for this, we need to know the kinetic energy and electron-electron coulomb energy operators for the interacting electron system. One the one hand, these terms are unknown. On the other hand, it is extremely difficult to find numerical solutions. Finite-temperature density functional theory provides an alternative and cost-effective way. It shows the electronic potential energy is a unique functional of the electronic density.

It simplify the problem since there are only three spatial coordinates in the electronic density. In order to find the thermal ground state electronic density, the Kohn-Sham scheme maps the interacting electron system onto a fictitious non-interacting system, where they share the same thermal ground state electronic density. The solution of Kohn-Sham equation leads to the electronic potential energy.

3. Although the iron atom in the solid phase below or close to the Curie temperature and at ambient pressure contains a non-vanishing magnetic moment, whether the liquid iron at a high temperature like 3500 K is still magnetic is uncertain. But if it does, we need to consider spin dynamics that occur in the paramagnetic state. However, the finite-temperature density functional theory and its spin-polarized extension fail to describe the magnetism at high temperature as they fail to capture the longitudinal and transversal fluctuation.

Since the aim is to build a conceptual framework, the mathematical vigour and proof are not pursued. Any reader being interested in mathematical physics may refer to the excellent book by [START_REF] Eschrig | The fundamentals of density functional theory[END_REF] for more details.

Deriving classical molecular dynamics

For a system with N nuclei and N e electrons, The starting point of our discussion is the time-dependent Schrodinger equation,

i ∂Ψ({R} N , {r} Ne , t) ∂t = HΨ({R} N , {r} Ne , t), (2.1)
where {r} Ne and {R} N denote a set of electron and nuclei coordinates, is the reduced Planck constant, t is the time and H is the total Hamiltonian operator the effective potential energy governing the ions' motion in small band gap system E0({R}N ) the potential energy governing the ions' motion in large band gap system Pi the momentum of the ith nucleus defined as,

H = T e + T i + V e-e + V e-i + V i-i , (2.2)
where T e and T i are the kinetic energy operator of electrons and ions, respectively, V e-e , V i-i and V e-i are the electron-electron, electron-nuclear and nuclear-nuclear Coulomb interaction operator, respectively. Here the electron coordinates do not include the spin degrees of freedom meaning spatial orbitals are the same in their spin-up and spin-down electrons so that the band occupancy can be up to 2.

As the computational cost increases as the square of the number of coordinates, we need to avoid treating a large amount of coordinates simultaneously. To solve this problem, we can employ the Born-Oppenheimer (BO) approach [START_REF] Born | Dynamical theory of crystal lattices[END_REF] to separate the fast electronic motion from the slow nuclear motion and to divide the total Hamiltonian into two parts,

H fast = T e + V e-e + V e-i + V i-i (2.
3)

H slow = T i . (2.4)
The solution of Eq. 2.3 at fixed nuclei configurations {R} N is given by,

H fast Φ l ({r} Ne ; {R} N ) = E l ({R} N )Φ l ({r} Ne ; {R} N ), (2.5)
where Φ l ({r} Ne ; {R} N ) and E l ({R} N ) are resulted eigenfunction and eigenvalue, here and thereafter the parametric dependence of any function on the nuclear coordinates is indicated by the symbol after the semicolon, and l denotes the electronic state. We stress that the parametric dependence means that

dΦ l ({r} Ne ;{R} N dR i is always zero but dE l ({R} N )
dR i may and may not be zero, where i denotes the ith nucleus.

If Eq. 2.5 is known for all nuclei configurations, we can expand the total wavefunction Ψ as,

Ψ({R} N , {r} Ne , t) = ∞ l=0 Φ l ({r} Ne ; {R} N )χ l ({R} N , t), (2.6)
where χ l ({R} N , t) is the expansion coefficient. The insertion of Eq. 2.6 into Eq. 2.1 followed by exploiting the orthonormalization condition of yields a coupled differential equations (see Eq. ( 2.6) in [START_REF] Marx | Ab initio molecular dynamics: basic theory and advanced methods[END_REF], which is challenging to handle. Therefore, some approximations need to be made to simplify the problem.

The large band gap system

For a system at a temperature that is much lower than its band gap, we can ignore the mixing between different electronic levels in Eq. 2.6 and use only the first term to expand the total wavefunction Ψ,

Ψ({R} N , {r} Ne , t) = Φ 0 ({r} Ne ; {R} N )χ 0 ({R} N , t). (2.7)
where Φ 0 ({r} Ne ; {R} N ) is the ground state wavefunction, defined in Eq. 2.5. It means electron always stays in its ground state (l = 0). Plugging the right-hand side term of both Eq. 2.7 and Eq. 2.5 into Eq. 2.1 leads to,

i ∂χ 0 ({R} N , t) ∂t = (H slow + E 0 ({R} N )χ 0 ({R} N , t) (2.8)
where E 0 ({R} N ) is ground state energy (l = 0) at the fixed nuclei configuration {R} N as defined in Eq. 2.5 when l = 0. If we treat the nuclei as classic particles and neglect the quantum effects, Eq. 2.8 can be transformed into,

M i dR 2 i dt 2 = -∇ i E 0 ({R} N ), (2.9)
where M i is the mass of the ith nucleus [START_REF] Marx | Ab initio molecular dynamics: basic theory and advanced methods[END_REF].

Since the de Broglie thermal wavelength of an iron atom is much less than the average atomic distance in fluid iron, it may be reasonable to negelect the nuclei quantum effects [START_REF] Hansen | Theory of Simple Liquids[END_REF]. Approximating Ψ({R} N , {r} Ne , t) with only one term in Eq. 2.7 and neglecting the quantum effects of nuclei inevitably cause errors, which are supposed to be small. We have to live with it in order to find a numerical solution.

2.2.2

The small band gap and metallic system

For a system that has a small band gap or is metallic like iron at high temperature (T ) , we can expand Ψ in a similar way as Eq. 2.7 but needs to take into account the thermal excitation of electrons explicitly. The electronic subsystem cannot be described by a pure state but by a mixed state (Eq. 2.14). We can employ the similar procedure as in Section 2.2.1, which leads to,

M i dR 2 i dt 2 = -∇ i F 0 ({R} N ), (2.10)
F 0 ({R} N ) = - 1 β ln Tr exp (-βH fast ) (2.11)
where β = 1/k B T , k B is the Boltzmann's constant and F 0 ({R} N ) is the effective potential energy of the electron subsystem at fixed nuclei configurations {R} N including the free energy of the electrons and the coulombic attraction within nuclei.

The details on the derivation of Eq. 2.11 can be found in [START_REF] Zwanzig | Transition from quantum to" classical" partition function[END_REF].

Eq. 2.11 suggests that electrons response instantaneously to the ionic motions and reach the thermal equilibrium [START_REF] Alavi | Ab initio molecular dynamics with excited electrons[END_REF]. We can also re-write Eq. 2.10 in a more compact way [START_REF] Zwanzig | Transition from quantum to" classical" partition function[END_REF] that will be used in Chapter 4,

Z = 1 N !h 3N exp -β N i=1 P 2 i 2M i + F 0 ({R} N ) N i dR i dP i . (2.12)
where Z is the canonical partition function, N is the total number of ions, and P i is the momentum of the ith nucleus. In Eq. 2.12, the electronic contribution to the partition function has been treated in a quantum-mechanical way, while the nuclei has been treated in a classic way under the action of the free energy of the electron system (Eq. 2.11).

Summary

In summary, this section has presented a practical way to solve Eq. 2.1 by separating the electronic motions from the nuclear motions. If we treat nuclei as classic particles and neglect their quantum effects, the nuclei follow Newton's second law of motion (Eq. 2.10 and 2.9). The force is obtained from the position derivative of the electronic potential energy. For the metallic system like iron, the electronic potential energy should include the contribution of the excited states (Eq. 2.11). In order to maintain the temperature of the system at a desired value (T ), the Nosé thermostat [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF] is used. The conserved quantity becomes [START_REF] Wentzcovitch | Energy versus free-energy conservation in first-principles molecular dynamics[END_REF], .13) where K i the kinetic energy of ions, U thermo and K thermo are the potential and kinetic energy of the thermostat, respectively. It should be noted that [START_REF] Wentzcovitch | Energy versus free-energy conservation in first-principles molecular dynamics[END_REF] suggest F 0 ({R} N ) can be calculated at an arbitrary temperature. However, we have shown in Eq. 2.11 that the electron subsystem must stay at the simulated temperature (T ). For the sake of simplicity, we omit hereafter the parametric dependence of Φ l ({r} Ne ; {R} N ) on the ionic position {R} N . We also neglect the dependence of any function on independent variables like {R} N and {r} Ne . For instance, F 0 ({R} N ) will be written as F 0 . the probability of the system being found in the state

E total = F 0 ({R} N ) + K i + U thermo + K thermo . ( 2

Finite-temperature density functional theory

ψ l F [Γ]
the Helmholtz free energy of the electronic subsystem as described by Γ Γ0

the density matrix to give the minimal F = F0 for the electron subsystem {Φ 0l } a complete orthonormal basis set to form Γ0, and equal to Φ l ({r}N e ; {R}N ) p 0l the probability of the system being found in the state Φ 0l n the electron density S the entropy of an interacting electron system νext the external potential due to electron-ion interactions Ve-i Ω the grand potential energy in the grand canonical ensemble n0 the ground state electron density that gives F0 Te the kinetic energy of a non-interacting electron system Ss the entropy of a non-interacting electron system J [n] the classic Coulomb energy for a non-interacting electron system Fxc [n] the exchange-correlation contribution to the free energy ν ks

the Kohn-Sham external potential φ l an arbitrarily complete orthonormal basis set in the single particle Hilbert space f l the probability of the system being found in the state φ l

Finite-temperature canonical-ensemble theory

In Eq. 2.11, the evaluation of F 0 requires information of all eigenvalues and eigenstates of H fast , which needs a tremendous computational cost. The minimum principle may provide an alternative method that is easier to implement. Before describing this theorem, we need to introduce some terminology. For the canonical ensemble considering here, the electron subsystem is described by a density matrix, which is characterized by a probability distribution over all pure states that have the same particle number N e ,

Γ = ∞ l=1 p l |ψ l ψ l | , p l ≥ 0, ∞ l p l = 1 (2.14)
where the Dirac bra-ket notation is used, {ψ l } forms a arbitrarily complete orthonor-mal basis set in the N e particle Hilbert space and is not necessarily equal to the eigenstates of H fast , and p l is the probability of the system being found in the state |ψ l . Then we can define the Helmholtz free energy F as,

F [Γ] = F [{p l , ψ l }] = ∞ l=1 p l ( 1 β lnp l + < ψ l |T e + V e-e + V e-i |ψ l >) + V i-i . (2.15)
The quantum statistical mechanics has shown in the thermal equilibrium, the Helmholtz free energy is at its minimum (see Chapter 3 in [START_REF] Parr | Density functional theory of atoms and molecules[END_REF],

F [Γ] ≥ F [Γ 0 ] = F [{p 0l , Φ 0l }] = F 0 (2.16)
where Φ 0l are the eigenstates of H fast at the fixed configuration {R} N , and p 0l is defined as,

p 0l = exp(-βE l ) ∞ i exp(-βE l ) , (2.17)
where E i is the eigenvalues of H fast . Eq. 2.16 sets up the variational principle to find F 0 . However, the computational complexity is very similar to the direct calculation of eigenvalues and eigenstates of H fast .

A more cost-effective method is density functional theory (DFT), which allows to perform the minimization over the electron density which depends on only three spatial coordinates. For a system with non-degenerate ground state, the Hohenberg-Kohn (HK) theorem [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF] proves there is one-to-one correspondence between the external potential, the resulted non-degenerate ground state, and the associated ground state density. Therefore, any ground state observable including the ground state energy is a unique functional of electron density. The extension of this theorem to high temperature was pioneered by [START_REF] Mermin | Thermal properties of the inhomogeneous electron gas[END_REF], just one year after Hohenberg and Kohn published their seminal work. For the discussion of finite-temperature density functional theory, we need to define ensemble N -representable electron density associated with a density matrix (Γ) defined in Eq. 2.14 as,

n(r) = N e ∞ l=0 p l d 3 r 2 d 3 r 3 • • • d 3 r N ψ * l (r, r 2 , . . . , r N )ψ l (r, r 2 , . . . , r N ).
(2.18) and rewrite the electron-nuclear Coulomb interaction operator in Eq. 2.5 as a functional of the electron density,

V e-i = d 3 rν ext (r)n(r).
( 2.19) where ν ext is the external potential due to electron-ion interactions. [START_REF] Mermin | Thermal properties of the inhomogeneous electron gas[END_REF] proved that in the grand canonical ensemble at a given temperature and chemical potential, no two external potentials can lead to the same equilibrium density. This suggests the grand potential is a unique functional of electron density. Based on this fact, we can define a universal functional (F L [n]) of the electron density which is independent of the external potential and the quantity

Ω[n] = drν ext n(r) + F L [n] (2.20)
is at its minimum and equal to the grand potential Ω 0 when n(r) is the equilibrium density n 0 (r). [START_REF] Mermin | Thermal properties of the inhomogeneous electron gas[END_REF] also mentioned about these arguments can be adapted to the canonical ensemble with a few minor changes, which has been given in the Chapter 3 of [START_REF] Parr | Density functional theory of atoms and molecules[END_REF]. With the help of the Lieb functional, we can contruct the density functional theory for the minimum principle in Eq. 2.15 as,

F 0 = min n F [n] = min n (F L [n] + V e-i [n]) + V i-i (2.21)
where the minimization is performed over all ensemble N -representable electron densities defined in Eq. 2.18, and F L [n] is a universal density functional independent of the external potential ν ext ,

F L [n] = T e [n] + V e-e [n] + S[n] = min {ψ l ,p l }→n ∞ i p l ( ψ l | T e |ψ l + ψ l | V e-e |ψ l + 1 β lnp l ) (2.22)
where {ψ l , p l } → n is a constrained search over all sets {ψ l , p l } that can generate n and S[n] is the entropy of the interacting electron system. For more details on the Lieb functional and N -representable electron density, please refer to Chapter 2 in [START_REF] Engel | Density functional theory: An Advanced Course[END_REF]. If the definition domain of F [n] includes densities that integrate to a fractional number of particles and δF δn exits, the minimum principle in Eq. 2.21 indicates it is possible to determine the ground state density by a variational equation,

δ δn(r) F [n] -µ( d 3 rn(r) -N e ) = 0 (2.23)
where the Lagrange multiplier µ is used to impose the constraint,

d 3 rn(r) = N e .
(2.24)

Kohn-Sham formulation

We have built a varitional principle to find the F 0 based on Eq. 2.23. However, we still have no information on how to calculate ψ l | T e |ψ l and ψ l | V e-e |ψ l . The Kohn-Sham formulation [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] employs a mapping from the full interacting system with V e-e = 0 onto a fictitious non-interacting (V e-e = 0) system where the electrons move within an effective potential ν ks ,

ν ext V e-e =0 ← ---→ n 0 V e-e =0 ← ---→ ν ks , (2.25)
where we have assumed interacting ensemble N -representable electron density is also non-interacting ensemble N -representable. The Kohn-Sham method is exact since the non-interacting system yields the same ground-state density as the real system. We can then transformed Eq. 2.21 into,

F 0 = min n(r) (T e [n] + V e-e [n] + S[n] + V i-i )(interacting electron system) = min n(r) (T s [n] + J[n] + F xc [n] + S s [n] + V i-i )(non-interacting electron system) (2.26)
T s is the kinetic energy of non-interacting electron system, S s is the entropy of non-interacting electron system, J[n] is the classic Coulomb energy by,

J[n] = dr dr n(r )n(r) |r -r| , (2.27)
and F xc [n] is the exchange-correlation contribution to the free energy that includes the difference of the kinetic energy, entropy and Coulomb energy between the interacting electron system and the non-interacting electron system,

F xc [n] = (T [n] + S[n]) -(T s [n] + S s [n]) + (V e-e [n] -J[n]).
(2.28)

For the non-interesting system, we can express the electron density n with respect to the single particle states as,

n(r) = ∞ l=0 f l φ * l (r)φ l (r), f l ≥ 0, ∞ l=0 f l = 1 (2.29)
where {φ l } forms a arbitrarily complete orthonormal basis set in the single particle Hilbert space and and f l is the probability of the system being found in the state |φ l .

Plugging Eq. 2.29 into Eq. 2.26 leads to,

F 0 = min {φ l ,f l } ∞ l=0 [f l φ l | - 1 2 ∇ 2 |φ l + 1 β (f l lnf l + (1 -f l )ln(1 -f l ))] + dr dr n(r )n(r) |r -r| + F xc [n] (2.30)
Based on the conclusion in Eq. 2.16, the minimization is obtained if φ l are the eigenfunctions

[- 1 2 ∇ 2 + ν ks (r)]φ l (r) = ε l φ l (r) (2.31)
and f l is

f l = 1 1 + exp[β(ε l -µ]) (2.32)
where ε l is the eigenvalues and µ is the chemical potential calculated from the constraint,

∞ l=0 f l = N e (2.33)
and ν ks is defined as,

ν ks (r) = ν ext (r) + n(r ) |r -r | + δF xc [n] δn .
(2.34) Eq. 2.31 have to be solved self-consistently. And the resulted equilibrium density can be inserted into Eq. 2.26 to obtain F 0 (see Chapter 9 in [START_REF] Parr | Density functional theory of atoms and molecules[END_REF]).

In the finite-temperature density functional theory, the physically meaningful quantity is F 0 and its analytic derivative with respect to the atomic position and cell volume from which the force and electronic pressure are calculated. It is questionable to use the quantity F 0 -S s as the internal energy for the interacting electron system since S s is the electronic entropy of the fictitious non-interacting electron system. Besides, we cannot determine the heat capacity for the metallic system based on the fluctuation-dissipation theorem since FT-DFT provides no information on the fluctuation of the electronic internal energy.

Ab initio spin dynamics

Spin-polarized DFT at zero temperature

Magnetism comprises a large variety of phenomena that can be characterized by quantities such as the magnetic moment, the magnetic order and the ordering temperature and so on. Although there is no single theoretical approach to all magnetic phenomena, the spin-polarized version of DFT allows the access to the magnetic moment and the magnetic order at zero temperature. It can also serve to extract other quantities like exchange interactions that can serve as input for other theoretical approaches. In this section, we will present how to extend the Kohn-Sham formulation described in Section 2.3.2 to the spin-polarized case at zero temperature.

The basic variables of spin density functional theory are the scalar electronic density n(r) and the vector of the magnetization density m(r). von [START_REF] Von Barth | A local exchange-correlation potential for the spin polarized case. i[END_REF] have extended the DFT concept to spin-polarized systems by defining the Hohenberg-Kohn-Sham spin density functional as,

E 0 [n(r), m(r)] = min n(r),m(r) [T e [n(r), m(r)] + V e-e [n(r), m(r)] + V i-i ] (2.35)
It means the internal energy is at a minimum when n(r) and m(r) are the equilibrium electron and magnetization density, repectively. Here we choose the internal energy E 0 instead of the free energy F 0 since the electronic entropy at zero temperature is zero and E 0 = lim T →0 F 0 . Then we can map the interacting electron system to a non-interacting electron system and express n(r) and m(r) with respect to the single particle states. The details are not given here, which can be found in [START_REF] Bihlmayer | Density-functional theory of magnetism[END_REF]. In this way, we are able to derive a spin-polarized version of Kohn-Sham equations without a magnetic field,

- 1 2 ∇ 2 + n(r ) |r -r | + ν ext (r) 1 0 0 1 + δE xc δn(r) × φ ↑ l φ ↓ l = ε l × φ ↑ l φ ↓ l , (2.36)
where ↑ and ↓ denote the spin-up and spin-down states, and the density matrix is defined as,

n α,β (r) = ∞ l=1 φ * ,α l (r)φ β (r), α, β =↑ or ↓ . (2.37)
It can be decomposed into a scalar and vectorial part, correponding to the charge and magnetization density,

n(r) = 1 2 n(r) 1 0 0 1 + σm(r) = 1 2 n(r) + m z (r) m x (r) -im y (r) m x (r) + im y (r) n(r) -m z (r) ,
(2.38) where σ is the Pauli matrices. And the magnetic moment is defined as,

M I = Ω I m(r)dr, e I = M I |M I | , (2.39)
where Ω I is a sphere centered at Ith atom and e I is the orientation of the local magnetic moment.

Spin dynamics at high temperature

We can use the method described in the last section to simulate magnetic systems and determine their ground state accurately. But at high temperature, we need to consider the magnetic excitations consisting of the Stoner single particle and collective excitation [START_REF] Kaul | Phase transitions and finite temperature magnetism: Experiment and analysis[END_REF]. The single particle excitation is due to the transfer of an electrons-hole pair in the band of opposite spins. The collective excitation is caused by a spin flip from atom to atom and can be further divided into the transversal and longitudinal fluctuation [START_REF] Kaul | Phase transitions and finite temperature magnetism: Experiment and analysis[END_REF]. They would change the orientation and magnitude of the local magnetic moment from one site to the other, respec-tively. When the temperature reaches the Curie point or above, the transversal and longitudinal fluctuation destroy the long-range magnetic order and the ferromagnetic system would become paramagnetic. However, the individual atom may still contain a non-vanishing magnetic moment. We can explain the presence of spin dynamics or magnetic excitations at high temperatures by analogy with atom vibrations in a perfect lattice. For the N V T ensemble in thermodynamic equilibrium, the Helmholtz free energy should be at its minimum. At high temperatures, atom vibration increases the internal energy but brings in extra entropy due to the increased disorder. The overall effect is the entropy gain is larger than the cost of the increasing internal energy. Similarly, although spin dynamics increases the internal energy, the system would obtain more entropy due to the magnetic disorder.

Electronic structure calculations within DFT can be viewed as a modern extension of the Stoner-type description of magnetism [START_REF] Abrikosov | Recent progress in simulations of the paramagnetic state of magnetic materials[END_REF]. It is capable of reproducing ground-state magnetic properties of 3d transition metals accurately, whose magnetism comes from the imbalance between the number of spin-up and spin-down itinerant 3d electrons. However, DFT fails to describe the magnetism at finite temperature and significantly overestimates the Curie temperatures of iron by a factor of five due to its inability of treating the transversal and longitudinal fluctuation [START_REF] Abrikosov | Recent progress in simulations of the paramagnetic state of magnetic materials[END_REF].

To find a scheme to capture these fluctuations at high temperatures, we start again with the time-dependent Schrodinger equation but need to include the spin degrees of freedom explicitly, then precede to separate different dynamic time scales. Following this, [START_REF] Antropov | Ab initio spin dynamics in magnets[END_REF] partition the total Hamiltonian into a slow-motion part that relates to the spin density and nuclear motion. By introducing the rigid spin approximation and neglecting the longitudinal fluctuation, the time evolution of the spin density is described by a rotation of the magnetization density inside the atomic spheres. The remaining electronic Hamiltonian at fixed atomic positions and the local magnetic moment direction plays the role of the potential energy for the movement of nuclei and the magnetic moment direction. Therefore, we can write the equation of motion for the direction of the local magnetic moment and nuclei [START_REF] Ma | Dynamic simulation of structural phase transitions in magnetic iron[END_REF] as,

Z = 1 N !h 3N exp -β N I=1 P 2 I 2M I + F 0 ({R} N , {e} N ) N i dR i dP i de i .
(2.40) where e i is a unit vector of the local magnetic moment associated with ith atom. F 0 ({R} N , {e} N ) is the effective potential energy at fixed {R} N and {e} N that can be tacked by finite-temperature constrained spin-polarized density functional theory [START_REF] Ma | Constrained density functional for noncollinear magnetism[END_REF]. It means we need to extend the theory in Sec. 2.4.1 to include the electronic entropy and replace E 0 [n(r), m(r)] with F 0 [n(r), m(r)]. In addition, the band occupancy needs to be chosen based on Eq. 2.32. Band structure calculations suggest the magnetization density in Fe is well localized and each iron atom could be associated with a local magnetic moment that behaves in a Heisenberg-like manner [START_REF] Abrikosov | Recent progress in simulations of the paramagnetic state of magnetic materials[END_REF]. Therefore, the rigid spin approximation might be reasonable at least for iron. However, we cannot apply this method to Ni since it does not develop a well localized magnetic moment and the longitudinal fluctuation plays a dominant role [START_REF] Abrikosov | Recent progress in simulations of the paramagnetic state of magnetic materials[END_REF].

Both experimental [START_REF] Waseda | Atomic distribution and magnetic moment in liquid iron by neutron diffraction[END_REF] and theoretical studies [START_REF] Lichtenstein | Finite-temperature magnetism of transition metals: An ab initio dynamical mean-field theory[END_REF] have suggested the liquid Fe is in the paramagnetic state with a nonvanishing local magnetic moment. Although we may use the method described above to study the effect of spin dynamics, it is too expensive to do so. We need to further simplify the problem into a manageable size. For this, we re-consider the argument of dynamic timescales. [START_REF] Abrikosov | Recent progress in simulations of the paramagnetic state of magnetic materials[END_REF] suggest the spin de-coherence time is 10 -14 -10 -15 s which is very close to the characteristic timescale of electrons (10 -15 s). Therefore, we may apply the adiabatic approximation to the orientation of the local magnetic moment as well and re-write Eq. 2.40 as,

Z = 1 N !h 3N exp -β N I=1 P 2 I 2M I + F 0 ({R} N ) N i dR i dP i . (2.41)
where we have integrated out the magnetic degrees of freedom compared to Eq. 2.40, and F 0 ({R} N ) is defined as,

F 0 ({R} N ) = - 1 β {e}m exp(-βF l ({e} m , {R} N )) (2.42)
And the force on the {R} N is,

f = m p m f m ({e} m , {R} N ), p m = exp[-F 0 ({e} m , {R} N )] ∞ i exp[-F 0 ({e} i , {R} N )] , (2.43)
which can be calculated by the self-consistent KKR-CPA method [START_REF] Pindor | Disordered local moment state of magnetic transition metals: a self-consistent KKR CPA calculation[END_REF]. This equation is equivalent to the disordered local momentum picture introduced by [START_REF] Gyorffy | A first-principles theory of ferromagnetic phase transitions in metals[END_REF]. However, the KKR-CPA method has difficulties to be integrated into the plane wave density functional theory code. Instead, [START_REF] Alling | Effect of magnetic disorder and strong electron correlations on the thermodynamics of crn[END_REF] develop the magnetic sampling method to obtain the force by averaging over a few magnetic configurations with randomly distributed spin-up and spin-down local magnetic moments,

f = Nc i=1 f i ({e} i , {R} N ) N c . (2.44)
Although Eq. 2.44 differ slightly from Eq. 2.43, [START_REF] Alling | Effect of magnetic disorder and strong electron correlations on the thermodynamics of crn[END_REF] have shown for a large amount of magnetic configurations these two formula are equivalent. In present work, we will employ the magnetic sampling method to study the effect of the paramagnetic state on the physical properties of liquid iron at high temperatures and also discuss the magnetic state on its critical point in Chapter 5.

Approximate paramagnetism by either non-magnetism or ferromagnetism?

One might wonder whether it is appropriate to approximate paramagnetism by non-magnetism or ferromagnetism as the application of Eq. 2.41 to study the paramagnetic phase is expensive. We can solve this problem from the perspective of the Helmholtz free energy of the electron-ion system (not just the electronic free energy!). If the magnetic fluctuation in the paramagnetic phase is substantial, an extra magnetic entropy would be induced. Then neither the non-magnetic state nor the ferromagnetic phase shares the same Helmholtz free energy as the paramagnetic phase.

Despite this, it is conventional and quite popular to model the paramagnetic state by the ferromagnetic state for iron-related materials. Then an empirical correction for the magnetic entropy is added to the free energy of the ferromagnetic state. The widely used term for the magnetic entropy is k B ln(m + 1), where m is the average magnitude of the local magnetic moment. It is derived from the quantum mechanical entropy per magnetic atom (k B ln(2S + 1)), where for a given total spin angular momentum S its projection on the z-axis can take 2S + 1 values. As discussed by [START_REF] Khmelevskyi | Longitudinal integration measure in classical spin space and its application to first-principle based simulations of ferromagnetic metals[END_REF], the generalization from the quantum spin space to the classic spin space is problematic. For a local magnetic moment with the magnitude of m, it can point to more than m + 1 directions. [START_REF] Tsuchiya | Spin transition in magnesiowüstite in earth's lower mantle[END_REF] used the term k B ln(n(2S + 1)) to calculate the paramagnetic entropy of magnesiowustite, where n is the electronic configuration degeneracy whose value is 3 for the high-spin phase and 1 for low-spin phase. It is unclear why the electronic configuration degeneracy enters the magnetic entropy term. To show this, in the high-spin magnesiowustite, there are six electrons for Fe 2+ ; two of them occupy the two e g orbitals and three of them the three t 2g orbitals. The remaining electron would occupy one of the three t 2g orbitals. We denote the occupancy as p 1 , p 2 , p 3 , where the sum of p 1 + p 2 + p 3 is one. In each case, the quantum magnetic entropy is still k B ln(2S + 1). Therefore, the total quantum magnetic entropy is

S total = p 1 × k B ln(2S + 1) + p 2 × k B ln(2S + 1) + p 3 × k B ln(2S + 1) = (p 1 + p 2 + p 3 ) × k B ln(2S + 1) = k B ln(2S + 1).
(2.45)

Although one electron can occupy one of the three t 2g orbitals in the high-spin state of magnesiowustite, the total quantum magnetic entropy is k B ln(2S + 1). Considering the great uncertainty in the magnetic entropy, we should avoid estimating the free energy of the paramagnetic phase by giving an empirical correction to the free energy of the ferromagnetic phase.

Introduction

Once differentiated, the core of planets and planetesimals is dominated by liquid or solid iron, alloyed with nickel and various lighter elements [START_REF] Hirose | Composition and State of the Core[END_REF].

Because of its obvious geophysical significance, considerable effort was put to determine both theoretically and experimentally the phase diagram of iron [START_REF] Alfe | The melting curve of iron at the pressures of the earth's core from ab initio calculations[END_REF]Campbell, 2016;Caracas, 2016;[START_REF] Tateno | The Structure of Iron in Earth's Inner Core[END_REF] up to Earth's inner core conditions (around 360 GPa and 6000 K) and beyond. Recently, a complete set of equations of state (EOS) was proposed, covering 7 -30 g/cm 3 densities and 10,000 -1,000,000 K temperatures [START_REF] Sjostrom | Quantum molecular dynamics of warm dense iron and a five-phase equation of state[END_REF].

Much less is known about the thermodynamic and thermophysical properties of iron below 1 GPa and above its melting point (around 1850 K). [START_REF] Hixson | Sound speed and thermophysical properties of liquid iron and nickel[END_REF] have measured the liquid iron density at 0.2 GPa up to 4000 K. [START_REF] Grosse | The densities of liquid iron and nickel and an estimate of their critical temperature[END_REF] measured the liquid-vapor equilibrium density up to the boiling point at 1 bar and 3160 K. In order to obtain the critical point, a large extrapolation must be made. The first method is to employ empirical equations of state with several adjustable parameters, which can be determined from available experimental data [START_REF] Fortov | Shock waves and equations of state of matter[END_REF][START_REF] Medvedev | Wide-range multiphase equation of state for iron[END_REF]; the second one is to use the law of rectilinear diameter [START_REF] Grosse | The densities of liquid iron and nickel and an estimate of their critical temperature[END_REF]. However, it is unclear whether these extrapolations work at high temperature where no experimental data are available.

Indeed, the regime of low densities and high temperatures, still not yet well characterized, is typical for the after-shock state of proto-planetary cores occurring in the aftermath of catastrophic events such as giant impacts. The Earth's Moon formed after such a giant impact between the proto-Earth and Theia, an astronomical body whose most commonly accepted size is that of Mars [START_REF] Asphaug | Impact Origin of the Moon?[END_REF]Canup, 2004b). Hydrodynamic impact simulations show that it results in the formation of a disk [START_REF] Canup | Forming a Moon with an Earth-like Composition via a Giant Impact[END_REF][START_REF] Canup | The University of Arizona space science series[END_REF][START_REF] Ćuk | Making the moon from a fast-spinning earth: a giant impact followed by resonant despinning[END_REF] or a synestia [START_REF] Lock | The origin of the moon within a terrestrial synestia[END_REF]. The disk might be iron-depleted, producing a small Moon's core. However, results of these simulations heavily rely on available EOS. An experimental result on iron found the shock pressure required for vaporization when compressed from ambient conditions and then decompressed to 1 bar to be around 507 (+65, -85) GPa [START_REF] Kraus | Impact vaporization of planetesimal cores in the late stages of planet formation[END_REF], lower than previous estimates of 887 GPa [START_REF] Pierazzo | A reevaluation of impact melt production[END_REF]. This implies that the cores of a large number of the planetesimals from the late stage of accretion largely vaporized during the impacts [START_REF] Kraus | Impact vaporization of planetesimal cores in the late stages of planet formation[END_REF].

In order to assess whether the core of the planets undergoes significant vaporization during a giant impact, we employ ab initio molecular-dynamics simulations to explore iron over a wide density region encompassing the critical point (CP) and the Hugoniot lines of the shocked cores. In addition, we characterize the structural and transport properties of iron including diffusion coefficients, viscosity and thermal conductivity in the low-density region, which may provide a better constraint on the behavior of iron in the proto-lunar disk.

Simulation details

Ab initio molecular dynamics

Ab initio molecular dynamics simulations were performed by VASP code (Kresse and Furthmüller, 1996;Kresse and Furthmüller, 1996). We used the N V T ensemble, where N = 108 atoms, V is the volume of simulated system, and the temperature T is maintained by the Nosé thermostat [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF]. A Verlet algorithm is used to integrate the Newton's equation of motion with a timestep of 1 fs. The total simulation time at each temperature and density condition is at least 10 ps. The interatomic forces were calculated by employing the projector augmented wave (PAW) method [START_REF] Blöchl | Projector augmented-wave method[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] within the framework of finite-temperature density functional theory [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF][START_REF] Mermin | Thermal properties of the inhomogeneous electron gas[END_REF]. The PBE formalism [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] for the exchange correlation term and the valence electron configurations (3d 7 4s 1 ) for the pseudopotential were used. The partial occupancy for the electronic calculation were calculated using a Fermi smearing scheme with a width at simulated temperature. The energy cut-off for the plane-wave basis set was set to 550 eV. The break condition for the electronic self-consistent loop was 10 -4 eV. The number of electronic bands was adapted to the temperature conditions such as to cover the entire spectrum of the fully and partially occupied states and to include enough non-occupied bands. The convergence of the pressure tensor and the energy are on the order of a few percent when compared to a grid of 4×4×4 k-points generated by Monkhorst-Pack scheme [START_REF] Monkhorst | Special points for brillouin-zone integrations[END_REF]) and a kinetic energy cut-off of 850 eV.

It should be noted that both experiments [START_REF] Waseda | Atomic distribution and magnetic moment in liquid iron by neutron diffraction[END_REF] and theoretical simulations [START_REF] Lichtenstein | Finite-temperature magnetism of transition metals: An ab initio dynamical mean-field theory[END_REF] suggest that liquid iron might be in a paramagnetic state. As discussed by [START_REF] Marqués | Ab initio study of the structure and dynamics of bulk liquid Fe[END_REF], spin-polarized MD simulations yield an inherent long-range ferromagnetic order. In order to avoid such residual magnetic state, we decided to perform non-spin-polarized simulations to approximate the paramagnetic state of liquid iron at low pressure and high temperature. This is the mean field approximation of the paramagnetic state, even if it neglects the spin fluctuations whose effects will be discussed in Chapter 5.

Construction of the spinodal line

During the simulations at low temperatures, with decreasing density, the pressure reaches a local minimum. This marks the liquid spinodal point: the minimal density at which the liquid is stable. At densities lower than the spinodal, the liquid is unstable and cavitation occurs [START_REF] Speedy | Stability-limit conjecture. an interpretation of the properties of water[END_REF]. Under further expansion, the pressure starts to increase and the local maximum marks the gas spinodal: the maximum density at which the gas is metastable. Between the gas and liquid spinodal den-sities, neither gas nor liquid can exist as a single phase, but rather they co-exist. This is similar to the van der Waals model. In order to fit the pressure-density curves we employ a simple third-order polynomial function. This method has been successfully used by other theoretical studies on super-cooled silicon [START_REF] Vasisht | Liquid-liquid critical point in supercooled silicon[END_REF]. Spinodal lines with negative pressure have been reported in experiments [START_REF] Green | Water and Solutions at Negative Pressure: Raman Spectroscopic Study to -80 Megapascals[END_REF], classic MD simulations on the metastable extension of liquid water [START_REF] Poole | Phase behaviour of metastable water[END_REF], and first-principles MD on the metastable extension of liquid silicon [START_REF] Zhao | Phase behavior of metastable liquid silicon at negative pressure: Ab initio molecular dynamics[END_REF].

Structural analysis

The short-range order in the fluid is revealed by the pair distribution function (g F e-F e (r)),

g F e-F e (r) = V N (N -1) i j =i δ(r -r ij ) (3.1)
where r is the correlation length, r ij is the distance between atoms i and j, N is the total number of iron atoms in the simulations, V is the total volume of the system, and is the ensemble average.

The mean-square displacements

The mean-square displacements (MSD) are defined as,

M SD(τ ) = [r i (t 0 + τ ) -r i (t 0 )] 2 (3.2)
where the values are averaged over the total number of atoms and time origins (t 0 ) and τ represents a sliding time window spanning a portion of the trajectory. The asymptotic slope of the mean square displacement with respect to time yields the diffusion coefficient D in the long-time limit,

D = 1 6 lim τ →∞ M SD(τ ) τ .
(3.3)

Velocity autocorrelation function

We calculate the normalized velocity autocorrelation function as,

Φ(τ ) = v(0)v(τ ) v(0)v(0) , (3.4)
where v(t) is the velocity of a particle at time t. The self-diffusion coefficient D is given by the time integral of the velocity autocorrelation function by,

D = lim t→∞ 1 3 t 0 dτ v(0)v(τ ) (3.5)

Entropy calculations

The Fourier-transform of the velocity auto-correlation function yields the total movement of the atoms in the fluid defined as,

F (ν) = ∞ 0 Φ(τ ) cos 2πντ dτ (3.6)
The entropy can then be obtained by integrating over the vibrational part of this spectrum, in the same way as we do for solids. However, Eq. 3.6 captures not only the agitation of the atoms but also their diffusion. The latter is zero in solids, which allows us to directly obtain the entropy; but for fluids by definition it is finite and positive, and thus must be removed from the spectrum of Eq. 3.6. For this we employ the two-phase thermodynamic method [START_REF] Lin | The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: Validation for the phase diagram of Lennard-Jones fluids[END_REF] to decompose the total spectrum of Eq. 3.6 into a diffusive, gas-like part and a purely vibrational solid-like part,

F (ν) = (1 -f g )F s (ν) + f g F s (ν) (3.7)
where f g is the gas-like fraction. The entropy stemming from the gas-like and the solid-like parts is obtained using the hard sphere model and the harmonic oscillator model respectively. This method gives a reasonable estimation of entropy for pure liquid metals [START_REF] Desjarlais | First-principles calculation of entropy for liquid metals[END_REF]. We verify again our implementation and use the same parameters as [START_REF] Desjarlais | First-principles calculation of entropy for liquid metals[END_REF] and conduct simulations at 0 GPa and 1800 K for liquid iron. We obtain a value for the entropy of 11.05 k B /atom, compared to 12.00 k B /atom in [START_REF] Desjarlais | First-principles calculation of entropy for liquid metals[END_REF] . The discrepancy comes from the magnetic entropy, which is estimated to be about 1 k B /atom [START_REF] Desjarlais | First-principles calculation of entropy for liquid metals[END_REF]) and which we did not include in our calculation.

Viscosity

The viscosity is given by the famous Green-Kubo formula,

η = V k B T ∞ 0 P αβ (t)P αβ (0) dt (3.8)
where V is the total volume, k B is the Boltzmann constant, T is the temperature, P αβ stands for the off-diagonal element of the stress tensor and t denotes the correlation time.

Electrical and thermal conductivity

The electrical and thermal conductivity were calculated based on 20 configurations, equally spaced by 500 fs using the Kubo-Greenwood formalism [START_REF] Greenwood | The Boltzmann Equation in the Theory of Electrical Conduction in Metals[END_REF][START_REF] Kubo | Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems[END_REF]. This was successfully applied to the Earth's outer core [START_REF] Pozzo | Thermal and electrical conductivity of iron at Earth's core conditions[END_REF]. The Brillouin zone was sampled with the Baldereschi point [START_REF] Baldereschi | Mean-value point in the brillouin zone[END_REF]. A test with a k-points mesh of 4 × 4 × 4yields results within 1 % difference to using only the Baldereschi point.

We calculate the dynamic Onsager coefficients using the Kubo-Greenwood formula [START_REF] Pozzo | Thermal and electrical conductivity of iron at Earth's core conditions[END_REF],

L lm (ω) = (-1) l+m 2πe 2 2 3m 2 e ωV k n i,j=1 3 α=1 w(k)[f (ε i,k ) -f (ε j,k )] × | Ψ j,k |∆ α |Ψ i,k | 2 [ε j,k -µ] l-1 [ε i,k -µ] m-1 δ(ε j,k -ε i,k -ω) (3.9)
where e and m e are the electron charge and mass respectively, is the reduced Plank's constant, V is the cell volume and n is the number of bands used in the simulations, α sum runs over the three spatial directions, Ψ i,k and ε i,k are the Kohn-Sham wavefunctions and associated eigenvalues at the corresponding k-point, f (ε i,k ) is the Fermi factor. The δ function is represented by a Gaussian with width equal to the average spacing between the eigenvalues. w(k) is the weight of the particular k-point in the Brillouin zone.

The electrical (σ(ω)) and thermal (k(ω)) conductivity are obtained as,

σ(ω) = L 11 (ω), (3.10)
and

k(ω) = 1 e 2 T (L 22 (ω) - L 12 (ω)L 21 (ω) L 1 1ω
).

(3.11)

The σ 0 and k 0 are the respective values of σ(ω) and k(ω) in the limit ω → 0. In order to verify our implementation, we conducted simulations at 328 GPa and 6350 K for liquid iron using exactly the same parameters as [START_REF] Pozzo | Transport properties for liquid silicon-oxygen-iron mixtures at Earth's core conditions[END_REF]. At these conditions, the density of liquid iron is 12.95 g/cm 3 . We obtain an agreement better than 1 % for both types of conductivity, when compared to the values reported in [START_REF] Pozzo | Transport properties for liquid silicon-oxygen-iron mixtures at Earth's core conditions[END_REF].

Results and discussion

The critical point

We perform first-principles molecular dynamics simulations in the 3000-15000 K temperature range and densities below 8 g/cm 3 . This regime is characteristic of the aftermath conditions of giant impacts (Canup, 2004a). We compute the pressure dependence of the density along several isotherms (Fig. 3.1). We use a third-order polynomial expansion of the pressure as a function of density to identify the liquid spinodal and the position of the critical point, as detailed in the methodology. For iron, we identify a liquid spinodal point for all isotherms up to 9000 K. Along this latter isotherm the minimum pressure corresponding to the liquid spinodal is obtained at 2.40 g/cm 3 . At 9000 K, we extend the simulations towards even lower densities, which allows us to observe also a maximum along the pressure-density curve. This corresponds to the gas spinodal, lying at 1.85 g/cm 3 . Starting with the Pressure (GPa)

Figure 3.1: Variation of pressure as a function of density for iron along several isotherms. Along a given isotherm below the critical temperature, with volume expansion, the pressure may decrease to reach negative values. These negative pressures indicate the presence of hydrostatic tension in the system. According to the classic nucleation theory [START_REF] Karthika | A Review of Classical and Nonclassical Nucleation Theories[END_REF], the first-order transitions need to overcome energy barriers due to the surface energy, which prevents the formation of the thermodynamic stable phase. Therefore, this stage is thermodynamically metastable but mechanically stable. The minimum of the pressure marks the liquid spinodal (solid symbols). Joining the spinodal points yield the spinodal line (the black solid line). At densities lower than that of the liquid spinodal the pressure starts to increase until it reaches a maximum, which marks the gas spinodal, as shown in the left inset figure. At densities between the two spinodal lines a twophase mixture coexists. Because of technical computing limitations we compute the gas spinodal only at temperatures close to the critical one, that is 9000 K. Above the critical temperature (9350 K) the pressure decreases continuously with decreasing density, but does not show any minima or maxima. We obtain the critical point to be in the range 1.85-2.40 g/cm 3 and 9000-9350 K (black empty rectangle). The right inset shows comparisons of the critical point, between our estimate and the ones inferred from experiments [START_REF] Fortov | Shock waves and equations of state of matter[END_REF][START_REF] Medvedev | Wide-range multiphase equation of state for iron[END_REF].

9350 K isotherm the pressure varies monotonically without any local minimum or maximum; this is characteristic of the supercritical state. Therefore, the position of the CP is bracketed by the two spinodal lines, which intersect in the CP itself, and by the last isotherm with minima and maxima and the first isotherm with monotonical pressure variation. For iron, using the results of our simulations we predict that the CP lies in the 1.85-2.40 g/cm 3 , and 9000-9350 K range (Fig. 3.1). These values correspond to pressure of 4-7 kbars. The coordination number calculated by counting the number of atoms within the first minimum of g(r) from a central atom. There is a weak dependence of the coordination number with the temperature. However, its density dependence is very strong below 7.75 g/cm 3 , then saturate at a value of 13 up to 6000 K and 13.3 g/cm 3 [START_REF] Alfè | Structure and dynamics of liquid iron under Earth's core conditions[END_REF][START_REF] Vočadlo | First principles calculations on crystalline and liquid iron at Earth's core conditions[END_REF], which is the condition at Earth's inner-outer core boundary.

Static

Several typical pair distribution functions (PDF) calculated from the MD simulations are shown in Fig. 3.2(a). As the periodic boundary condition is used, the maximum correlation length is limited to half of the cell length. All PDFs start at zero, and continue as such up to a certain distance that defines the exclusion radius. This is caused by the strong repulsive force between any two atoms that prevents atoms from staying too close to each other. In the iron simulations, the exclusion radius depends weakly on the density and temperature; it decreases from 1.9 Å at 3000 K and 7.75 g/cm 3 to 1.8 Å at 12000 K and 0.37 g/cm 3 indicating a decline of repulsive force over this temperature and density regime.

The first peak in PDF, usually considered as a good approximation for the average bond distance in the fluid, occurs at 2.5 Å, with a height of 2.0 at 3000 K and 7.75 g/cm 3 . These values remain relatively constant over a large range of temperature and density conditions, but the peak broadens with an increasing of temperature due to thermal activation. Previous ab initio molecular dynamics simulations at 1850 K and around 6.95 g/cm 3 also reported the main peak position of the PDFs at around 2.5 Å [START_REF] Marqués | Ab initio study of the structure and dynamics of bulk liquid Fe[END_REF].

The first minimum of the PDF marks the end of the first coordination sphere. At 3000 K and 7.75 g/cm 3 , this lies at 3.3 Å and remains relatively constant over a large range of temperature and density conditions. We observe a second peak (the second coordination shell) at 4.5 Å, which again broadens with an increasing of temperature.

Above 7500 K and below 2.70 g/cm 3 , the positions of the first minimum and of the second peak become less obvious. And at 12000 K and 0.37 g/cm 3 , the PDF decays steadily to a value close to one without an apparent minimum point. Due to the small cell size used here, we could not observe the long-distance behaviour in PDF beyond 6 Å.

The integral of the PDF up to the first minimum gives the coordination numbers (CNs). The Fe-Fe coordination as a function of density and temperature has been shown in Fig. 3.2(b). We observe a weak dependence of the coordination number with temperature. For example, at a density of about 6 g/cm 3 , the coordination number decreases from about 11 at 3000 K down to about 9.5 at 7500 K, and by extrapolation down to about 8 at 12000 K. However, density has a more pronounced effect on the changes in coordination number. Along the 3000 K isotherm the coordination number increases from about 9.5 at about 6 g/cm 3 up to 13.5 at 8 g/cm 3 . At 3000 K and 7.75 g/cm 3 (about 5 GPa in pressure) the coordination is around 13.6, indicating a very close-packed liquid. Actually the value of 14 seems to be a saturation value for coordination in liquid iron at high pressures and densities, i.e. at 6000 K and 13.3 g/cm 3 , the conditions at inner-outer core boundary inside the Earth, as was shown in other previous ab initio MD studies [START_REF] Alfè | Structure and dynamics of liquid iron under Earth's core conditions[END_REF][START_REF] Vočadlo | First principles calculations on crystalline and liquid iron at Earth's core conditions[END_REF]. The change in slope around 7.75 g/cm 3 implies that probably the compression mechanism above this density is likely to be caused by the change in the second coordination shell. With decreasing density along the same isotherm, the Fe-Fe coordination decreases to around 10.8 at a density of 5.81 g/cm 3 . In the supercritical region, at 12000 K, the Fe-Fe coordination drops from 6.3 at 2.7 g/cm 3 density to 1.6 at 0.37 g/cm 3 . A change in slope of the density-coordination dependence occurs around 1.46 g/cm 3 along 12000 K. At the critical point, the Fe-Fe coordination number is around 6.

Speciation

The coordination number as obtained from PDF only gives an average information on the number of atoms that surrounds a central iron atom. The distribution of the bond lengths and the lifetimes of atomic clusters are still missing. Consequently, we also use a geometric criterion to analyse the interatomic connectivity at every configuration generated during our simulations. The geometric cut-off length for two atoms to be considered connected, i.e. bonded, is chosen to be the radius of the first coordination sphere that is the first minimum of PDF. Then the size of any given iron coordination polyhedron is defined as the number of atoms within the range of the first minimum of PDF. Once the individual atoms forming a given coordination polyhedron are identified, we monitor the time they obey the geometrical criterion for bonding to define the polyhedron lifetime. Fig. 3.3(a) shows the coordination of atoms at different conditions. In general, we observe, as expected, a depolymerization of the fluid with decreasing temperature. At 3000 K, from 7.97 g/cm 3 to 6.78 g/cm 3 , the iron is highly polymerized with few or no isolated Fe atoms. With further decrease of density, we observe the presence of isolated Fe atoms with an increasing abundance from 0.02 at 6.78 g/cm 3 to 0.15 at 5.81 g/cm 3 indicating a gradual depolymerization. Along the 4000 K isotherm, the abundance of isolated iron atoms increases rapidly from 0.07 at 6.78 g/cm 3 to 0.50 at 5.14 g/cm 3 . Above 7500 K and below 5.46 g/cm 3 , we observe a nearly constant concentration of 0.85 of isolated iron atoms, irrespective of the temperature. The lifetime of isolated Fe atoms is shown in Fig. 3.3(b); in general, the increase of temperature and the decrease of density result in rising their lifetime. In the supercritical regime, at 12000 K, the lifetime of isolated atoms increases linearly with decreasing the density. At the lowest density of 0.37 g/cm 3 calculated in the present study, the lifetime is around 450 femtoseconds. We also observed some small Fe 2 dimers starting to appear from 5000 K and 6.60 g/cm 3 , as illustrated in Fig. 3.3(c). There is a weak temperature dependence of the Fe 2 abundance, with a major role played by the density. The abundance of Fe 2 clusters increases considerably in the 5.46 to 4.32 g/cm 3 density range. Below 5.46 g/cm 3 density, the abundance of Fe 2 fluctuates around 0.11 until 1.08 g/cm 3 , then increases again to 0.14 at further decreasing density to 0.37 g/cm 3 . We only find a small amount of iron trimers (Fe 3 ) at 12000 K. The lifetime of Fe 2 and Fe 3 clusters (Fig. 3.3(d)) is below 100 femtoseconds, shorter than that of Fe. Our results suggest the critical iron fluid is mainly atomic with the Fe dimers being the next most abundant species. Velocity autocorrelation function a 7.75 g/cm 3 6.6 g/cm 3 5.81 g/cm 3 5.81 g/cm 3 5.81 g/cm 3 5.81 g/cm 3 4.57 g/cm 3 3.86 g/cm 3 3.86 g/cm 3 3.86 g/cm 3 2.98 g/cm 3 1.46 g/cm 3 0.73 g/cm 3 0.37 g/cm 3 ideal gas Intensity 3000 K, 7.75 g/cm 3 3000 K, 6.6 g/cm 3 3000 K, 5.81 g/cm 3 4000 K, 5.81 g/cm 3 5000 K, 5.81 g/cm 3 7500 K, 5.81 g/cm 3 7500 K, 4.57 g/cm 3 7500 K, 3.86 g/cm 3 8750 K, 3.86 g/cm 3 10000 K, 3.86 g/cm 3 4(a). At 3000 K and 7.75 g/cm 3 , we observe a broad peak from 0-400 cm -1 , which is caused by breathingtypes vibrations of the first coordination shell.

Velocity autocorrelation function

As a time-dependent correlation function, the velocity autocorrelation function (VAF) not only reveals the underlying dynamical processes like diffusional and vibrational motions operating in an atomic or molecular system at a microscopic scale, but also shows a direct connection with the macroscopic properties like the diffusion coefficient. These facts offer a good reason to study VAF, whose straightforward calculation details are provided in Section 3.2.5. We plot the overall trend from 3000 K and 7.75 g/cm 3 to 10000 K and 3.86 g/cm 3 in Fig. 3.4(a). At 7.65 g/cm 3 and 3000 K, VAF decays to reach a first minimum at a value of 80 fs. The formation of this first minimum is generally attributed to the cage effect formed by the nearest neighbors, which exert a restoring opposite force on the central atom when this encounters the cage during its vibration. There are two more minima at larger correlation times. Globally, the three observed minima represent different vibrational modes. With further decreasing density to 5.81 g/cm 3 at 3000 K, as the atoms become increasingly spaced, only two minima are found, and their positions are shifted to longer correlation times (i.e. low frequency) due to a decrease of the interaction strength (also see Fig. 3.5). With increasing temperature at fixed density of 5.81 g/cm 3 , the minimum becomes less obvious and disappears at 7500 K where a single exponential decay is present, indicating the cage structure can no longer hold the Fe atoms due to the large thermal velocity at high temperature. The VAF decay is weakly dependent on temperature, but appears at relatively similar densities.

To illustrate the dependence of the VAF with the density, we show its variation along the 12000 K isotherm in Fig. 3.4(b). The decay of VAF becomes faster at high densities since collisions occur more frequently due to the relatively short distances between atoms, which break correlations. The decay time of VAF at 12000 K and 0.37 g/cm 3 is around 2500 fs, suggesting that at these conditions, only 12 independent configurations are sampled during a 30 ps simulation. The situation would become even worse when the temperature and density are below 12000 K and 0.37 g/cm 3 since the correlation time increases with decreasing temperature and density. [START_REF] Reed | Hybrid molecular dynamics: an approach to low density simulations[END_REF] found the statistical inefficiency, which represents the average time needed to sample two statistically uncorrelated configurations during an MD simulation, is almost linearly scaling with the inverse of density. For the low-density phase, an intrinsic problem is that there are two times scales: the inter-collision time and collision duration. As the collision duration is very short, the atoms move in space for a very long time before the next collision occurs. A common method to speed up the correlation time is to use a larger time step. But a very large time step would cause missing collisions and thus decreasing accuracy of sampling of the configurational space. [START_REF] Reed | Hybrid molecular dynamics: an approach to low density simulations[END_REF] proposed a hybrid algorithm that combines the time-driven simulation with event-driven simulation, which is frequently used in the hard-sphere models. This may provide a promising solution to the time scale problems. However, many questions need to be answered before its practical implementations combined with density functional theory will become available, especially on how to define a collision radius. The inset shows a log-log zoom-in plot for short times. The inset outlines a slope change that separates the two characteristic regions of the diffusion: ballistic, for times less than about 100 fs, and diffusive, at higher times.

Diffusion

The mean square displacements (MSD) measure the average distance that atoms move during a reference time window. We show in Fig. 3.6 four examples that cover the whole studied temperature range from 3000 to 12000 K and density range from 7.75 to 0.37 g/cm 3 . The total diffusion process can be divided into two stages, as seen by the change of slope in the log-log plot in the inset of Fig. 3.6. The first part, on the order of about 100 fs, is the ballistic part dominated by the free particle motion before collision. The second part is dominated by the diffusion of the particles; in this region the square displacement linearly scales as the function of time, where the slope of the MSD yields the diffusion coefficient. Figure 3.7: The temperature and density dependence of the self-diffusion coefficients of fluid iron. The solid circles and black crosses represent diffusivities determined from velocity autocorrelation function and mean-squared displacement, respectively. The green squares and black triangles denote several self-diffusion coefficients of liquid iron at high densities and temperatures previously reported in the literature [START_REF] Alfè | Structure and dynamics of liquid iron under Earth's core conditions[END_REF][START_REF] Posner | Structural changes and anomalous self-diffusion of oxygen in liquid iron at high pressure[END_REF].

We also calculate the diffusion coefficient by integrating the velocity autocorrelation function, as shown in Fig. 3.7. The diffusivity increases by two orders of magnitude from about 1.0 × 10 -8 m 2 s -1 at 7.75 g/cm 3 and 3000 K to about 1.3 × 10 -6 m 2 s -1 at 0.37 g/cm 3 and 12000 K. At high densities and low temperatures, the diffusion coefficients correspond to distances travelled by the iron atoms on the order of about 2.5 cm in 24 hours, which would ensure highly efficient chemical equilibration close to the surface of the iron droplets. At low densities and high temperatures, the diffusion coefficients correspond to displacements on the order of 0.5 m in 24 hours, which would contribute significantly to the chemical equilibration inside large-scale fluid volumes. As expected, the increasing temperature and decreasing density would increase the diffusivity. There is an excellent agreement between the values of the diffusion coefficients obtained using the two methods (the slope of the MSD and the integration of the VAF). Our values also compare well with previous theoretical studies in liquid iron [START_REF] Alfè | Structure and dynamics of liquid iron under Earth's core conditions[END_REF][START_REF] Posner | Structural changes and anomalous self-diffusion of oxygen in liquid iron at high pressure[END_REF], at higher pressures and lower temperatures. [START_REF] Assael | Reference Data for the Density and Viscosity of Liquid Aluminum and Liquid Iron[END_REF], which is close to our 3000 K isotherm. The black solid circles are the viscosity at 3000 K and 7.75 g/cm 3 determined from the Stoke-Einstein relation, which is in a fair agreement with the prediction by the Green-Kubo method.

Viscosity

The viscosity of iron was computed previously only in its high-density liquid state [START_REF] Alfè | Structure and dynamics of liquid iron under Earth's core conditions[END_REF][START_REF] De Wijs | The viscosity of liquid iron at the physical conditions of the Earth's core[END_REF] and applied to the dynamics of the liquid outer core of the Earth. We obtain the viscosity from the self-correlation of the shear components stress tensor. As a general rule this method requires much longer simulations times than are needed to obtain the other properties. We can improve the quality of the results by averaging over all the three shear components of the stress tensor. Fig. 3.8(a) displays the stress autocorrelation function (SAF) for fluid iron at 3000 K and 7.75 g/cm 3 . After 0.15 ps, SAF decays to zero with a still large noise, about 5% compared with the value at t = 0. The inset in Fig. 3.8(a) shows the viscosity integral as a function of time. After 0.4 ps, the shear viscosity fluctuates within the range 3.5 × 10 -3 -5 × 10 -3 Pa s.

We present the variation of viscosity as a function of density along several isotherms in Fig. 3.8(b). The values of viscosity at conditions typical to the outer parts of the protolunar disk, i.e. low density and high temperature, are in the range 0.5 × 10 -3 to 4 × 10 -3 Pa s. They show a decreasing trend with decreasing density, but due to the large scatter in values, we cannot identify any reliable trend with temperature. The values that we obtain for conditions specific to the protolunar disk are one order of magnitude smaller than in the Earth's outer core, estimated to be on the order of 2 × 10 -3 -15 × 10 -3 [START_REF] Alfè | Structure and dynamics of liquid iron under Earth's core conditions[END_REF][START_REF] De Wijs | The viscosity of liquid iron at the physical conditions of the Earth's core[END_REF]. [START_REF] Assael | Reference Data for the Density and Viscosity of Liquid Aluminum and Liquid Iron[END_REF] compiled viscosity data for liquid iron from 1850 K and 7.0 g/cm 3 to 2500 K and 6.40 g/cm 3 . We find the viscosity at 3000 K and 6.40 g/cm 3 of 2.1 × 10 -3 Pa s to be comparable to the corresponding value at 2500 K and 6.40 g/cm 3 of 2.276 × 10 -3 Pa s [START_REF] Assael | Reference Data for the Density and Viscosity of Liquid Aluminum and Liquid Iron[END_REF].

At 3000 K and 7.75 g/cm 3 we perform a further check and compute the viscosity using the Stoke-Einstein relations, which connects diffusion coefficient to estimate the viscosity:

η = k B T 2πaD , (3.12)
where η is the viscosity, k B is the Boltzmann constant, T is temperature, D is the diffusion coefficient and a is the effective atomic diameter chosen to be the radius of the first peak in the radial distribution function from Fig. 3.2. This method has been used in previous work to determine the viscosity of liquid iron at Earth's core conditions [START_REF] De Wijs | The viscosity of liquid iron at the physical conditions of the Earth's core[END_REF]. By applying this method, the estimate viscosity is 2.8 × 10 -3 Pa s, where a fair agreement with Green-Kubo method is achieved.

Electrical and thermal conductivity

We used the Kubo-Greenwood formalism [START_REF] Greenwood | The Boltzmann Equation in the Theory of Electrical Conduction in Metals[END_REF][START_REF] Kubo | Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems[END_REF] to determine the electrical and thermal conductivities. In this method, the electron-ion collisions are accurately described, while the electron-electron scattering and is neglected and its effect remains unclear [START_REF] Dufty | On the Kubo-Greenwood model for electron conductivity[END_REF]. All calculated electrical, thermal conductivity and the Lorentz constant of pure fluid iron as a function of density and temperature are shown in Fig. 3.9. Along 3000 K isotherm, the electrical conductivity decreases from 1.2 × 10 6 Ω -1 m -1 at 7.76 g/cm 3 to 0.8 × 10 6 Ω -1 m -1 at 5.81 g/cm 3 , while the thermal conductivity declines from 80 Wm -1 K -1 to 40 Wm -1 K -1 . Along other isotherms, the electrical and thermal conductivity decrease monotonically upon a decreasing density as well.

At a density that is beyond 6.19 g/cm 3 , the electrical conductivity decreases with temperature, as expected for a typical metal. The increasing temperature would generally enhance the electron-ion collisions causing a decrease of the relaxation time and electrical conductivity. Below 6.19 g/cm 3 , the trend of electrical conductivity as a function of temperature is reversed, where the electrical conductivity increases with temperature. It indicates the fluid iron at these conditions is more like a semiconductor, where the temperature dependence may be as a result of the increasing concentration of carriers due to thermal excitation. The observed change of temperature dependence around 6.19 g/cm 3 suggests a possible metal-non-metal transition and has also been observed in the other expanded metals like aluminium [START_REF] Recoules | Electrical conductivity of hot expanded aluminum: Experimental measurements and ab initio calculations[END_REF].

Our calculated Lorentz number (Fig. 3.9(c)) almost keeps constant around 2.44 × 10 -8 W Ω K -2 over the whole temperature and density range in this study. Compared with [START_REF] Korell | Paramagneticto-Diamagnetic Transition in Dense Liquid Iron and Its Influence on Electronic Transport Properties[END_REF], it is no surprise that a very good agreement is reached as the same method and DFT package are used.

Hugoniot lines

The behavior of materials under shock waves can be described using the Rankine-Hugoniot equations. These equations relate the density, pressure, and internal energy after shock to the initial state by,

E -E 0 = (P + P 0 )(V 0 -V )/2 (3.13)
where E, P , V are the internal energy, pressure and volume, respectively. And the 0 subscript denotes the initial state. The MD simulations that we performed at various isotherms contain all the information needed to build the Hugoniot EOS (see Fig. We have calculated the principal Hugoniot line with a initial condition of 300 K and 7.874 g/cm 3 , which compares well with theoretical simulation results by [START_REF] Sjostrom | Quantum molecular dynamics of warm dense iron and a five-phase equation of state[END_REF] and shock experiments.

We consider two representative initial states. The first case has iron at 1 GPa and 1500 K, conditions similar to what we could expect to have in small planetesimals.

For these conditions that we call warm Hugoniot the EOS intercepts the iron melting curve at 130 GPa. Previous shock experiments [START_REF] Chen | High Pressure and High Temperature Equationof-State of Gamma and Liquid Iron[END_REF] on facecentered cubic iron with an initial condition of 1570 K and 1 bar show that the Hugoniot intercepts the iron melting curve at 80 GPa. In the second case we consider the initial state at 40 GPa and 4000 K, which may be representative for the state of the core in Mars-sized impactors (Canup, 2004a). At these conditions iron is already molten. Fig. 3.11(a) shows the computed Hugoniot lines for these two cases. During shocks the temperature can easily reach thousands of degrees and the pressures hundreds of GPa. These would be typical conditions for the core state during giant impacts.

The actual amount of vaporization after an impact depends on entropy achieved. The entropy at the boiling point was estimated at 15.84 k B /atom at 3100 K and 1 bar [START_REF] Kraus | Impact vaporization of planetesimal cores in the late stages of planet formation[END_REF]. If the peak shock conditions during impact exceed this entropy value, then the onset of vaporization may take place and part of the shocked material can vaporize upon release and cooling [START_REF] Ahrens | Shock melting and vaporization of lunar rocks and minerals[END_REF].

We compute the entropy of the liquid iron along the two Hugoniot lines and at the spinodal points (Fig. 3.11) as a function of temperature from the vibrational spectra (see Section 3.2.6). For peak shock conditions at 15000 K, we estimate that entropy can reach 19.1 k B /atom and 18.6 k B /atom along the warm and the hot Hugoniot curves respectively. At these conditions the entropy is high enough to result in partial vaporization of the iron core. Our results show that above 7500 K, the entropy along the warm Hugoniot is less than along the hot Hugoniot. The entropy difference between these two Hugoniot lines is relatively small (0.5 k B /atom) in the 7500 K to 15000 K range. If we relate entropy to the peak shock pressure, based on the computed entropy along the warm Hugoniot line, we find that the shock pressure required to reach the onset of vaporization upon release and cooling is 312 GPa. This is less than previous estimates of 390 GPa [START_REF] Kraus | Impact vaporization of planetesimal cores in the late stages of planet formation[END_REF]. Along the hot Hugoniot, the onset vaporization pressure is 365 GPa. This is only slightly higher than that of warm Hugoniot.

Vaporization of small planetesimals

The onset of core vaporization can easily be reached in case of impacts of small planetesimals, like the ones that might have occurred either during the first stages of formation of the solar system, or during the late veneer. As a model example we consider a differentiated planetesimal with a mantle made of enstatite (MgSiO 3 ) and a core made of iron; we set the core-mantle boundary at 1 GPa and 1500 K [START_REF] Raymond | Building the terrestrial planets: Constrained accretion in the inner Solar System[END_REF]. We approximate the shock wave as a planar wave [START_REF] Melosh | Planetary Surface Processes[END_REF] traveling through the two layers. However, this yields a simplified estimate of the peak pressure and does not thoroughly describe the pressure distributions in these bodies. When the impact occurs, shock waves travel through the silicate layers of the two bodies. At the core-mantle boundary, because of the density contrast between silicates and iron, the shock wave is partly reflected, going backward into the mantle, and partly transmitted, going forward into the core. Assuming a steady shock in a model MgSiO 3 -based mantle [START_REF] Militzer | Equation of state calculations of hydrogen-helium mixtures in solar and extrasolar giant planets[END_REF] and in the iron core, the impedance match method allows us to determine the properties of the reflected wave in the mantle and the transmitted wave in the core [START_REF] Forbes | Shock Wave Compression of Condensed Matter[END_REF]. Fig. 3.12 illustrates the propagation of the shock wave according to this model through the [START_REF] Sjostrom | Quantum molecular dynamics of warm dense iron and a five-phase equation of state[END_REF]) that starts at ambient conditions and the melting curve of iron [START_REF] Bouchet | Ab initio equation of state of iron up to 1500 GPa[END_REF] are shown for reference only. For a smaller impactor like Moon, we approximate the initial state conditions as 1 GPa and 1500 K; for a large impactor like Mars, the initial state is set at 40 GPa and 4000 K. The Hugoniot lines cross because the gains in temperature and pressure is not linear with respect to changes in initial conditions. The shaded area represents estimated temperature gradients ranges in the metallic cores of the different objects involved in the impact [START_REF] Antonangeli | Toward a mineral physics reference model for the Moon's core[END_REF][START_REF] Hirose | Composition and State of the Core[END_REF][START_REF] Stewart | Mars: A New Core-Crystallization Regime[END_REF]. (b) Computed entropy along the two Hugoniot lines. The star indicates the experimentally estimated entropy of boiling liquid iron at 3000 K, marked also by the dashed line. .12: Graphical representation of the impedance match analysis for the core of the small impactor. The shock wave generated at the moment of the impact travels through the impactor's mantle assumed to be homogeneous. The shock state at the core-mantle is given by the MgSiO 3 principal Hugoniot line (green line) estimated based on previous ab initio simulations [START_REF] Militzer | Equation of state calculations of hydrogen-helium mixtures in solar and extrasolar giant planets[END_REF]. More recent experimental Hugoniot data points are represented by crosses [START_REF] Fratanduono | Thermodynamic properties of MgSiO 3 at super-Earth mantle conditions[END_REF]. At the core-mantle boundary the shock wave splits in two opposite waves. The one traveling forward enters the core. The final state in the core is given by the intersection between the re-shocked MgSiO 3 Hugoniot (red line) and the iron Hugoniot (blue line) with an initial state at 1 GPa and 1500 K (the warm Hugoniot in Fig. 3.11). The inset shows the fraction of iron that would vaporize from the corresponding impactor's core as a function of the impact velocity.

planetesimal.

As the shock proceeds through the mantle of the impactor this can lead to partial fragmentation. During this process the mantle fragments and can detach from the core leaving behind bare fragments of shocked core. In the post-shocked state, fragments of planetesimal core without mantle confinement can undergo an isentropic release into vacuum. In this case, if the entropy reached during the shock is high enough, the core may partly vaporize; otherwise it will remain liquid and accrete to the impacted body, or escape gravitationally and eventually crystallize.

For heads-on collisions of small planetesimals, the impact velocity required to onset vaporization is around 11.5 km/s. With an impact velocity around 13.5 km/s, peak pressure and temperature reach 450 GPa and 15000 K (Fig. 3.12). At these conditions about 22 % of iron would vaporize. In contrast, oblique impacts greatly reduce the peak pressure, due to a sin(θ) factor where θ is the obliquity (Pierazzo and Melosh, 2000a). For the maximum frequency impact angles of 45 • (Pierazzo and Melosh, 2000b), the velocity threshold to onset vaporization increases to 15.3 km/s. With mean impact velocities at 14.5 km/s and median impact angles at 40 • , around 70% of the impacts of N-body simulations [START_REF] Raymond | Building the terrestrial planets: Constrained accretion in the inner Solar System[END_REF] yield velocities larger than our threshold. This suggests that core vaporization is a common process during planetary formation. During the collision with the Earth's mantle, the core of the incoming planetesimal can be efficiently mixed into the molten silicate pond locally produced by the impact itself or into a pre-existing larger magma ocean. If such impacts happen after the Earth's core formation, this process would then increase the amount of highly siderophile elements that is seen today trapped and dispersed into the Earth's mantle.

Vaporization during giant impacts

In the case of giant impacts, the geometry effect plays an important role in controlling the shock peak conditions. As the validity of the impedance matching method is limited to the impacts where the lateral dimension of the impactor is small compared to the distance the shock wave has propagated [START_REF] Melosh | Planetary Surface Processes[END_REF], it has only a limited applicability. However, the entropic and pressure criteria for vaporization still hold. For impacts with Mars-sized bodies, because of the hotter initial state of their cores, our simulations suggest a vaporization pressure of only 312 GPa. This is again smaller than previous estimations by [START_REF] Kraus | Impact vaporization of planetesimal cores in the late stages of planet formation[END_REF] suggesting that even more iron will be vaporized than previously thought. However, the amount of iron that can be vaporized depends also on the local pressure conditions as the process of the impact itself takes its due course. As the predicted pressure thresholds for vaporization can be easily reached, a large amount of iron receives enough entropy to vaporize. The entropy threshold can even be easier exceeded due to the entropy gain after the first and secondary shocks and the conversion of gravitational potential energy to internal energy [START_REF] Carter | The Energy Budgets of Giant Impacts[END_REF][START_REF] Nakajima | Melting and mixing states of the earth's mantle after the moon-forming impact[END_REF]. Once again during this process the confinement of core fragments by the surrounding .13: Isentropic release of the iron core starting from various peak conditions attained during large impacts. The shocked material gains entropy that is conserved during the release. If this entropy is larger than that of the boiling point then partial vaporization may occur. We compute the entropy along the spinodal and the Hugoniot lines at different temperatures, which allows us to map the entropy increase during various impact scenarios. We find that the entropy of the critical point, at 21.8 k B /atom, is reached for impacts with peak pressures of 605 GPa in the case of small planetesimals or for peak pressures of 825 GPa in the case of Marssized objects. Peak pressures of 312 GPa are enough to provide entropy higher than the entropy at boiling of iron at 1 atm, i.e. 15.84 k B /atom [START_REF] Kraus | Impact vaporization of planetesimal cores in the late stages of planet formation[END_REF]. These conditions can be easily exceeded during giant impacts between the proto-Earth and Theia, but could also be reached in almost half of the impacts with planetesimals during the late veneer. However, the liquid-vapor dome is reached only if the density, and hence the pressure, is allowed to decrease sufficiently. This can happen if the mantle is stripped away when fragments of the core are allowed to decompress without the mantle confinement. mantle may prohibit the isentropic expansion and thus the vaporization (Fig. 3.13). But if during the impact parts of the mantle are detached from fragments of the core [START_REF] Nakajima | Melting and mixing states of the earth's mantle after the moon-forming impact[END_REF], then during the isentropic release these fragments will then undergo partial vaporization. Part of the vapor will remain in the outer part of the disk and eventually condense to form the Moon's core while the rest will fall into the central body and mix into the magma ocean. [START_REF] Kendall | Differentiated planetesimal impacts into a terrestrial magma ocean: Fate of the iron core[END_REF] suggest fully mixing between impactor's core and proto-Earth Magma Ocean can be achieved for iron blobs that are less than 100 km across. Consequently, core fragmentation, promoted by partial vaporization during release, will enhance equilibrium and/or mixing between the impactor's core and the molten silicates, on a larger degree than the previous estimations of hydrodynamic simulations, which generally predicted the impactor's core directly merge into the Earth's core [START_REF] Canup | Forming a Moon with an Earth-like Composition via a Giant Impact[END_REF][START_REF] Ćuk | Making the moon from a fast-spinning earth: a giant impact followed by resonant despinning[END_REF]. Then the mixing process can easily explain the recent W-isotope data, which require at least 30 % core-mantle equilibration in the aftermath of the giant impact [START_REF] Nimmo | Tungsten isotopic evolution during latestage accretion: Constraints on Earth-Moon equilibration[END_REF][START_REF] Rudge | Broad bounds on Earth's accretion and core formation constrained by geochemical models[END_REF][START_REF] Touboul | Tungsten isotopic evidence for disproportional late accretion to the earth and moon[END_REF].

Conclusions

In this chapter, we perform ab initio molecular dynamics to determine the position of the critical point of iron, and to characterize the fluid iron over a wide density and temperature range. Based on our calculations, we predict the critical point of iron to be in the 9000-9350 K temperature range and 1.85-2.40 g/cm 3 density range, corresponding to a pressures range of 4-7 kbars. We find that the lowdensity fluid is highly depolymerized, its structure dominated by isolated atoms and dimers. The coordination number increases strongly with density below 7.75 g/cm 3 , then saturates at a value of 14 that persists at least up to 6000 K and 13.3 g/cm 3 , which is the condition at the Earth's inner-outer core boundary. The viscosity of iron in the outer parts of the protolunar disk is also extremely low, on the order of 10 -3 Pa s, which is one order of magnitude smaller than the value in the Earth's liquid outer core [START_REF] De Wijs | The viscosity of liquid iron at the physical conditions of the Earth's core[END_REF]. A low viscosity implies a limited role in the global energy budget of the disk, a term oftentimes neglected in magneto-hydrodynamic simulations, and a possible presence of turbulence as the primary mechanism of transporting mass in the disk. Finally, the computed electronic and thermal conductivities decrease with the density, suggesting a gradual reduction of the metallic character. We do not observe any discontinuity in the density dependence of the electric conductivity. The only visible effect is that the electric conductivity falls more rapidly as the density decreases below 5.14 g/cm 3 .

The determination of the Hugoniot lines and our estimations of the amounts of entropy gained during giant impacts show that the core of Theia underwent partial vaporization. This would easily explain the recent W-isotope data which requires at least 30% core-mantle equilibration in the aftermath of the giant impact (Nimmo

Introduction

In the previous chapter, we used the spinodal line to estimate the position of the critical point of iron. As discussed by [START_REF] Binder | Beyond the van der waals loop: What can be learned from simulating lennard-jones fluids inside the region of phase coexistence[END_REF], the spinodal line is due to the presence of an interface between the liquid and vapor phases as a result of the finite size effect. In this chapter, we present the Gibbs ensemble (GE) Monte Carlo (MC) method, which was first developed by [START_REF] Panagiotopoulos | Direct determination of phase coexistence properties of fluids by monte carlo simulation in a new ensemble[END_REF] and is able to determine the liquid-vapour equilibrium without an explicit interface between two phases. In this method, we simulate two independent phases by two boxes at the same time and allow them to exchange volume and particles to reach the chemical equilibrium. Following this way, it directly mimics the macroscopic equilibrium at a microscopic level but avoids the interface issues allowing for a more reliable determination of the phase equilibrium. At this stage, we limit ourselves to the N V T version of the Gibbs ensemble where the number of particles (N ), the total volume (V ) of two boxes and the temperature (T ) is fixed. As this method relies on a reasonable number of particle exchanges to achieve the same chemical potential, it is not very useful for studying equilibria involving high-density phases.

In the classical Gibbs ensemble scheme [START_REF] Panagiotopoulos | Direct determination of phase coexistence properties of fluids by monte carlo simulation in a new ensemble[END_REF], the potential energy evaluated from a force field is used in the acceptance criterion for a trial move. This can be extended to include ab initio simulations, where the internal energy of the electron-ion system plays a similar role as the potential energy in the classical scheme. In the pioneering work of [START_REF] Mcgrath | Toward a monte carlo program for simulating vapor-liquid phase equilibria from first principles[END_REF], the GEMC method has been combined with KS-DFT -often referred as ab initio (aiGEMC) -to improve its accuracy. Ever since different groups have applied the aiGEMC technique to calculate the liquid-vapor equilibrium line of water (McGrath et al., 2006a,b;[START_REF] Schienbein | Liquid-vapor phase diagram of rpbe-d3 water: Electronic properties along the coexistence curve and in the supercritical phase[END_REF]Marx, 2018), methanol, methane (McGrath et al., 2011) and argon [START_REF] Goel | Predicting vapor liquid equilibria using density functional theory: A case study of argon[END_REF]. All above materials are insulators and have a critical point that is below 700 K. Due to the relatively large band gaps, the thermal excitation of electrons in these systems is not significant, while it may play an important role in determining the phase equilibrium of metallic systems at high temperature. In order to capture the thermal excitations of electrons, it is necessary to make full use of FT-KS-DFT which explicitly includes the contribution of excited electronic states. Because of the extra electronic entropy term in FT-KS-DFT framework, a great care should be given to the electronic contribution as the excited states also play a role in the acceptance ratio. This problem is of considerable conceptual importance, and thus warrants a brief outline of statistical mechanics derivations in order to find an appropriate energy term that can be used in the acceptance criterion for a MC trial move. Here we chose the Monte Carlo method due to its simplicity for practical implementations. But the generalization to Gibbs ensemble molecular dynamics is straightforward (e.g. [START_REF] Palmer | Molecular dynamics implementation of the gibbs ensemble calculation[END_REF].

Statistical mechanics of the Gibbs ensemble

Partition function in the Gibbs ensemble

We start with Eq. 2.12, then perform the integral over the momenta {P } N analytically and renormalize the coordinates {R} N by the factor L ≡ V 1/3 leading to the standard configurational canonical N V T partition function Z config (N, L, T ) [START_REF] Balian | From Microphysics to Macrophysics: Methods and Applications of Statistical Physics[END_REF][START_REF] Huang | Statistical Mechanics[END_REF]:

Z(N, V, T ) = V N N !Λ 3N dB N i=1 ds i exp (-βF 0 (N, L, {s} N )) ≡ V N N !Λ 3N dB Z config (N, L, T ), (4.1)
where Λ dB = h/ √ 2πM k B T is the de Broglie wavelength, and s i ∈ [0; 1] ∀i, F 0 (N, L, {s} N ) is an effective potential energy that is equal to F 0 ({R} N ) as defined in Eq. 2.11. The Gibbs ensemble [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF][START_REF] Panagiotopoulos | Phase equilibria by simulation in the gibbs ensemble[END_REF][START_REF] Panagiotopoulos | Direct determination of phase coexistence properties of fluids by monte carlo simulation in a new ensemble[END_REF]) is a particular case of the canonical ensemble, made of two sub-systems of volumes V 1 and V 2 containing respectively N 1 and N 2 particles. The two sub-systems are in contact with the same thermostat and can exchange volume and atoms, but the sum N 1 + N 2 = N and V 1 + V 2 = V are conserved. This construction ensures that all the intensive quantities are equal when equilibrium is reached: temperature (by construction), pressure (P 1 = P 2 ) and chemical potential (µ 1 = µ 2 ). The corresponding partition function for the Gibbs ensemble is written as the product of the partition functions of each sub-system by Eq. ( 4.1) and sum over all volumes and particle distributions:

Z Gibbs = Ntot N 1 =0 Vtot 0 dV 1 Z 1 (N 1 , V 1 , T ) × Z 2 (N 2 , V 2 , T ) = Ntot N 1 =0 Vtot 0 dV 1 V N 1 1 V N 2 2 N 1 !N 2 !Λ 3(N 1 +N 2 ) dB ×Z config 1 (N 1 , L 1 , T ) × Z config 2 (N 2 , L 2 , T ). (4.2)
From Eq. 4.2, we can see that one microstate in the Gibbs ensemble is specified by six parameters

(N 1 , L 1 , {s} N 1 , N 2 , L 2 , {s} N 2 )
. It is also interesting to note that the partition function written as in Eq. 4.2 includes the possibility for one box to be empty or vanishing volume. Although counter-intuitive, these states, if explored by the GEMC algorithm, must be taken into account [START_REF] Smit | Computer simulations in the gibbs ensemble[END_REF].

Monte Carlo in the canonical ensemble

A direct computation of the integral appearing in Eq. ( 4.1) or Eq. (4.2) is in general not possible, but most of the time we only need to estimate ensemble averages such as:

A = 1 Z config (N, L, T ) N i=1 ds i A(N, L, {s} N ) exp (-βF 0 (N, L, {s} N ) . (4.3)
For such calculations specific methods have been developed among which is the famous MC importance sampling algorithm introduced by [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF]. As such approaches are fairly well-known, we will only give the key features for the sake of completeness and to introduce the notations that we will need when presenting our numerical results.

The Metropolis algorithm is a Markov chain Monte Carlo method that can draw a sequence of random samples from the probability distribution of

ρ({s} N ) = exp{-βF 0 (N, L, {s} N )} Z config (N, L, T ) (4.4)
without the need to know Z config (N, L, T ). For a Markov process, we need to define a transition probability P ({s } N |{s} N ), which is the probability of transitioning from the current state {s} N to any other state {s } N . The choice of P ({s } N |{s} N ) needs to meet the detailed balance,

ρ({s} N )P ({s } N |{s} N ) = ρ({s } N )P ({s} N |{s } N ), (4.5)
to assure that the Markov process asymptotically reaches a unique stationary distribution such that,

lim n→∞ P n ({s} N |{s } N ) = ρ({s} N ), (4.6)
which is independent of the initial state ({s } N ). If the initial state is not sampled from ρ({s} N ), we need to equilibrate the system and discard the initial part of simulations. [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] separate P ({s } N |{s} N ) into two parts as,

P ({s } N |{s} N ) = Q({s } N |{s} N )A({s } N |{s} N ) (4.7)
where Q({s } N |{s} N ) is the conditional probability of proposing a state {s } N given {s} N and chosen to be symmetric meaning Q({s

} N |{s} N ) = Q({s} N |{s } N ), A({s } N |{s} N ) is the acceptance ratio by A({s } N |{s} N ) = min(1, ρ({s } N ) ρ({s} N ) ). (4.8)
Then the general procedure for each MC step is:

• compute the potential energy (F 0 (N, L, {s} N )) of the system in the old state {s} N ,

• make a move proposal based on Q({s } N |{s} N ),

• compute the potential energy (F 0 (N, L, {s } N )) of the system in the new state {s } N ,

• compute the energy change,

• accepte or reject the move proposal based on A({s } N |{s} N ).

Monte Carlo in the Gibbs ensemble: implementation and setup

In the GEMC method, we simulate two independent phases by two boxes simultaneously and allow them to exchange volume and particles. To realize this, three types of trial moves ( relax both systems to the pre-defined temperature with an acceptance ratio of

A j ({s } N j |{s} N j ) = min(1, exp -βF 0 (N j , L j , {s } N j ) exp -βF 0 (N j , L j , {s} N j ) ), (4.9)
where we need to keep the volume (V 1 , V 2 ) and number of particles (N 1 , N 2 ) in both boxes unchanged,

• a random volume rearrangement (V 1 , V 2 → V 1 , V 2 )
to equate the pressure of these two boxes with an acceptance ratio of

A(V 1 , V 2 |V 1 , V 2 ) = min(1, V 1 V 1 N 1 V 2 V 2 N 2 e -β∆(U 1 +U 2 ) ), (4.10)
where ∆U j (j = 1, 2) is

∆U j = F 0 (N j , L j , {s} N j ) -F 0 (N j , L j , {s} N j ), (4.11)
and

V 1 + V 2 = V 1 + V 2 = V
• a transfer of a randomly chosen particle between the two sub-systems (N 1 , N 2 → N 1 , N 2 ) to equate the chemical potential with an acceptance ratio of

A(N 1 , N 2 |N 1 , N 2 ) = V 2 V 1 N 1 N 2 e -β∆(U 1 +U 2 ) , (4.12)
where ∆U j (j = 1, 2) is

∆U j = F 0 (N j , L j , {s} N j ) -F 0 (N j , L j , {s} N j ), (4.13)
and

N 1 + N 2 = N 1 + N 2 = N .
The derivation of the acceptance ratio can be found in [START_REF] Panagiotopoulos | Direct determination of phase coexistence properties of fluids by monte carlo simulation in a new ensemble[END_REF]. Here I show an example of how to derive Eq. 4.12. The relative weight of two microstates in the Gibbs ensemble is:

ρ(N 1 , L 1 , {s} N 1 , N 2 , L 2 , {s} N 2 ) ρ(N 1 , L 1 , {s} N 1 , N 2 , L 2 , {s} N 2 ) = V 1 N 1 V 2 N 2 N 1 !N 2 !Λ 3(N 1 +N 2 ) dB exp -βF 0 (N 1 , L 1 , {s} N 1 exp -βF 0 (N 2 , L 2 , {s} N 2 V N 1 1 V N 2 2 N 1 !N 2 !Λ 3(N 1 +N 2 ) dB exp (-βF 0 (N 1 , L 1 , {s} N 1 ) exp (-βF 0 (N 2 , L 2 , {s} N 2 )
.

After plugging in

N 1 = N 1 , N 2 = N 2 , {s} N 1 = {s} N 1 , {s} N 2 = {s} N 2
, it can be reduced to:

ρ(N 1 , L 1 , {s} N 1 , N 2 , L 2 , {s} N 2 ) ρ(N 1 , L 1 , {s} N 1 , N 2 , L 2 , {s} N 2 ) = V 1 V 1 N 1 V 2 V 2 N 2 e -β∆(U 1 +U 2 ) (4.15)
where ∆U 1 and ∆U 2 are defined in Eq. 4.11.

The implementation of the Markov chain described in the previous section is relatively straightforward. We define three probabilities η Displ. = 0.5, η Vol. = 0.25, and η Part. = 0.25 for the choice among the different moves for each cycle, such as η Displ. + η Vol. + η Part. = 1. We define δs max as the maximal (normalized) particle displacement, which can be different for the two boxes, chosen such as to reach an acceptance ratio for moves of type (i) around 50% for each sub-system, and δV max as the maximal volume exchange, similarly optimized for an acceptance rate of about 50%. The detailed values are compiled in Tab. 4.2.

In particular for sodium we used as initial conditions a liquid box of 17 Å with 80 atoms and an empty vapor box of 25 Å for the simulations below 2000 K, and a liquid box of 18 Å with 80 atoms and an empty vapor box of 18 Å at 2000 K.

One accepted configuration at 2000 K served as starting point for the simulations at 2100K, 2200 K, 2300K and 2400 K. At 2500 K, we start with two 40-atoms boxes of 17 Å. In addition, we have performed an extra simulation at 2000 K with 120 atoms to estimate the finite-size effect, which is found to be small at this condition (see Figure 4.4). All the simulation boxes are cubic and the dimensions above define the edge of the box. The equilibrium values and their uncertainty were calculated using the autocorrelation technique [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF]. The error bars reported in the following sections are one-sigma error bars.

Setup of the DFT calculations

We compute the energy of the two sub-systems at each MC step using first-principles calculations in the projector augmented wave (PAW) method [START_REF] Blöchl | Projector augmented-wave method[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] of the DFT in the VASP (Kresse and Furthmüller, 1996;Kresse and Furthmüller, 1996) implementation. We employ the Generalized-Gradient Approximation in the Perdew-Burke-Ernzerhof formalism [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] for the exchange correlation term. We treat the 3s 1 as valence electron configurations for the PAW pseudopotentials. The partial occupancies for the electronic calculation are calculated using a Fermi-Dirac smearing scheme with a width corresponding to the nuclear temperature. The energy cut-off for the plane-wave basis set was set to 400 eV. The break condition for the electronic self-consistent loop was 10 -4 eV. The number of electronic bands was adapted to the temperature conditions such as to cover the entire spectrum of the fully and partially occupied states and to include enough non-occupied bands. The Brillouin zone was sampled with the Baldereschi point [START_REF] Baldereschi | Mean-value point in the brillouin zone[END_REF]. A test with a grid of 2 × 2 × 2 k-points yields results within 0.5% difference in energy. At each MC step, for each ionic configuration we compute the effect potential energy F 0 as defined in Eq.2.11 including the free energy of the electrons in a coulombic potential due to the ions and the ion-ion interaction. Our implementation is done outside of the VASP package (Kresse and Furthmüller, 1996), which is called only as an energy routine.

To compare with the Gibbs ensemble method, we also performed molecular dynamics simulations of the liquid phase in the DFT framework at several densities along the 2000 K isotherm. In order to stay consistent with the MC calculations we used the same DFT parameters for the MD simulations. The temperature was kept constant thanks to a Nosé thermostat [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF]. We used a fixed volume cell containing 80 atoms. The time step was set to 2 fs for a total duration of 20 ps.

Results and discussion

Stationary state and equilibrium

At the very beginning of the simulation, there is a net particle flux from the liquid box to the vapor box since the latter is initially empty (Figure 4.2). The driving force is the difference in chemical potentials. Because of the random character of the acceptance of moves along the Markov chains, particles from the gas box may also be transferred to the liquid box. This is captured by fluctuations of the density in each of the two boxes. After about 20000-50000 attempted moves, both boxes have reached a stationary state and the equilibrium is achieved. At equilibrium, for a temperature of 2000 K the cell length of both boxes fluctuates around 18 Å. The pressure in the liquid phase is 0.0 ± 0.2 kbar and 0.08 ± 0.03 kbar for the vapor phase. The liquid phase has 76 atoms on average but the vapor phase only 4. The stationary state thus seems to correspond to a thermal, dynamical and chemical equilibrium.

The acceptance ratios reported in Tab. 4.2 are very satisfactory with typical values between 10 and 75 % for each move. The 1200 K simulation shows a very low acceptance rate for the particle exchange because the temperature is very low compared to the energy barrier. Since the acceptance rate is non zero, the results are still reliable since a long enough simulation has been run. Figure 4.2 shows the evolution of different quantities as a function of the MC step at 2000 K. Table 4.1 lists all the values of the thermodynamic quantities.

Liquid-vapor equilibrium

As can be seen from Figure 4.3, at low temperature there is a clear distinction between the low-density and the high-density phase that corresponds to the liquid and vapour phase. It is then possible to determine the average thermodynamic quantities of both the vapor and the liquid phases by averaging over each distinct distribution. We perform an extra simulation at 2000 K with 120 atoms to estimate the finite-size effect on the liquid-vapour equilibrium density.

2

At 2400 and 2500 K, only one phase is present. The density, pressure, effective energy of gas and liquid phase indicate the average properties in each box.

Table 4.2: Acceptance ratio for the different moves in the present ab initio Gibbs ensemble study at 1200, 1500, 1800, 2000, 2100, 2200, 2300, 2400 and 2500 K. The maximum displacement δs max is in reduced coordinates of [0,1). et al., 1966) (black stars). The blue line is a fit of the equilibrium line using the scaling law and the red line is the law of rectilinear diameter with A = 0.80 ± 0.02 and B = 0.07 ± 0.01 parameters (see text for more details).

As temperature increases the two peaks become less and less separated, eventually preventing for a clear difference between the two phases. As we approach the critical temperature, the simulations have a higher probability of switching identity or even having two phases at once in the same simulation box. The latter is due to a comparable magnitude of the surface tension effect and the entropy contribution as already observed by [START_REF] Smit | Computer simulations in the gibbs ensemble[END_REF]. It results in the appearance of three peaks in the density distribution plot. In order to better quantify the density of the three phases (gas, liquid and the mixed phase), we fitted the three peaks by three Gaussian functions (solid black line in Figure 4.3) as suggested by [START_REF] Smit | Computer simulations in the gibbs ensemble[END_REF]. The center of the Gaussian is assumed to be the average density of each phase and its width is the standard deviation entering in the determination of the uncertainty. The low density peak corresponds to the gaseous phase and the high density peak to the liquid phase. The middle peak, close to the average density of the two boxes is the mixed phase, and is disregarded in the liquid-vapor equilibrium analysis. At 2300 K, the fluctuations become extremely large and we thus decided to only show the results in the density plot for reference but not to use them for the fit since they offer too loose a constraint on the critical point. For 2400 and 2500 K both boxes reach a very similar equilibrium and seem identical. This means that these conditions are above the critical point.

Based on the equilibrium densities, we can plot a vapor-liquid coexistence curve as shown in Figure 4.4. In general we obtain a good agreement compared to ex-perimental data available in the literature [START_REF] Dillon | Measurement of Densities and Estimation of Critical Properties of the Alkali Metals[END_REF]. At 1200 K, the saturated liquid density is 20% lower than that of experiments, which may be due to the choice of the exchange correlation functional that has already been noted in previous studies of argon [START_REF] Goel | Predicting vapor liquid equilibria using density functional theory: A case study of argon[END_REF]. A precise direct determination of the critical point is made difficult by the finite size of our system since the correlation length is expected to tend to infinity at the critical point; this cannot be captured within our small simulation cells. The critical point may however be approximated using the law of rectilinear diameter [START_REF] Rowlinson | Liquids and liquid Mixtures[END_REF]:

ρ L + ρ V = 2 ρ c + B 1 - T T c (4.16)
together with the scaling law [START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF])

ρ L -ρ V = A 1 - T T c β , (4.17)
where ρ L and ρ V are the densities of the coexisting liquid and vapor, at a given temperature T. The fitted parameters are the critical temperature T c and density ρ c , and the two constants A and B. β is the critical exponent, which is fixed here at 0.326 [START_REF] Wilding | Critical-point and coexistence-curve properties of the lennard-jones fluid: A finite-size scaling study[END_REF], as for other three dimensional systems [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF][START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF].

We obtained the critical point by applying the scaling law to all data points above 2000 K. Using our data we obtain our best fit with A = 0.87±0.1 and B = 0.19±0.03.

The critical point lies at 2338 ± 108 K and 0.24 ± 0.03 g/cm 3 . We stress here that the density values at 1800 K are compatible with the extrapolation of the scaling law, providing confidence in our fit. Our theoretical critical temperature is slightly lower than the experimental value at 2573 ± 171 K [START_REF] Dillon | Measurement of Densities and Estimation of Critical Properties of the Alkali Metals[END_REF], and the critical density is similar to the experimental value of 0.21 ± 0.02 g/cm 3 . We want to underline that both our values and the experimental ones are the result of extrapolations, as in both cases it is too challenging to obtain equilibrium data in the very vicinity of the critical point. With this in mind, the good agreement that we obtain between our calculations and experiments confirms the suitability of the ab initio Gibbs ensemble method for the determination of accurate coexistence curves. We also include the Clausius-Clapeyron plot of the saturated vapor pressure and density as a function of the inverse temperature in Figure 4.6. We obtain a nice affine behavior for both the logarithm of the density and of the pressure on these plots. We also have a relatively good agreement with the experimental density data.

In order to check our Gibbs ensemble results, we also performed a set of ab initio molecular dynamics simulations in the canonical ensemble at 2000 K for different densities. This allows us to analyse the structure of the liquid and to determine the corresponding spinodal point [START_REF] Speedy | Stability-limit conjecture. an interpretation of the properties of water[END_REF]. a clear minimum close to 0.48 g/cm 3 ; this is the liquid spinodal point. This is the smallest density at which the liquid is metastable. Above this density the fluid is homogeneous, as shown for example in the insets of Figure 4.5. At lower densities, in the unstable branch, bubbles form proving that the liquid becomes unstable. The density of the spinodal point is close yet lower than the equilibrium density of 0.50 ± 0.02 g/cm 3 predicted by the Gibbs ensemble method. We thus have a full consistency between these two completely different methods ensuring the reliability of the Gibbs ensemble method.

Structure of the liquid

We compare the structure of the liquid as we obtain it using the MC Gibbs ensemble and the MD approach. We analyse the radial distribution function (RDF)(also see Section 3.2.3) and stress here that the calculation of the RDF in the Gibbs ensemble is performed on a series of snapshots, and the number of particles and volume of each phase fluctuate. Figure 4.7 shows the RDF at several temperatures as extracted from our Gibbs ensemble simulations. The main peak lies around 3.5 Å. Along the vapor-liquid equilibrium line, the position of the peak changes slightly and broadens due mostly to the temperature effects. The spherical integration of the RDF from 0 to its first minimum gives the coordination number. At 2000 K and 0.52 g/cm 3 the agreement of the RDF as obtained in the MC and in the MD simulations is excellent. As the two methods start from different initial configurations and use different paths, they give a remarkable consistent outcome as they explore the configurational space. While this is not an absolute proof that we achieved ergodicity, it strongly suggests that our simulations are satisfying it. [START_REF] Dillon | Measurement of Densities and Estimation of Critical Properties of the Alkali Metals[END_REF]. The dashed line on the right graph denotes the pressures from fitting vapor density to the ideal gas law. As expected, a good agreement with calculated pressure is achieved in the low density range. In the high density region, a deviation is observed.

Conclusions

We have implemented the ab initio Gibbs ensemble algorithm and performed a series of simulations to compute the liquid-vapor equilibrium line and the critical point of sodium. We emphasize the electronic contribution at finite temperature, which is essential for metallic systems, but not always clearly explained in the literature. The effective nuclear potential energy defined in Eq. 2.11 should be used in the acceptance rule for a MC trial move in the Gibbs ensemble. We demonstrated that our simulations reached a mechanical and chemical equilibrium and the calculated phase coexistence curve and critical point of sodium are in a good agreement with the experimental results. The comparison of our results with molecular dynamics also showed very good consistency. Therefore, we confirm the reliability and validity of the ab initio Gibbs ensemble method.

Chapter 5

Ab initio Gibbs ensemble study of the liquid-vapor equilibrium and critical point of iron

Introduction

Iron, as the major constitutive material of the metallic core in a telluric planet and as a basic engineering material, has received and continues to receive considerable attention from academics across multiple disciplines. Significant efforts were made for studying the thermodynamic and thermophysical properties of its liquid and solid phases for a wide range of temperature and pressure conditions by either experimental or theoretical methods. These information have been used, on the one hand, to design new technological devices, and on the other hand, to infer the physical state and investigate the thermodynamic evolution of planetary cores. Recent developments in the field of extreme shock experiments and giant impact simulations led to an increased interest in obtaining a complete thermal equation of state (EOS) that can describe all vapor, liquid and solid phases of iron (Canup, 2004a). However, only a few studies have been performed to examine the vapor-liquid equilibrium (VLE) and the critical point (CP) of iron [START_REF] Assael | Reference Data for the Density and Viscosity of Liquid Aluminum and Liquid Iron[END_REF][START_REF] Grosse | The densities of liquid iron and nickel and an estimate of their critical temperature[END_REF][START_REF] Hixson | Sound speed and thermophysical properties of liquid iron and nickel[END_REF]. The scarceness of data inhibited building an accurate set of EOS for iron that spans a wide-enough density and temperature region to cover all the different phases of interest. These limitations cause problems in predicting the chemical and dynamic consequences in a giant impact simulation. Therefore, it is imperative to close the knowledge gap of the VLE and CP of iron.

However, it is still challenging to determine the VLE and CP experimentally above 4000 K. Computer simulations might provide an alternative option and sometimes are the only viable way. The Gibbs ensemble Monte Carlo (GEMC) method, coupled to finite-temperature density functional theory, has just opened up the possibility of studying the VLE and CP for a metallic system with a reasonable computational cost. Its first application to sodium has shown an excellent agreement with available experimental results (see Chapter 4). Unfortunately, we cannot directly apply this technique to iron due to the additional complication from the high-temperature magnetism.

Both experiments [START_REF] Waseda | Atomic distribution and magnetic moment in liquid iron by neutron diffraction[END_REF] and theoretical simulations [START_REF] Lichtenstein | Finite-temperature magnetism of transition metals: An ab initio dynamical mean-field theory[END_REF] suggest that the liquid iron might be in a paramagnetic state with a non-vanishing local magnetic moment at ambient pressure. Therefore, we need to take into account the magnetic excitation or the magnetic degrees of freedom, including the transversal fluctuation, which changes the orientation of the local magnetic moment, and the longitudinal fluctuation that modifies the magnitude of the local magnetic moment [START_REF] Kaul | Phase transitions and finite temperature magnetism: Experiment and analysis[END_REF]. We must remind the readers that a spin-polarized DFT calculation with a fully self-consistent electronic optimization to find the electronic ground state will only yield the size and direction of the local magnetic moment that corresponds to the zero temperature on the magnetic degrees of freedom [START_REF] Gambino | Longitudinal spin fluctuations in bcc and liquid fe at high temperature and pressure calculated with a supercell approach[END_REF]. As pointed out by [START_REF] Abrikosov | Recent progress in simulations of the paramagnetic state of magnetic materials[END_REF], the variation in the size and orientation of the local magnetic moment observed in the spin-polarized ab initio molecular dynamics simulations is only a direct impact of atom vibrations.

Several new approaches have been developed to study the magnetic excitation. Among them, the accurate one is based on the dynamic mean-field theory [START_REF] Lichtenstein | Finite-temperature magnetism of transition metals: An ab initio dynamical mean-field theory[END_REF]. Because of the computational limitations and poor tractability, this technique is still limited to a small system size that is insufficient to properly describe the structural disorder of a fluid, even metallic like iron. To reduce the complexity to a manageable extent, we first note that each iron atom can be associated with a well-defined local magnetic moment that behaves in a Heisenberg-like manner, although 3d itinerant electrons cause the magnetism [START_REF] Abrikosov | Recent progress in simulations of the paramagnetic state of magnetic materials[END_REF]. It suggests that the longitudinal fluctuation might play a limited role and is thus often neglected. We can further simplify the problem by considering that the magnetic fluctuations adiabatically follow the atomic motions. This treatment is a special case of the disordered local moment molecular dynamics (DLM+MD), where the paramagnetic state is realized by the statistical average of many different random magnetic configurations known as the magnetic sampling method [START_REF] Alling | Effect of magnetic disorder and strong electron correlations on the thermodynamics of crn[END_REF]. In this way, we treat the atomic vibrations and the transversal fluctuations on a similar footing. However, the incorporation of the magnetic degrees of freedom into the aiGEMC method is too expensive. Fortunately, as we shall see later, the paramagnetic state can be well approximated by the ferromagnetic state in terms of the liquid-vapor equilibrium densities.

In the present study, we employed the aiGEMC method to determine the liquidvapor equilibrium line and the critical point of iron in the non-magnetic (NM) and ferromagnetic (FM) state. Then we used DLM+MD method to determine the equilibrium density for the liquid iron from 1800-7500 K at 0 GPa, where the paramagnetic (PM) state is treated by the magnetic sampling method.

Simulation details

In the aiGEMC method, we simulate two independent phases (liquid and vapor) by two boxes simultaneously but allow them to exchange volume and particles to achieve the chemical equilibrium. We employed finite-temperature density functional theory [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF][START_REF] Mermin | Thermal properties of the inhomogeneous electron gas[END_REF] as implemented in the Vienna Ab initio Simulation Package (VASP) (Kresse and Furthmüller, 1996;Kresse and Furthmüller, 1996) to compute the effective potential energy for the Monte Carlo acceptance criteria. The Perdew-Burke-Ernzerhof functional [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] and the projector augmented wave (PAW) potential [START_REF] Blöchl | Projector augmented-wave method[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF]) with 3d 7 4s 1 as valence electrons were used.

The partial occupancies for the electronic structure calculation were calculated using a Fermi-Dirac smearing scheme with a width corresponding to the simulation temperature. The energy cut-off for the plane-wave basis set was set to 400 eV and the Brillouin zone was sampled with the Baldereschi point [START_REF] Baldereschi | Mean-value point in the brillouin zone[END_REF]. For the ab initio molecular dynamics simulations, we used the same DFT parameters as the MC calculations to be consistent. The temperature was controlled by the Nosé thermostat [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF]. The timestep of the aiMD simulation was set to 1 fs for a total simulation time of at least 3 ps. For the paramagnetic phase, the transversal fluctuation is taken into account by the magnetic sampling method. We exploited 40 magnetic configurations at every MD step to obtain the average force exerted on atoms, where the initial local magnetic moments for the electronic minimization were set to be orientated in a randomly spin-up or down direction but with a zero total magnetic moment in the simulation cell [START_REF] Alling | Effect of magnetic disorder and strong electron correlations on the thermodynamics of crn[END_REF]. For the ferromagnetic phase, we reinitialized the system to be ferromagnetic at every MD step in order to prevent the electronic subsystem from getting stuck in a local minimum during the electronic minimization. In this study, we used 80 atoms for both aiGEMC and aiMD simulations. But an extra aiGEMC simulation with 108 atoms at 7500 K has been performed to estimate the finite size effects. First, we approximate the paramagnetic state by a nonmagnetic simulation. We plot the density evolution as a function of MC step in Fig. 5.1(a). In order to verify that the simulation results are independent of the initial condition, including the starting volume and atomic positions for each box and the total volume of the two boxes, we perform two aiGEMC simulations at 3500 K. The first simulation employs two identical boxes, each with 40 particles and the same atomic position. The second simulation starts with one empty box and one box with 80 atoms. However, these two simulations have different total volumes.

Results and Discussion

In the former case, we observe an instantaneous phase separation, i.e. all the atoms move into one box, leaving the other one empty (Fig. 5.1(a)). When the atoms are transferred into one fully populated box, its volume also changes to accommodate this amount of atoms. Even if one box becomes empty, there is still a volume exchange between the two boxes. Eventually, the simulation converges to 10.11 3 Å 3 for the liquid box, which corresponds to 7.18 ± 0.03 g/cm 3 density. In the second case, we did not observe any particle transfer between the two boxes, but the volume exchange is still present. Its final volume is 10.15 3 Å 3 , which corresponds to 7.10 g/cm 3 ± 0.05 density. Consequently, we observe that the two simulations yield the same density for the liquid iron phase, suggesting that the value of the density we obtain using this procedure is reliable.

As we approach the critical point, we observe that the density of the equilibrium liquid drops with increasing temperature. We also observe the appearance of the first atoms in the vapor box at 5000 K, and the equilibrium vapor density increases with increasing temperature. At 7500 K, we perform another simulation with 108 iron atoms to estimate the finite-size effects, which is found to be small.

We plot the unnormalized density probability distribution functions in Fig. 5.

1(b).

There are two well-separated peaks in the density histogram below 7500 K, which corresponds to the liquid and vapor phase. At 8500 K, there is a significant fluctuation in the density evolution since the temperature is very close to the critical temperature. At 9500 K, the two boxes frequently exchange their identities, and a single peak is present, indicating the simulated temperature is above the critical temperature.

Coming back to the simulations at 3500 K, the extrapolated density of the vapor phase at this temperature should be on the order of 10 -4 g/cm 3 . It implies that we need to exploit a simulation cell of at least 200 3 Å 3 . Due to the limitation on the computation facilities, we could not further increase the simulation cell beyond the 70 3 Å 3 which is used in the present study at 3500 K. However, an empty vapor box should not affect the accuracy of the liquid equilibrium density and position of the critical point. The empty box without any atom acts like a pressure buffer and sets the pressure on the liquid phase to zero. Although the real equilibrium pressure at 3500 K is not zero if a sufficiently large simulation cell is employed, it should be very close to zero. As we will explain in Section 5.3.4, a variation of even 1 GPa in pressure only changes the density of the liquid by less than 1% at 3500 K. Therefore, the zero pressure density for the liquid phase is an excellent approximation to the equilibrium liquid density. In a second step, we take into consideration the magnetism and perform aiGEMC simulations using the ferromagnetic state for iron. We plot the density evolution as a function of MC step and corresponding unnormalized probability distribution functions in Fig. 5.2. We recognize a very similar trend in the FM-aiGEMC simulations compared to the NM-aiGEMC simulations. At 3500 K, we start with one empty box and a full box with 80 atoms. We do not observe any particle transfers between the two boxes, but the volume exchange is still present. The simulation yields the liquid density of 6.47 ± 0.04 g/cm 3 . With increasing temperature, the equilibrium liquid density falls. At 5000 K we start to observe a few atoms in the vapor box, and the equilibrium vapor density increases with increasing temperature.

From the density histogram we can see two well-separated peaks in the density histogram below 7500 K and a significant fluctuation in the density evolution above 7500 K as the temperature is very close to the critical temperature. At 8000 K, we identify three peaks in the density histogram for the FM-aiGEMC simulation, where the middle peak is due to the comparable magnitude of the surface tension effect and the entropy contribution as already mentioned by [START_REF] Smit | Computer simulations in the gibbs ensemble[END_REF]. At 9500 K, the two boxes frequently exchange their identities, and we only observe a single peak, suggesting the simulated temperature is above the critical temperature.

The critical point

The liquid-vapor equilibrium line for iron obtained from the aiGEMC simulations is plotted in Fig. 5.3. As the correlation length tends to be infinity near the critical point, a direct determination of the critical point from finite simulation cells is difficult. However, we might approximate it by applying the scaling law (Eq. 4.17) and the rectilinear law (Eq. 4.16) to the equilibrium liquid and vapor densities above 7500 K. This method has been successfully applied to sodium and the resulted critical point compare well with experimental results (see Chapter 4). For the ferromagnetic case, we obtain the critical point at 1.64 ± 0.23 g/cm 3 and 8690 ± 141 K. For the nonmagnetic case, the critical point lies at 1.41 ± 0.17 g/cm 3 and 8701 ± 100 K. The nearly identical critical point for iron in the nonmagnetic and ferromagnetic state indicates the effect of magnetic state of iron is limited. We extrapolate the saturated vapor pressure to the critical density and obtain the critical pressure of 2.0 kbar and 2.4 kbar for iron in the nonmagnetic and ferromagnetic state, respectively. The critical point obtained in the present study agrees very well with [START_REF] Medvedev | Wide-range multiphase equation of state for iron[END_REF], which lies at 1.638 g/cm 3 and 8763 K. In their study, a semi-empirical equation of state for iron with several adjustable parameters was developed and these parameters were determined by fitting to available experimental data. Considering the scarceness of experimental data for iron in the low-density regime, the excellent agreement is unexpected.

The bulk modulus

If neglecting the difference between the NPT ensemble and the Gibbs ensemble with a continuous change in the number of atoms for each box, we can estimate the bulk (a) The liquid-vapor equilibrium line for iron obtained from the aiGEMC simulations. The blue and green color denote simulations for iron in the nonmagnetic and ferromagnetic state, respectively. The solid circles represent the resulted liquid or vapor densities from the aiGEMC simulations, while empty triangles are densities determined by the law of rectilinear diameter. The blue and green line is a fit of the equilibrium liquid and vapor densities to the scaling law and the law of rectilinear diameter. The predicted critical point is represented by the solid stars. (b) The Clausius-Clapeyron plot of the logarithm of the saturated vapor pressure obtained from the aiGEMC simulations as a function of temperature. [START_REF] Hixson | Sound speed and thermophysical properties of liquid iron and nickel[END_REF][START_REF] Grosse | The densities of liquid iron and nickel and an estimate of their critical temperature[END_REF][START_REF] Beutl | Thermophysical properties of liquid iron[END_REF] aiGEMC-NM aiGEMC-FM Figure 5.5: The zero pressure density up to 7500 K obtained from ab initio MD simulations. For comparison, we also plot the equilibrium liquid iron density derived from aiGEMC simulations and experimentally determined liquid iron density at around 0 GPa [START_REF] Beutl | Thermophysical properties of liquid iron[END_REF][START_REF] Grosse | The densities of liquid iron and nickel and an estimate of their critical temperature[END_REF][START_REF] Hixson | Sound speed and thermophysical properties of liquid iron and nickel[END_REF]. modulus for the liquid phase based on the pressure-density fluctuation (Fig. 5.4). The bulk modulus of the liquid iron phase in the NM state decreases from 250 GPa at 3500 K to 10 GPa at 7500 K, while in the FM state, it falls from 95 GPa at 7500K to 8 GPa at 7500 K. The large discrepancy at 3500 K is mainly due to the difference in the liquid iron density, which is 7.18 g/cm 3 in the NM state compared to 6.47 g/cm 3 in the FM state. The large bulk modulus for the liquid iron suggests a variation of 1 GPa in the pressure only changes the density by at most 1% at 3500 K and 10% at 7500 K. Therefore, we can safely use the density at 0 GPa to approximate the density at the liquid-vapor equilibrium density since the difference is insignificant. We could not give a reliable estimation of the bulk modulus at 8500 K due to a substantial density fluctuation. However, our previous MD simulations at 8100 K and 8750 K (see Fig. 3.1) suggest the bulk modulus should be very small since a pressure deviation of 2 kbar changes the density by 1 g/cm 3 .

The liquid density at zero pressure from ab initio MD simulations

In order to estimate the density of liquid iron at zero pressure, we perform ab initio MD simulations in the NVT ensemble to extract the average pressure at two or three densities, from which the zero pressure density is obtained by a linear interpolation.

For the liquid iron in the FM state, we notice that the electronic optimizations in the ab initio MD simulation are often trapped in a local minimum. We show an example of such simulations starting with two different initial configurations at 3500 K and 6.63 g/cm 3 in Fig. 5.6. In both simulations, the liquid iron was initialized to be ferromagnetic, and the wavefunction and charge density extrapolation scheme as implemented in VASP was employed. From Fig. 5.6(a) and (d), we can see that there is a large energy drift, and the total magnetization of the system decreases steadily to zero and then oscillates around it. To reveal the origin of the energy drift, we select several snapshots and re-calculate the effective potential energy as defined in Eq. 2.13 by re-initializing the system to be ferromagnetic. The resulted potential energy is found to be lower than the one using the extrapolation scheme, and both the effect potential energy and the total magnetization fluctuates around a constant value.

The energy difference might be caused by the different initial trial wavefunction employed as a starting point for the electronic optimization. As there are many local minima in the potential energy surface for a spin-polarized calculation, the electronic minimization algorithm may be trapped in such a minimum that is related to the initial conditions. Since the electronic minimization fixes the magnetic degrees of freedom at 0 K, the only physically meaningful magnetic configuration is the one with the lowermost effective potential energy. In this case, it is the ferromagnetic state. Therefore, in the present study we re-initialize the system to be ferromagnetic at every MD step in order to maintain the system to be ferromagnetic. We also stress that the simulations after 1 ps shown in Fig. 5.6(c) and (f) do not represent a paramagnetic state due to the lack of spin dynamics, even if the total magnetization is around zero.

We employ the DLM + MD method to evaluate the effect of paramagnetism on the liquid density at zero pressure. In order to determine the number of magnetic configurations needed to converge the force exerted on the atoms, we choose two snapshots from FM-aiMD simulations at 1850 K and 7.22 g/cm 3 and generate 1000 magnetic configurations by randomly assigning the initial magnetic moments in the up and down direction but with a zero total magnetic moment. In the DLM+MD scheme, we allow the electronic minimization and the spin dynamics is introduced by the average over many magnetic configurations. Fig. 5.8 and Fig. 5.9 show the force exerted on one randomly selected atom, pressure, the local magnetic moment of one randomly selected iron atom, the total magnetization, and their moving averages over magnetic configurations. We find that the force and pressure converged to 5% or better after employing more than 40 magnetic configurations. Therefore, we choose 40 magnetic configurations for the PM-aiMD simulations. It is also interesting to see that the selected atom in the first snapshot has a local magnetic moment of about 2.8 µ B which is independent of the initial magnetic configuration, while in the second snapshot, it is more spread in a range of 0-2 µ B . The difference may be caused by the variation in the local atomic environment.

The zero pressure densities of the liquid iron from aiMD simulations are plotted in Fig. 5.5. For iron in the nonmagnetic and ferromagnetic state, the determined densities are in a good agreement with the values determined by the aiGEMC simulations with a difference that is less than 10%. In addition, the zero pressure density for iron in the paramagnetic state are very close to the values for iron in the ferromagnetic state, but is larger than the densities for iron in the nonmagnetic state.

In order to understand this trend, we show the distribution of the local magnetic moment for iron atoms from FM-aiGEMC, FM-aiMD, and PM-aiMD simulations in Fig. 5.10. As expected, the magnetic moment distribution from FM-aiGEMC matches well with FM-aiMD simulations in all temperature conditions considered here. With increasing temperature, a small number of atoms flip their magnetic moment direction to be negative, and more particles become nonmagnetic. For the PM-aiMD simulation at 1850 K, there are two distinct peaks centered at 2.3 µ B and -2.3 µ B . Starting from 5000 K, the third peak at around 0 µ B appears and its height increases with temperature, suggesting more atoms become nonmagnetic as well. The average size of the magnetic moment as a function of temperature is displayed in Fig. 5.11. It falls with increasing temperature due to the competing effect of the density, electronic Fermi smearing, and thermal motion. The decreasing density would cause the strong localization of the 3d electrons and boost the local magnetic moment, while the latter two factors reduce it. The average magnetic moment in FM is larger than PM, resulting in a larger effective atomic volume and thus a slightly lower zero pressure density. The decreasing average magnetic moment with temperature also reduces the difference of the zero pressure density between the nonmagnetic state and the magnetic state. Therefore, it is reasonable to expect that the magnetic state only slightly affect the critical point.

Comparing to the liquid iron densities determined by experiments at around zero pressure [START_REF] Beutl | Thermophysical properties of liquid iron[END_REF][START_REF] Grosse | The densities of liquid iron and nickel and an estimate of their critical temperature[END_REF][START_REF] Hixson | Sound speed and thermophysical properties of liquid iron and nickel[END_REF], our NM-aiMD and NM-aiGEMC simulations considerably overestimate the liquid density by 1.2 g/cm 3 (18%) at 1850 K, while the aiGEMC and aiMD calculations with iron in the FM or PM state overestimate the liquid density by 0.4 g/cm 3 (5%) at 1850 K. The improved agreement given by spin-polarized simulations compared to non-spin-polarized simulations confirms that the magnetism plays an important role in determining the physical properties of liquid iron at high temperature. At 3500 K, the liquid iron density in the ferromagnetic or paramagnetic state from our simulations is 0.4 g/cm 3 denser than the experimental value, although the experimental data are scattered above 3000 K. The discrepancy might be caused by the PBE exchange-correlation functional. Previous DLM-MD simulations for bcc and fcc iron at experimental densities from 300 K to 1662 K reported a negative pressure of 7 GPa, indicating PBE functional overestimate the zero pressure density [START_REF] Alling | Strong impact of lattice vibrations on electronic and magnetic properties of paramagnetic fe revealed by disordered local moments molecular dynamics[END_REF]. More experiments, especially at above 3500 K, are also needed in order to have a better insight into the origin of the discrepancy.

The structure of the fluid

We analyze the fluid structure by the radial distribution function (RDF) extracted from our aiGEMC and aiMD simulations. For the liquid phase, the main peak lies at 2.4 Å and almost keeps constant from 3500 K to 8500 K. For the vapor phase, we only show two cases at 7500, 8000 and 8500 K, where the vapor simulation box contains more than ten atoms that allow the computation of EDF. As expected, we fail to identify any local structure in the vapor phase. Our results suggest that the liquid iron in the NM state has a very similar structure to that in the FM and PM states. At 3500 K, the agreement of the RDF as obtained by the aiGEMC and the aiMD simulations is excellent. As the two methods start from different initial configurations and sample the configurational space with different paths, the remarkable consistency strongly suggests that our simulations achieved ergodicity, although this is not a formal proof. The RDF's spherical integration from 0 to its first minimum gives the coordination number, which drops from 13 at 3500 K to 7 at 8500 K for the liquid phase.

Conclusions

We have combined the aiGEMC method and the aiMD method to determine the liquid-vapor equilibrium and the critical point of iron in the nonmagnetic, ferromagnetic and paramagnetic states. The paramagnetic phase is treated by the magnetic sampling method. By comparing the liquid density at zero pressure from 1850 K to 7500 K, our study shows the FM liquid iron phase is a good proxy to the PM liquid iron phase in terms of liquid-vapor equilibrium densities. We also found the difference of the zero pressure density between the nonmagnetic state and the magnetic state drops with temperature due to the decreasing average magnetic moment. We predict the critical point lies at 1.69 g/cm 3 and 8814 K for the magnetic phase, and at 1.38 g/cm 3 and 8697 K for the nonmagnetic phase. The improved agreement with available experimental results given by spin-polarized simulations compared to non-spin-polarized simulations suggests we need to consider the magnetism in the simulations. d)) as defined in Eq. 2.13, the effective potential energy (solid squares in (b) and (e)) as defined in Eq. 2.11 and the total magnetization (solid lines in (c) and (f)) of the system as function of time from the aiMD simulations starting with two different initial configurations at 3500 K at 6.63 g/cm 3 . In both simulations, the system was initialized to be ferromagnetic at the first MD step and the wavefunction and charge density extrapolation scheme as implemented in VASP was employed for the following MD steps. We selected several snapshots and performed the electronic minimization by re-initializing the system to be ferromagnetic. The resulted effective potential energy and the total magnetization (empty circles in (b) and (c)) is lower and higher than the one calculated with the wavefunction extrapolation scheme, respectively. It suggests that for the spin-polarized simulations, the results of the electronic minimization depends on the initial trial wavefunction or charge density. The constant of motion (a) and the total magnetization (b) of the system as function of time for FM-aiMD simulations at 3500 K and 6.50 g/cm 3 represented by the blue line, and 6.41 g/cm 3 represented by the green line, respectively. In both simulations, the system was initialized to be ferromagnetic at every time step. .12: Radial distribution function (g(r)) for the liquid phase at 3500, 5000, 7500 and 8500 K (a), and for the vapor phase at 7500 K and 8500 K (b), as computed with the aiGEMC simulations. The shaded areas correspond to our estimate of the one-sigma uncertainty. The curves were shifted for readability. (c) Comparing the g(r) resulted from aiGEMC simulations to that from aiMD simulations for iron in the different magnetic states at 3500 K. For aiMD simulations, we show two radial distribution functions that are at a density higher or lower than the zero pressure density. (d) The coordination number for the liquid iron as a function of temperature.

Chapter 6

Conclusions and future work

Conclusions

The prevailing theory to explain the origin of the Moon is the giant impact hypothesis, where a Mars-sized impactor collides with the proto-Earth, and the Moon is subsequently accreted from the proto-lunar disk made of the ejected materials. As laboratory-scale experiments are not able to simulate such planetary-scale impacts, our understanding of the giant impact mostly comes from hydrodynamic simulations. However, the outcome of these simulations depends heavily on the available equations of state to describe the thermodynamic response of the constitutive materials of the proto-Earth and of the impactor to extreme shock waves.

Iron as a building block material of the terrestrial planets naturally received significant attention. But the major effort has been put to determine its phase diagram up to the Earth's core conditions (126-360 GPa and 3000-7000 K) and beyond. The studies of iron at low densities and high temperatures are still scarce, causing great uncertainty in estimating its liquid-vapor equilibrium densities and the critical point. Consequently, it prevents us from developing an accurate equation of state for iron. As an essential ingredient in the hydrodynamic simulations, the equation of state would affect the impact energy distribution and thus cause different geochemical and geodynamic consequences, such as the amount of iron vaporized, the extent of iron-silicates equilibration, and the depth of the magma ocean. This study focuses on the thermodynamic and thermophysical properties of iron in the low-density regime to better understand iron's behavior during giant impacts.

In the first attempt, we employed ab initio molecular dynamics to study the mechanically stable limit of the liquid iron under the hydrostatic tension (i.e., the spinodal line), which helps locate the position of the critical point. We found the critical point lies in the temperature range of 9000-9350 K and the density range of 1.85-2.40 g/cm 3 corresponding to a pressures range of 4-7 kbars. We also computed two Hugoniot lines starting from two realistic initial states (1GPa and 1500 K, 40 GPa and 4000 K) to model the behavior of the iron core in planetesimals and planets during the impacts. By comparing the entropy values, determined by the two-phase thermodynamic method, along these Hugoniot lines to that at the boiling point of 1 bar and 3150 K, we found that the shock pressure required to reach the onset vaporization upon release and cooling is 312 GPa and 365 GPa, respectively. The application of the impedance matching method yields the peak condition achieved for a planetesimal colliding with the Earth. We found at least 70 % of the impacts in the late veneer results in a peak pressure larger than the threshold, suggesting the partial core vaporization is readily achieved. This would help mix the highly siderophile elements into magma ponds or oceans and explain the excess amount of siderophile elements in the Earth's mantle [START_REF] Rubie | Formation of the Earth's Core[END_REF]. For giant impacts, the geometry effect plays an essential role in controlling the shock peak conditions, and the impedance matching method is no longer valid. However, considering the significantly lower threshold for vaporization than the previous estimation (887 GPa) which was developed based on an inaccurate equation of state, our study suggests that Theia's core underwent partial vaporization and previous hydrodynamic simulation underestimate the vapor production during giant impacts. After cooling down, the vaporized iron would condensate into a distribution of small droplets on the centimeter-level. It would enhance the iron-silicates equilibration and easily explain the recent W-isotope data that requires at least 30% core-mantle equilibration in the aftermath of the giant impact.

Additionally, we characterize the structural and transport properties of the fluid iron in the low-density region, including diffusion coefficients, viscosity, and thermal conductivity. They may provide a better constraint on iron's behavior in the protolunar disk. The diffusivity increases by two orders of magnitude from about 1.0 × 10 -8 m 2 s -1 at 7.75 g/cm 3 and 3000 K to about 1.3 × 10 -6 m 2 s -1 at 0.37 g/cm 3 and 12000 K. We found the low-density fluid is highly depolymerized and is mainly made of isolated atoms. The coordination number increases strongly with density below 7.75 g/cm 3 , then saturates at a value of 14 that persists at least up to 6000 K and 13.3 g/cm 3 , which is the condition at the Earth's inner-outer core boundary. The viscosity of iron in the outer parts of the proto-lunar disk is also extremely low, on the order of 10 -3 Pa s. The computed electronic and thermal conductivities decrease with the density, suggesting a gradual reduction of the metallic character.

To better constrain the liquid-vapor equilibrium and the critical point of iron, we implemented the Gibbs ensemble Monte Carlo method, which avoids the interface issues and allows for a more reliable determination of the phase equilibrium. The Gibbs ensemble method, coupled with the density functional theory, has successfully been applied to the insulating system, where the internal energy of the electron-ion system is used in the Monte Carlo acceptance criteria. However, we need to consider excited electronic states for the metallic system at high temperatures. Therefore, it is necessary to utilize the finite-temperature density functional theory to capture these excited states. By taking a step back to statistical physics, we found the effective potential energy that includes the electronic free energy should be used in the Monte Carlo acceptance criteria. Then we applied this method to sodium where several experimental data are available. The obtained liquid-vapor equilibrium and critical point are in good agreement with the experimental results, confirming our implementation's reliability and validity.

However, the application of this method to iron is not straightforward due to the extra complication from the paramagnetism of the fluid iron at high temperatures. We first performed ab initio Gibbs ensemble Monte Carlo simulations for iron in the nonmagnetic and ferromagnetic states. Based on the pressure-density fluctuation, we found the bulk modulus of the liquid iron from 3500 K to 7500 K is large enough that a variation of 1 GPa in pressure only changes the density by a few percent. Therefore, we can approximate the liquid-vapor equilibrium densities for the liquid branch with the densities at zero pressure in this temperature range. This is confirmed by our ab initio molecular dynamics simulation results. In order to evaluate the effect of the possible paramagnetism, we employed the disordered local moment molecular dynamics to investigate the liquid density at zero pressure, where the paramagnetic state is treated by the statistical average of many random magnetic configurations known as the magnetic sampling method. Our study shows the ferromagnetic liquid iron phase is a good proxy to the paramagnetic liquid iron phase in terms of liquid-vapor equilibrium densities. The improved agreement with available experimental results for the liquid iron density at around ambient pressure up to 4000 K given by spin-polarized simulations compared to non-spin-polarized simulations suggests we need to consider the magnetism in our simulations. Our preliminary results predict the critical point lies at 1.69 g/cm 3 and 8814 K for the magnetic iron phase, and at 1.38 g/cm 3 and 8697 K for the nonmagnetic iron phase. But we are running more simulations from 7500 K to 8500 K in order to accurately determine the critical point.

Future work

Within the time available, several problems could not be resolved. However, they are necessary to better understand the thermodynamic properties of iron in the low-density regime. It includes the following:

1. The role of exchange-correlation functional on the liquid-vapor equilibrium of sodium and iron. In this study, we employed the PBE-GGA exchangecorrelation functional. It is a general trend that the GGA functionals underbind atoms, leading to a slighly longer bond length. Therefore, underestimating the liquid sodium densities around zero pressure, as we observed in Chapter 4, meets our expectations. For bcc iron at zero temperature, PBE underestimates its density as well. However, we found PBE overestimates the zero pressure densities of the liquid iron at high temperatures, even including the magnetism. Therefore, future work is needed to clarify the reason for the inconsistent performance of PBE on different iron phases. The inclusion of dispersion interactions may change the phase equilibrium, which needs to be investigated in the future. However, we expect that the vapor may suffer a more significant impact since its bulk modulus is very small, and a tiny variation of pressure would shift the density considerably.

2. The effect of the longitudinal fluctuation of the magnetization on the liquidvapor equilibrium of iron. The density difference of 0.5 g/cm 3 at 3500 K between experimental results and disorder local magnetic moment molecular dynamics simulations might be partially caused by the longitudinal fluctuation that is not considered in our study. It would shift the zero pressure densities to a lower value and improve the agreement with the experimental results. However, there is no consistent and cost-effective way to incorporate this effect in DFT calculations. Therefore, more works need to be done on improving the theoretical treatment for paramagnetism at high temperature.

Figure 1

 1 Figure 1.1: The dynamic evolution of the Moon. (a) The lunar orbital elements at present. The inclination of the lunar orbit to the ecliptic plane is 5 • . The mutual inclination of the Moon's orbit to the Earth's equator varies from 18 • to 28 • . (b)The evolution of the mutual inclination as a function of radius. The mutual inclination must be at least 12 • when the Moon forms. However, the giant impact theory predicts the moon formed at the Earth's equator predicting an initial value of 0 • . This discrepancy can be reconciled with the Moon-Sun resonance[START_REF] Touma | Resonances in the early evolution of the earth-moon system[END_REF], the disk-satellite resonance[START_REF] Ward | Origin of the moon's orbital inclination from resonant disk interactions[END_REF] or the Laplace plane transition(Ćuk et al., 2016), which can excite the mutual inclination to a high value (see the black dashed line in b). This figure is reproduced with permission from[START_REF] Touma | Evolution of the earth-moon system[END_REF]. (c)[START_REF] Ćuk | Making the moon from a fast-spinning earth: a giant impact followed by resonant despinning[END_REF] have shown the Earth's rotational rate slows down due to the Moon-Sun resonance which can be used to drain away the excess angular momentum. They suggest an initial high angular momentum at the time of the giant impact is dynamically feasible. However, subsequent studies using different tidal models found too much[START_REF] Wisdom | Early evolution of the earth-moon system with a fastspinning earth[END_REF] or too little[START_REF] Ward | Analytical model for the tidal evolution of the evection resonance and the timing of resonance escape[END_REF] angular momentum be removed from the Earth-Moon system, raising the question whether the Moon-Sun resonance is still a viable mechanism. This figure is reproduced with permission from[START_REF] Ćuk | Making the moon from a fast-spinning earth: a giant impact followed by resonant despinning[END_REF].

  Photosphere Convection

Figure 1

 1 Figure 1.3: (a) The phase diagram of the condensed iron up to 360 GPa and 8000 K. The melting curve was determined by the resistance-heated diamond-anvil cell technique[START_REF] Sinmyo | Melting curve of iron to 290 gpa determined in a resistance-heated diamond-anvil cell[END_REF]. The pressure at the core-mantle boundary (CMB) is around 130 GPa. The pressure at the inner-core boundary (ICB) is about 330 GPa. The empty red square is the melting temperature at ICB determined by the ab initio method[START_REF] Alfe | Temperature of the inner-core boundary of the earth: Melting of iron at high pressure from first-principles coexistence simulations[END_REF]. (b) The liquid-vapor equilibrium line in the temperature-density space. The black solid line and red dotted line represents the liquid-vapor equilibrium line estimated by[START_REF] Fortov | Shock waves and equations of state of matter[END_REF] and[START_REF] Grosse | The densities of liquid iron and nickel and an estimate of their critical temperature[END_REF], repectively. The blue point denotes the density of liquid iron at the boiling point of 1 bar. The magenta shared region is the experimental conditions done by[START_REF] Hixson | Sound speed and thermophysical properties of liquid iron and nickel[END_REF]. (c) The liquid-vapor equilibrium line in the temperature-pressure space. The supercritical state is represented by the shaded region.

Figure 3 . 2 :

 32 Figure 3.2: The static structure of hot fluid iron. (a) Pair distribution function for iron at selected temperatures and densities. The color of solid lines represents different temperatures, as shown in the legend of (b). The unit of those numbers is g/cm 3 . (b)The coordination number calculated by counting the number of atoms within the first minimum of g(r) from a central atom. There is a weak dependence of the coordination number with the temperature. However, its density dependence is very strong below 7.75 g/cm 3 , then saturate at a value of 13 up to 6000 K and 13.3 g/cm 3[START_REF] Alfè | Structure and dynamics of liquid iron under Earth's core conditions[END_REF][START_REF] Vočadlo | First principles calculations on crystalline and liquid iron at Earth's core conditions[END_REF], which is the condition at Earth's inner-outer core boundary.

Figure 3

 3 Figure3.3: Speciation of isolated clusters of hot fluid iron at low densities and high temperatures, representing a good approximant to the structure of the gaseous iron. In (a) and (b), the relative abundance and the average lifetimes of isolated iron atom are displayed. In (c) and (d), we have shown the abundance and lifetimes of iron dimer (Fe 2 ) represented by the circles, and iron trimer (Fe 3 ) in squares. Other clusters are not shown as their abundance are below 0.005. In the supercritical regime (10000 K and 12000 K), the iron fluid is dominated by long-lived monoatomic and diatomic clusters; at low temperatures we retrieve more familiar polymerized melts.

Figure 3 . 4 :

 34 Figure 3.4: (a) The overall trend of velocity autocorrelation function from 3000 K and 7.75 g/cm 3 to 10000 K and 3.86 g/cm 3 . (b) Velocity autocorrelation function along the supercritical isotherm 12000 K as a function of density.

Figure 3 . 5 :

 35 Figure 3.5: Vibrational spectrum of fluid iron calculated by Fourier transform of the velocity autocorrelation function shown in Fig. 3.4(a). At 3000 K and 7.75 g/cm 3 , we observe a broad peak from 0-400 cm -1 , which is caused by breathingtypes vibrations of the first coordination shell.

Figure 3 . 6 :

 36 Figure3.6: The mean-squared displacement (MSD) of fluid iron at selected densities and temperatures. The inset shows a log-log zoom-in plot for short times. The inset outlines a slope change that separates the two characteristic regions of the diffusion: ballistic, for times less than about 100 fs, and diffusive, at higher times.

  et al. (2019) Li et al. (2020) 3000K, Alfe et al. (2000) 6000K, Alfe et al. (2000)7000K,Alfe et al. (2000)

Figure 3

 3 Figure 3.8: Viscosity of fluid iron varying as a function of temperature and density. (a) Stress autocorrelation function for fluid iron at 3000 K and 7.75 g/cm 3 calculated by averaging over three off-diagonal stress tensor and normalizing by dividing the value at time zero. The inset shows the integral of the average stress autocorrelation function as a function of time. The shaded areas denote the fluctuation of the viscosity. (b) The overall overview of viscosity of fluid iron as a function of temperature and density. The black square represents an experimental value at 2500 K and 6.40 g/cm 3[START_REF] Assael | Reference Data for the Density and Viscosity of Liquid Aluminum and Liquid Iron[END_REF], which is close to our 3000 K isotherm. The black solid circles are the viscosity at 3000 K and 7.75 g/cm 3 determined from the Stoke-Einstein relation, which is in a fair agreement with the prediction by the Green-Kubo method.

  Figure 3.9: Electrical conductivity (a), thermal conductivity (b) and Lorentz number (c) of the fluid iron as a function of density and temperature. The Lorentz number is defined as L = k/σT . The results from our work are represented as solid circles colored based on temperature.The theoretic results at around 3500 K and 6.22 g/cm 3 are from[START_REF] Korell | Paramagneticto-Diamagnetic Transition in Dense Liquid Iron and Its Influence on Electronic Transport Properties[END_REF] shown as black squares. The black dashed line in (c) represented the ideal Lorentz number of 2.44.

Figure 3

 3 Figure 3.10: (a) Pressure as a function of density along several isotherms for liquid iron. (b)We have calculated the principal Hugoniot line with a initial condition of 300 K and 7.874 g/cm 3 , which compares well with theoretical simulation results by[START_REF] Sjostrom | Quantum molecular dynamics of warm dense iron and a five-phase equation of state[END_REF] and shock experiments.

Figure 3

 3 Figure 3.11: Computed Hugoniot lines for iron starting from two realistic initial states. (a) Temperature -pressure plots for various impact scenarios. The principal Hugoniot[START_REF] Sjostrom | Quantum molecular dynamics of warm dense iron and a five-phase equation of state[END_REF]) that starts at ambient conditions and the melting curve of iron[START_REF] Bouchet | Ab initio equation of state of iron up to 1500 GPa[END_REF] are shown for reference only. For a smaller impactor like Moon, we approximate the initial state conditions as 1 GPa and 1500 K; for a large impactor like Mars, the initial state is set at 40 GPa and 4000 K. The Hugoniot lines cross because the gains in temperature and pressure is not linear with respect to changes in initial conditions. The shaded area represents estimated temperature gradients ranges in the metallic cores of the different objects involved in the impact[START_REF] Antonangeli | Toward a mineral physics reference model for the Moon's core[END_REF][START_REF] Hirose | Composition and State of the Core[END_REF][START_REF] Stewart | Mars: A New Core-Crystallization Regime[END_REF]. (b) Computed entropy along the two Hugoniot lines. The star indicates the experimentally estimated entropy of boiling liquid iron at 3000 K, marked also by the dashed line.

  Figure3.12: Graphical representation of the impedance match analysis for the core of the small impactor. The shock wave generated at the moment of the impact travels through the impactor's mantle assumed to be homogeneous. The shock state at the core-mantle is given by the MgSiO 3 principal Hugoniot line (green line) estimated based on previous ab initio simulations[START_REF] Militzer | Equation of state calculations of hydrogen-helium mixtures in solar and extrasolar giant planets[END_REF]. More recent experimental Hugoniot data points are represented by crosses[START_REF] Fratanduono | Thermodynamic properties of MgSiO 3 at super-Earth mantle conditions[END_REF]. At the core-mantle boundary the shock wave splits in two opposite waves. The one traveling forward enters the core. The final state in the core is given by the intersection between the re-shocked MgSiO 3 Hugoniot (red line) and the iron Hugoniot (blue line) with an initial state at 1 GPa and 1500 K (the warm Hugoniot in Fig.3.11). The inset shows the fraction of iron that would vaporize from the corresponding impactor's core as a function of the impact velocity.

  Figure3.13: Isentropic release of the iron core starting from various peak conditions attained during large impacts. The shocked material gains entropy that is conserved during the release. If this entropy is larger than that of the boiling point then partial vaporization may occur. We compute the entropy along the spinodal and the Hugoniot lines at different temperatures, which allows us to map the entropy increase during various impact scenarios. We find that the entropy of the critical point, at 21.8 k B /atom, is reached for impacts with peak pressures of 605 GPa in the case of small planetesimals or for peak pressures of 825 GPa in the case of Marssized objects. Peak pressures of 312 GPa are enough to provide entropy higher than the entropy at boiling of iron at 1 atm, i.e. 15.84 k B /atom[START_REF] Kraus | Impact vaporization of planetesimal cores in the late stages of planet formation[END_REF]. These conditions can be easily exceeded during giant impacts between the proto-Earth and Theia, but could also be reached in almost half of the impacts with planetesimals during the late veneer. However, the liquid-vapor dome is reached only if the density, and hence the pressure, is allowed to decrease sufficiently. This can happen if the mantle is stripped away when fragments of the core are allowed to decompress without the mantle confinement.

cFigure 4 . 1 :

 41 Figure 4.1: Three different Monte Carlo moves in the Gibbs ensemble method. (a) A random displacement of a randomly chosen particle in each box. (b) A random volume rearrangement. (c) The transfer of a randomly chosen particle between the two sub-systems. (see main text for more details)

Figure 4 . 2 :

 42 Figure 4.2: Evolution of the number of particles (N ), cell length (L), pressure (P ), internal energy (E) and effective potential energy (U ) as a function of Monte Carlo steps in the vapor and in the liquid phase for the simulation at 2000 K.

Table 4 . 1 :

 41 Thermodynamic averages for the liquid and vapor phases for the ab inito Gibbs ensemble simulations. The uncertainties correspond to the one-sigma error bar.

Figure 4

 4 Figure 4.3: (A) Density versus Monte Carlo steps for the Gibbs ensemble simulations at 1200, 1500, 1800, 2000, 2100, 2200, 2300, 2400 and 2500 K from top to bottom, respectively. (B) Corresponding unnormalized probability distribution functions.

Figure 4 . 4 :

 44 Figure 4.4: Liquid-vapor equilibrium of Na obtained from the ab initio Gibbs ensemble simulations (blue circles) and its comparison with the experiment[START_REF] Dillon | Measurement of Densities and Estimation of Critical Properties of the Alkali Metals[END_REF] (black stars). The blue line is a fit of the equilibrium line using the scaling law and the red line is the law of rectilinear diameter with A = 0.80 ± 0.02 and B = 0.07 ± 0.01 parameters (see text for more details).

Figure 4 .

 4 5 shows the variation of pressure as a function of the density along the 2000 K isotherm. This curve exhibits

Figure 4 .

 4 Figure 4.5: Pressure evolution during the isothermal volume expansion for sodium at 2000 K. The insets show the snapshots at 0.52 and 0.38 g/cm 3 , respectively.

Figure 4 .

 4 Figure 4.6: Clausius-Clapeyron plot of the logarithm of the saturated vapor pressure and density as a function of the inverse temperature. The solid black squares are experimental date[START_REF] Dillon | Measurement of Densities and Estimation of Critical Properties of the Alkali Metals[END_REF]. The dashed line on the right graph denotes the pressures from fitting vapor density to the ideal gas law. As expected, a good agreement with calculated pressure is achieved in the low density range. In the high density region, a deviation is observed.

Figure 4 . 7 :

 47 Figure 4.7: Radial distribution function (RDF) g(r) at 1200, 1500 and 2000 K in the liquid phase as computed with the Gibbs ensemble MC simulations (full lines) and with the MD one at 0.52 g/cm 3 , 2000 K (dotted line). The shaded areas correspond to our estimate of the one-sigma uncertainty. The curves were shifted for readability. The inset shows the coordination number as a function of the temperature.

  Contents 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.2 Simulation details . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 78 5.3.1 The liquid-vapor equilibrium line from ai GEMC simulations 78 5.3.2 The critical point . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.3.3 The bulk modulus . . . . . . . . . . . . . . . . . . . . . . . . 81 5.3.4 The liquid density at zero pressure from ab initio MD simulations 84 5.3.5 The structure of the fluid . . . . . . . . . . . . . . . . . . . . 87 5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

5. 3 . 1

 31 The liquid-vapor equilibrium line from ai GEMC simulations 5.3.1.1 The nonmagnetic case

Figure 5

 5 Figure5.1: (a) The evolution of density as a function of Monte Carlo step from the aiGEMC simulations for iron in the nonmagnetic state at 3500, 5000, 6500, 7500, 8500 and 9500 K from top to bottom, respectively. (b) Corresponding unnormalized probability distribution functions. An extra simulation has been performed at 3500 K to examine whether the simulation results depend on the initial condition. By running the aiGEMC simulation with 108 iron atoms at 7500 K, we found the finite size effects are small at this condition.

Figure 5 . 2 :

 52 Figure 5.2: (a) The evolution of density as a function of Monte Carlo step from the aiGEMC simulations with iron in the ferromagnetic state at 3500, 5000, 6500, 7500, 8500, and 9500 K from top to bottom, respectively. (b) Corresponding unnormalized probability distribution functions.

Figure

  Figure 5.3:(a) The liquid-vapor equilibrium line for iron obtained from the aiGEMC simulations. The blue and green color denote simulations for iron in the nonmagnetic and ferromagnetic state, respectively. The solid circles represent the resulted liquid or vapor densities from the aiGEMC simulations, while empty triangles are densities determined by the law of rectilinear diameter. The blue and green line is a fit of the equilibrium liquid and vapor densities to the scaling law and the law of rectilinear diameter. The predicted critical point is represented by the solid stars. (b) The Clausius-Clapeyron plot of the logarithm of the saturated vapor pressure obtained from the aiGEMC simulations as a function of temperature.

Figure 5 . 4 :

 54 Figure 5.4: (a) The pressure-density fluctuation in the aiGEMC simulations for iron in the nonmagnetic (NM) and ferromagnetic state (FM) at 3500, 5000, 6500 and 7500 K. (b) The bulk modulus of the liquid iron phase defined as ρ dP dρ , where dP dρ is derived by a linear fitting of the pressure-density data.In comparison, we also performed aiMD simulations for the liquid iron in the nonmagnetic, ferromagnetic and paramagnetic state (PM) to find the zero pressure density and estimate its bulk modulus.

Figure 5

 5 Figure5.6: The constant of motion ((a) and (d)) as defined in Eq. 2.13, the effective potential energy (solid squares in (b) and (e)) as defined in Eq. 2.11 and the total magnetization (solid lines in (c) and (f)) of the system as function of time from the aiMD simulations starting with two different initial configurations at 3500 K at 6.63 g/cm 3 . In both simulations, the system was initialized to be ferromagnetic at the first MD step and the wavefunction and charge density extrapolation scheme as implemented in VASP was employed for the following MD steps. We selected several snapshots and performed the electronic minimization by re-initializing the system to be ferromagnetic. The resulted effective potential energy and the total magnetization (empty circles in (b) and (c)) is lower and higher than the one calculated with the wavefunction extrapolation scheme, respectively. It suggests that for the spin-polarized simulations, the results of the electronic minimization depends on the initial trial wavefunction or charge density.

  Figure5.7: The constant of motion (a) and the total magnetization (b) of the system as function of time for FM-aiMD simulations at 3500 K and 6.50 g/cm 3 represented by the blue line, and 6.41 g/cm 3 represented by the green line, respectively. In both simulations, the system was initialized to be ferromagnetic at every time step.

  Figure5.8: The calculated force for one randomly selected atom in the snapshot, pressure, the local magnetic moment for the selected atom, and the total magnetization of the different magnetic configurations. The accumulated average of the generated configurations is shown with a solid line.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

Table 2 .

 2 1: A summary of symbols used in Section 2.2. Ne the total number of electrons in the system {R}N a set of nuclear coordinates {R0, R1, .., RN } {r}N e a set of electronic coordinates {r0, r1, .., rN e } Ψ({R}N , {r}N e , t) the total wavefunction describing the dynamics of electrons and nuclei H fast the Hamiltonian governing the fast degrees of motion in H H slow the Hamiltonian governing the slow degrees of motion in H Φ l ({r}N e ; {R}N ) the eigenfunction of H fast at fixed {R}N E

	H	the total Hamiltonian operator of the system
	Te	the kinetic energy operator of electrons
	Ti	the kinetic energy operator of ions
	Ve-e	the electron-electron Coulomb interaction operator
	Ve-i	the electron-ion Coulomb interaction operator
	Vi-i	the ion-ion Coulomb interaction operator
	N	the total number of nuclei in the system

l ({R}N ) the eigenvalue of H fast at fixed {R}N l the electronic state (For instance, l = 0 denotes the ground state E0({R}N )) χ l ({R}N , t) the expansion coefficient to expand Ψ({R}N , {r}N e , t) by {Φ l ({r}N e ; {R}N )} F0({R}N )

Table 2 .

 2

2: A summary of symbols used in Section 2.3.

F0

same as F0({R}N ) but omitting its dependence on {R}N Γ a density matrix to describe the electronic subsystem at high temperature {ψ l } an arbitrarily complete orthonormal basis set in the Ne particle Hilbert space p l

3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Chapter 1

The formation of the Moon