Structural and Thermodynamic Properties of Magnesium-rich Liquids at Ultrahigh Pressure

: We explore the structural properties of Mg, MgO and MgSiO 3 liquids from ab initio com-1


13
Magnesium-rich liquids are significant minerals in planetary science because they rep-14 resent the main constituent of magma oceans [1] that formed when iron and silicates phase 15 separated during the formation of Earth and super-Earth planets [2][3][4]. Sizable magma 16 oceans have been predicted to persist over long periods of time [5]. Understanding how 17 the properties magnesium-rich liquids change with pressure and temperature will thus 18 enable us to better constrain the models of formation and evolution of rocky planets [6]. 19 Magnesium-rich liquids also provide valuable information for inertial confinement fusion 20 (ICF) experiments, where materials are exposed to extreme conditions [7][8][9]] that trans-21 form liquids into warm dense matter and dense plasmas that are difficult to understand. 22 State-of-the-art laboratories, including the National Ignition Facility (NIF) at Lawrence 23 Livermore National Laboratory, the Omega laser at the University of Rochester, the Z 24 machine at Sandia National Laboratory, and the SLAC laboratory, regularly investigate 25 these conditions to explore matter in the high-energy density regime, exploring matter 26 in the high-energy density regime by probing materials at ultrahigh pressure and tem-27 perature conditions never explored before, but that are present at the interior of white 28 dwarf stars [10], rocky exoplanets [11,12], and giant planets like Jupiter [13,14]. While these 29 experiments have accessed the warm dense matter regime and were able to achieve fusion 30 conditions recently [15], theoretical methods have difficulties accessing this regime due to 31 the strong level of ionization, which makes interpretation of experimental measurements a 32 hard task. 33 The properties of magnesium-rich liquids are very interesting. Their Grüneisen pa- 34 rameter increases upon compression [16] and the heat capacity can increase beyond the 35 ideal gas limit upon ionization [17]. MgO, SiO 2 and MgSiO 3 are all insulators in the solid 36 phase but they all become modest electronical conductors in liquid form at high enough 37 temperature and pressure, which implies that super-Earth planets can generate magnetic 38 fields in their mantle [5,18]. The structure of liquid magnesiosilicates varies substantially 39 with compression [19][20][21]. Some melts, such as MgSiO 3 , are good solvents for water. The 40 solubility was found to increase with pressure, which has important consequences for the 41 Earth's mantle [22]. At high temperatures, it is also important to identify the thermal and 42 pressure ionization regimes [23], as they modify the properties of the liquid as it becomes 43 partially ionized. Recent ab initio calculations predicted mixed coordination numbers for 44 silica between the 6-fold coordinated pyrite-type phase and the 9-fold coordinated Fe 2 P-45 type phase at high pressure, as well as the formation of superoxides such as SiO 3 and SiO 6 , 46 which may exist in the mantle of super-Earth planets [24]. 47 In this study, we perform ab initio simulations of the liquids Mg, MgO, and MgSiO 3 48 using a combination of density functional theory molecular dynamics (DFT-MD) and path 49 integral Monte Carlo (PIMC), providing a detailed characterization of the structure of the 50 liquids.

52
Rigorous discussions of the PIMC [25][26][27] and DFT-MD [28][29][30] methods have been 53 provided in previous works, and the details of our simulations have been presented in 54 some of our previous publications [17, [31][32][33]. Following earlier publications on hydrogen 55 and helium, PIMC and DFT-MD simulations have been combined to study the properties 56 of materials with core electrons in the regime of warm dense matter such as lithium 57 fluoride [34], boron [35], aluminum [31], oxygen [36], silicon [37,38], hydrocarbons [38,39], 58 and superionic water [40,41]. However, the structure of magnesium liquids has not been 59 explored in detail in this regime of extreme conditions. We combine PIMC [42] and and DFT-60 MD simulations as implemented in the Vienna Ab initio Simulation Package (VASP) [43]  For DFT-MD simulations, we employ Kohn-Sham DFT simulation techniques as 64 implemented in the Vienna Ab initio Simulation Package (VASP) [43] using the projector 65 augmented-wave (PAW) method [44,45], and molecular dynamics is performed in the 66 NVT ensemble, regulated with a Nosé thermostat. The time step was adapted to the 67 density and the temperature, ranging from 0.16 to 0.44 fs for simulation times from 1000 68 to 16 000 time steps, to ensure a reliable estimation of the thermodynamic quantities. The 69 pseudopotentials used in our DFT-MD calculations freeze the electrons of the 1s orbital 70 (He-core), which leaves 10, 12, and 6 valence electrons for Mg, Si, and O atoms, respectively. 71 Exchange-correlation effects are described using the Perdew, Burke, and Ernzerhof [46] 72 (PBE) generalized gradient approximation (GGA). However, for elemental Mg, the provided 73 Mg PBE pseudopotential did not give proper results for high densities, so we switched to 74 the local density approximation (LDA). We proceed in a similar way with MgO, the highest 75 densities were also simulated using the LDA functional. As shown in Ref. [47], the choice 76 of the pseudopotential in ab initio simulations of Mg has very little effects on the computed 77 thermodynamic properties. We obtain a very good agreement between both functionals 78 for a number of densities. Electronic wave functions are expanded in a plane-wave basis 79 with a energy cut-off as high as 7000 eV in order to converge total energy. Size convergence 80 tests with up to a 65-atom simulation cell at temperatures of 10 000 K and above indicate 81 that pressures are converged to better than 0.6%, while internal energies are converged 82 to better than 0.1%. We find, at temperatures above 500 000 K, that 15-atom supercells 83 are sufficient to obtain converged results for both energy and pressure, since the kinetic 84 energy far outweighs the interaction energy at such high temperatures [31,48]. The number 85 of bands in each calculation was selected such that orbitals with occupation as low as 86 10 −4 were included, which requires up to 14 000 bands in an 15-atom cell at 2 × 10 6 K and 87 two-fold compression. All simulations are performed at the Γ point of the Brillouin zone, 88 which is sufficient for high temperature fluids, converging total energy to better than 0.01% 89 compared to a grid of k-points.

91
In this section, provide a detailed characterization of the magnesium-rich liquids that 92 we have obtained from our ab initio simulations of Mg [33], MgO [32], and MgSiO 3 [17], 93 spanning the condensed matter, warm dense matter and plasma regimes. Computations 94 were performed for a series of densities and temperatures ranging from 0.321-86.11 g cm −3 and 95 10 4 -10 8 K. The full range of our EOS data points is shown in temperature-density and 96 pressure-density space in Fig. 1, along with the shock Hugoniot curve of each material.  The thick dashed lines correspond to the shock Hugoniot curves that we derived for Mg, MgO, and MgSiO 3 , with initial densities ρ 0 = 1.736 g cm −3 , 3.570 g cm −3 , 3.208 g cm −3 , respectively. The full EOS for each material is available in Refs [17,23,32,49]. Shock Experiments on MgO from McCoy et al. [50], and Root et al. [51], and shock experiments on MgSiO 3 from Fratanduono et al. [52], are included for comparison.

97
The ideal mixing approximation has been shown to perform well for temperatures 98 above 10 5 K. [53], and the magnitude of nonideal mixing effects was found to be small in 99 this regime, leading to shock Hugoniot curves of MgO and MgSiO 3 that are reproduced 100 with sufficient accuracy by combining the EOSs of the elemental substances with the addi-101 tive volume rule. This concept was extended to other mixtures [49] and good agreement the 102 shock Hugoniot curves of H 2 O and CO 2 was found between laboratory measurements and 103 theoretical predictions based on the linear mixing approximation. However, this approxi-104 mation breaks down at lower temperatures, where chemical bonds play an increasingly 105 important role. These bonds change the structure of the liquid, modifying the atomic 106 coordination. Here we will study how the structure of the liquid changes with density 107 and temperature. We will demonstrate that coordination of magnesium ions is sensitive to 108 presence of silicon and oxygen. The atomic trajectories obtained from DFT-MD simulations can be used to study the 111 local structure of the liquids. Using the radial distribution function, defined by we can obtain a measure of the structure of the liquid, which depends on temperature and 113 density. Here, N α and N β are the total number of nuclei of type α and β, respectively that 114 are contained in the volume V, while r ij = r i − r j is the separation between nuclei i and j. 115 N β (r) is the total number of nuclei of species β within a sphere of radius r around a nuclei 116 of type α [54,55]. This function, g αβ (r), can be interpreted as the probability of finding a 117 particle of type α at a distance r from a particle of type β.
In Fig. 2, we compare radial distribution functions of liquid magnesium over a wide 119 range of temperatures and densities. For every density, we observe that the average distance 120 to the nearest neighbor, given by the location of the first maximum of the g(r), slightly 121 decreases with increasing temperature. At 6.89 g cm −3 , this distance shifts from 1.9 Å at 122 20 000 K to 1.8 Å at 100 000 K. This 5% decrease is caused by stronger collisions and higher 123 kinetic energy. As expected, the distance to the nearest neighbors depends strongly on 124 density. At the highest density explored (51.67 g cm −3 ), this distance decreased to 0.96 Å. 125 This value does not depend much on temperature. For the three highest densities in our simulations (ρ 34.44 g cm −3 ), the system freezes into an simple cubic structure at 20 000 K, which corresponds to the stable phase of Mg observed experimentally at pressures exceeding 1 TPa, the highest pressures ever reported in experiments of Mg to date [56] . The arrows indicate the location of an emerging intermediate local minimum. 126 We can also observe in Fig. 2 that the shape of the all g(r) curves at a low densities 127 is fairly similar. All curves have two well-defined maxima and two minima. However, 128 for densities of 19.37 g cm −3 and higher, the liquid becomes significantly more structured. 129 The curve develops an additional intermediate local maxima that result in a new local 130 minimum at ∼ 1.55 Å. As density increases, this intermediate maximum becomes more 131 pronounced. In our simulations at 20 000 K, the system freezes into the simple cubic phase 132 for densities of 34.44 g cm −3 and higher, which corresponds to pressures higher than 133 24381 GPa. This simple cubic phase has been observed experimentally at pressures of 134 1 TPa [56]. With crystal structure search methods, a series of other high-pressure structures 135 including body-centered cubic, face-centered cubic, and simple hexagonal and simple cubic 136 phases have been predicted for magnesium [57]. When this crystallization happens in our 137 MD simulations, the subtle features of the liquid g(r) functions become amplified. But for 138 temperatures of 30 000 K and above, we found that the system remained in a liquid state 139 for all densities. Still the local minimum at r min ∼ 1.55 Å persisted. We will discuss these 140 changes in terms of the atomic coordination number in the next section.

142
A further measure of the structure of the liquid is the coordination number, given by 143 where C αβ is defined to be the number of atoms of type β that are within a spherical region 144 of radius r min , centered at an atom of species α [54]. Here, we adopted the usual convention 145 that r min is the location of the first minimum of the radial distribution functions g αβ (r). 146 The integrated nucleus-nucleus pair correlation function, given by N β (r) in Eq.
(1), can be 147 employed to define a coordination number in Eq.
(2) when evaluated at the location of the 148 first minimum, that is, N β (r min ) = C αβ .

149
In Fig. 3, we show the radial distribution functions of liquid elemental magnesium at 150 30 000 K and the corresponding integrated nucleus-nucleus pair correlation functions, N(r). 151 The shaded areas in the top panel highlights the interval [0, r min ] over which the integration 152 in Eq. (2) is performed. As density increases, the position of the first local maximum shifts 153 to smaller distances and the peak becomes narrower, which indicates that the distance to 154 first neighbors is decreasing. r min , also decreases with increasing density, reducing the total 155 area below the curve and, hence, the coordination number.

156
As the atoms get closer, the liquid develops a new structure, which is reflected in 157 the intermediate local maximum that starts developing in the g(r) for densities higher 158 than 25.83 g cm −3 at 30 000 K, as we can see in Fig. 3, which corresponds to a pressure 159 of 13147 GPa. The typical environment of a magnesium atom, after this intermediate 160 local maximum develops, is shown at the top of Fig. 2, where a new shell of nearest 161 neighbors forms. The slope (dE/dρ) T is positive, therefore pressure ionization is likely to 162 take place [23]. As described in the previous section, a transition occurs in the liquid for 163 higher densities, where a new intermediate maximum appears, which abruptly decreases 164 the value of r min and, hence, the coordination number. 165 In Fig. 4, we show the resulting Mg-Mg coordination number in liquid magnesium 166 for a wide range of temperatures and densities. As we can observe in the figure, the 167 coordination number first increases from 14 to 18 at T = 20 000 K and then abruptly drops 168 to 7, consistent with the diagram shown in Fig. 2. Something similar occurs at T = 30 000 K, 169 where the coordination number drops from 18 to 8. This is an indication that the system 170 prefers a coordination similar to the simple cubic structure, where the coordination number 171 is 6. For higher temperatures, this intermediate maximum never develops for this range of 172 densities, but it is likely to appear for ρ > 60 g cm −3 . For T > 30 000 K, the coordination 173 number remains between 13 and 18 at all densities. g(r) 6.89 g/cm 3 8.61 g/cm 3 12.92 g/cm 3 19.37 g/cm 3 25.83 g/cm 3 43.06 g/cm 3   the eigenvalues towards lower energies. This broadening does not affect the electrons in 182 the conduction band (3s) in a significant way. However, at 21.53 g cm −3 , the broadening 183 of the other bands is significant and almost closes the gap between the 2s and 2p bands. 184 According to Ref. [33], these conditions fall in the regime of pressure ionization where the 185 effects of thermal ionization are expected to be small. The ELF depicts location of electronic 186 charge in the between void between the nuclei, which is the typical electride behavior that 187 results from the repulsion of core electrons. For the simple cubic phase of Mg, this has 188 recently been reported by a experimental-theoretical study [56]. With ab initio methods, 189 this has been predicted to occur in many solid structures at high pressure [59]. In Fig. 5, we 190 show that also the Mg liquid exhibits electride behavior, which shares similarities to earlier 191 predictions for liquid iron [60]. In the solid phase (right), the electronic charge (yellow pockets) is localized in the voids between the nuclei, which depicts the typical electride behavior that has been predicted with ab initio method to occur in many solid structures including K [58, 61,62] and Mg [56]. Here we also find the Mg liquid to exhibit electride behavior (left), which has been predicted to occur in liquid Fe [60]. The Mg-Mg coordination number is much higher in liquid, elemental Mg (∼14) than in 202 liquid MgO. In the presence of oxygen, the Mg-Mg coordination number decreases from 8 203 to 4-6 with increasing density, as we observe in the right panel of Fig. 6. However, this lower 204 coordination is comparable to that of elemental Mg (∼7-10) if the density of that system is 205 high and the temperature is low, as we showed in Fig. 4. The Mg-O coordination number 206 was found to be around 7 for all temperatures and densities, consistent with previous 207 studies of MgO at high temperatures and densities [19]. Similarly, the O-O coordination 208 number remained approximately 6. This is smaller that the value found for the molecular 209 fluid GeO 2 liquid, where the average O-O coordination number was found to be around 9 210 and a Ge-O coordination number of 2 at 1500 K and low temperatures [63,64]. This is not 211 the case of MgO, which behaves as an atomic fluid. This is comparable to silica, which is 212 predicted to a have a mixed coordination between the 6-fold coordinated pyrite-type phase 213 and the 9-fold coordinated Fe 2 P-type phase at high pressure [24]. Therefore, coordination 214 numbers with oxygen between 6 and 9 are expected at these conditions.

215
In the case of MgSiO 3 , the Mg-O coordination number is also between 6 and 8, as we 216 can see in Fig. 7, but slightly larger than 8 at some densities. The Mg-Mg g(r) functions in 217 the left panel show that magnesium atoms are not correlated, with no clear signature of 218 a local minimum at most of the conditions that allows to identify a layer of first nearest 219 neighbors. When the identification is possible, the Mg-Mg coordination number in MgSiO 3 220 lies between 2 and 4. While the Mg-O coordination number seems to decrease with 221 density, the Si-O coordination number increases with increasing density, regardless of the 222 temperature, from 6 to 8. The O-O coordination number is larger in the MgSiO 3 liquid 223 than in the MgO liquid, reaching values between 7 to 9 in the latter case, while MgO 224 shows an average of 6 for most conditions. Regarding the distance to nearest neighbors, or 225 bond length, we find that the Mg-O bond length is larger in MgSiO 3 than in MgO. For a 226 comparable denisity of ∼25 g cm −3 , this distance is 1.3 Å in MgSiO 3 , while for MgO, this 227 distance is about 1.2 Å. At all conditions, the Si-O bond length is smaller than the Mg-O, 228 and both are larger than the O-O bond length. As in the case of MgO, there no signature of 229 an abrupt structure transition, as we observed in pure Mg.

231
We have studied the structural properties of magnesium liquids in a broad range of 232 temperatures and densities using ab initio simulations. We found evidence of a structural 233 transition in liquid Mg around 20 g cm −3 , where the emergence of a new intermediate 234 maximum in the radial distribution function leads to an abrupt decrease in the coordination 235 number with increasing density for low temperatures. The structural change in liquid 236 elemental magnesium is an indication of a transition to electride-type behavior, consistent 237 with recent experimental finding of electride phases of Mg at ultrahigh pressure [56].