Structure, equation of state, diffusion and viscosity of warm dense Fe under the conditions of giant planet core

Fe exists abundantly in the universe. In particular, the dynamical structures and transport properties of warm dense Fe are crucial for understanding the evolution and structures of giant planets. In this article, we present the ionic structures, equation of states, diffusion and viscosity of Fe at two typical densities of 33.385 g/cm$^3$ and 45 g/cm$^3$ in the temperature range of 1 eV and 10 eV, giving the data by the first principles calculations using quantum Langevin molecular dynamics (QLMD). Furthermore, the validation of Stokes-Einstein (SE) relation in this regime is discussed, showing the importance of choosing the effective atomic diameter. The results remind us of the careful usage of the SE relation under extreme conditions.


Introduction
The properties of complex materials such as Fe under extreme conditions are crucial for understanding the evolution of planets (earth, giant planets, extro-solar giant planets) [1] and stars such as sun [2,3]. More interestingly, the laser-driven and Z-pinch experiments can lead to very high densities up to hundreds of g/cm 3 and temperatures up to thousands of eV [4,5]. In particular, the temperatures of 1-10 eV with associated pressures of 10-1000 TPa are likely to exist in the interiors of massive planets [6,7]. Recent studies have shown the possible stable structures at these pressures [6,7], where face-centered cubic (fcc) structures are stable in the range of 7-21 TPa. However, the effect of temperatures on these structures is still not known, since the dynamical structures such as melting, diffusion and viscosity can be induced by temperatures, and the phases of Fe in earth-like exoplanets are likely to be liquid. In order to understand these behaviors, the experimental determination seems currently extremely expensive and difficult. Therefore, accurate simulations are required for the determinations of the structures and dynamics. The dynamical structures of Fe and its compounds from firstprinciples molecular dynamics (FPMD) at the physical conditions of the Earth's core have been studied widely [8,9,10], showing the complexity of the new phenomena in high energy density physics (HEDP). In particular, it is worthy to verify whether the Stokes-Einstein (SE) relation at so high densities is valid, since SE relation is one important bridge between dynamics and statistics, and usually adopted in the previous calculations of the viscosities from diffusion coefficients [8,9,10,11,12]. Up to now, the diffusion and viscosity can be obtained directly by molecular dynamics simulations with Kubo-Green relation [10,11,12,13]. Moreover, we can construct statistical models according to the simulation results, introducing easy applications in hydrodynamics. [14,15], However, molecular dynamics simulations are strongly dependent on the accuracy of many-body interatomic potentials, and there are usually some limitations for the statistical methods.
Cold stable fcc structure with its inducing liquid at high temperature is considered as an important phase in the giant planets such as Jupiter and exoplanets [6,7]. In this regime, the densities are from 33.9 g/cm 3 to 66.7 g/cm 3 , and the temperatures are from 1 eV to 10 eV. For these states, the melting behaviors and the transport properties are interesting and deserved to be studied carefully. Besides, much high density will introduce much different chemical bonds since the pressure larger than 1 Mbar will change the traditional chemistry dramatically [16,17,18,19,20]. Furthermore, when the pressure increases to 100 Mbar, the core electron charge density can be changed significantly [16,17,18]. Therefore, the transport properties such as equation of states (EOS), diffusion and viscosity will change a lot since these properties are related to the electronic structures. With respect to the methods of calculating these properties accurately, the first principles calculations are thus required because we did not understand their behaviors from semiclassical models such as average atom (AA) model [21,22] and Thomas-Fermi (TF) model or orbital free (OF) method [23,24,25], in which the orbital behaviors or chemical bonds information can not be described well.
Molecular dynamics combining finite-temperature density functional theory (DFT) [26] called quantum molecular dynamics (QMD) and quantum Langevin molecular dynamics (QLMD) have been successfully used in warm and hot dense matter, especially under very high pressures [8,10,18,19,20,27,28]. The validations of QMD and QLMD at high densities and temperatures have been verified a lot of times by comparing with experiments and interestingly reasonable results are found. In particular, QMD is considered as one of the most promising methods and has been successfully used to predict the transport properties of warm dense matter (WDM). [29,31] The biggest argument for QMD applications is that its limitation on the number of particles because of its computational cost at high temperature. However, QMD and QLMD can be efficient for the calculations of WDM under the conditions of giant planet core, and give relatively accurate transport properties using a few hundred particles. [29,30,31,32] In this article, we adopted QLMD to simulate the transport properties of dense Iron at 45 g/cm 3 from 1 eV to 10 eV. Except for, considering the density of 33.385 g/cm 3 at 100 eV along the Hugoniot curve [18], we also calculate the dynamical properties at 33.385 g/cm 3 up to 10 eV. Furthermore, the approximation of SE relation is discussed, showing the debatable validation under so dense regimes.

Methods
For the complexity of Fe at high density, both QMD with velocity rescaling method and QLMD simulations are tested firstly. In fact, at low temperature, the electronion collisions induced friction (EI-CIF) is small since the friction is proportional to the temperature. Therefore, QMD and QLMD are equivalent for the equilibrium properties at low temperature. Eventually, the advantage of QLMD is the higher computational efficiency than traditional QMD [33]. The electronic structures are solved based on finite temperature DFT implemented in the Quantum ESPRESSO package [34]. Considering the small atomic sizes at high density, the time step of 0.5 fs is used in order to keep the correct trajectories of ions. Pseudopotential (PP) is one of the key points in QMD and QLMD simulations. Here, we construct a new PP with 16 electrons in the valence and 0.9 atomic units for the radius cutoff within the generalized-gradient approximation (GGA) [35], which can promise the correctness with respect to the conditions we are studying. Using this PP, we can reproduce the bulk modulus properties of Fe as previous results calculated by full-electron calculations and experiments. More importantly [36,18], at high pressure, the pressure and band structures at the density of 33.9g/cm 3 and 48.23g/cm 3 are in good agreement with the pressure in Ref. [7] using the same PP [18]. Furthermore, the same PP has been successfully used to calculate the EOS of Fe in hot dense regime within a wide range of densities and temperatures [18], indicating the validation of PP in this work. Besides, 2000 time steps with a large convergent tolerance of 1.0×10 −4 in electronic structure calculations are adopted in order to reach the thermalization, and 10 ps time lengths with a small convergent tolerance of 1.0×10 −6 [33] are simulated to obtain the transport properties and thermal average. In order to accelerate the calculations, the Gamma point is only used for the representation of Brillouin zone. For comparing with the results from semiclassical methods, averaged atom molecular dynamics (AAMD) [22,37] is performed for all cases. When two atoms are close enough, their electronic distributions would be overlapped, inducing interactions between them. In order to describe this interactions, we should establish the interatomic potential. In AAMD method, the temperature-dependent pair potential, which can describe the contributions of the overlapped electronic distributions, is constructed from the AA electronic calculations, including the ionic correlations within the framework of AA model. Within this framework, the ions are moving on this pair potential, giving the ionic trajectories and ionic configurations at different times. From AAMD, we can basically use large number of particles and giving the convergent results corresponding to the size effects. The comparison between QLMD and AAMD has been reported before [19,33], showing their agreement at high temperatures. The details of AAMD can be found in Ref. [37].
The self-diffusion coefficient D can be obtained from two methods. One is from the trajectories by the mean-square displacement D R where the R i is the position of the i th particle. The other method is from the velocity autocorrelation function D V where V i is the velocity of the i th particle. Besides, the viscosity η can be obtained from the autocorrelation function of the off-diagonal component of the stress tensor so called Green-Kubo equation [10,38] where P 12 represents the averaged result for the five independent off-diagonal components of the stress tensor P xy , P yz , P zx , (P xx − P yy ) / 2, and (P yy − P zz ) / 2. Using this method, the self-diffusion coefficient, viscosity of warm and hot dense matter have been reported widely [8,9,10,11,12,29,39,40,41,42] within the framework of QMD, showing its validation in these extreme conditions. On the other side, the viscosity can also be calculated directly from the SE relation where a is an effective atomic diameter. This relation is statistically obtained from the Brownian motion of a macroscopic particle in liquid, but it is only an approximation for the atoms. If the validation of SE relation can be verified, the size of particles, the transport behaviors can be understood well [43]. In fact, for a Brownian particle, SE relation is equivalent to Eq. 2 and Eq. 3. In particular, for the dense matter, the diameters of atoms are very small, and the validation of SE relation should be examined very carefully.

Convergence tests
For calculating the transport properties from first principles, the number of atoms should be tested since diffusion and viscosity are strongly dependent on the system sizes. For this purpose, in QLMD simulations, we tested face-centered cubic (FCC) structures with 32, 108 and 256 atoms and body-centered cubic (BCC) structure with 54 atoms in a supercell at 45 g/cm 3 . The pressure and diffusion versus temperature relations with different number of atoms are shown in Fig. 1. It is concluded that the pressure is not much sensitive to the system sizes. In particular, the pressures are closer for different sizes when the temperatures are higher. The small difference at low temperatures may be caused by the long ranged orders with crystal structures. It is worth noting that the initial configurations of BCC and FCC structures are important for the calculation of diffusion in the solid phases, but the differences from different initial configurations disappear gradually with the increasing temperatures in the liquid phases (here above 5 eV).
According to the tests, it can be known that our systems are really convergent for the calculation of pressure and diffusion coefficients when we use 108 atoms with initial FCC structures. It should be noted that when the system becomes liquid (after the jump in diffusion coefficient), the results of 54 atoms, 108 atoms and 256 atoms are consistent. The difference at low temperatures would be from the different initial structures of BCC or FCC crystals. In fact, the small differences from different sizes can also come from the statistical errors for small sizes. Basically, this error can be overcome from the long time simulations.
For the viscosity, the results of the convergent tests are shown in Fig. 2. It is only useful for the liquids here since the viscosity is very large for solids. Here, the liquid phases are shown above 5 eV in our simulations. By comparison, we can know that the viscosity are much sensitive to the sizes. For small size of 32 atoms, the statistical errors are really large which is not likely to be compensated by increasing the simulation time. When the number of atoms is increased up to 108 and 256, the viscosities at different temperatures seem convergent. Therefore, we can safely use the 108 atoms for the calculations of pressure, diffusion, and viscosity in our cases. In fact, about one or two hundred atoms are usually used in QMD simulations for the transport properties [8,10,12,29,30,31], which has been shown to be convergent within reasonable errors.
For the effect of the particle number, we use large number of particles up to 5324 in AAMD calculations to verify the convergence, as shown in Fig. 3. The simulation time is up to 100 ps in each case in order to obtain the correct correlations of particles in these semiclassical calculations. It can be shown that although the number of particles can affect the convergence of diffusion and viscosity, all the differences are within 10% or smaller. When the number is reach to 4000, the diffusion and viscosity are almost convergent. In order to compare the results, we use 4000 particles in AAMD to calculate the transport properties below. It is worth noting that AAMD adopted pair potentials, which should be more sensitive to the size and argued for the viscosity since the unsymmetrical shear of the system can not be described well by pair potentials.
For the density of 33.385 g/cm 3 , we used 54 atoms with initial BCC structures to calculate their transport properties. The simulation details are also the same as those of 45 g/cm 3 . Basically, the solid phase of Fe at 33.385 g/cm 3 may be hexagonal closepacked (HCP). However, the determinations of the lattice constants of HCP at high temperatures are very difficult, and the dynamical properties such as the structures and melting behavior are much similar to the BCC structures [7,45]. Therefore, we used BCC structures as the initial configurations. In fact, from the tests of Fe at 45 g/cm 3 , it is known that the transport properties are independent on the initial configurations  when Fe melts. Furthermore, the effects of larger system sizes on calculated viscosities have been previously studied for models such as hard-sphere and Lennard-Jones liquids [44], showing meaningful values for µ even with only 32 atoms. Thus, the initial configurations and sizes in this work can promise the reasonable accuracy of our simulations.

Structural dynamics
The dynamical structures for dense matter are basic to understand the melting behaviors and the dynamics in planets core. However, dense liquid Fe is rarely studied before due to the lack of effective methods. Dense Fe can hold the solid phase or ordered structures at very high temperatures, form clusters and chemical bonds assisted by core electrons [18], and keep high melting temperature [45,46,47]. Here, the structures are investigated by looking at the radial distribution functions (RDF) g(r), as shown in Fig. 4 and Fig. 5.
From the RDF, we can know that when the temperature increases up to 5 eV, Fe exhibits liquid behaviors at 45 g/cm 3 , and above 3 eV at 33.385 g/cm 3 . We calculated the details of the dynamical structures near 5 eV and 3 eV for the two densities respectively, using the two phase approach (TPA) [45] with 216 atoms and 108 atoms to obtain the melting temperature. The melting temperatures are found to around 57000 K and 30000 K, within the errors of 200 K in our simulations. The same behaviors can be also found in the data of diffusion coefficients and viscosities, which are shown in Fig. 1 and Fig. 2 and will be discussed detailedly later. In this stage, the Fe at 45 g/cm 3 and 40000 K is also in solid phase in the Fe phase diagram [7].
We should note that the first peaks of the RDF move toward to the zero point, indicating the effective distance between ions decrease with the increasing temperatures. Therefore, we may use the positions of the first peaks in RDF as the effective atomic diameters when discussing about the SE relation later. If we use only the averaged ionic  radius, the temperature effect can not be included since the radius is only dependent on the type of elements and system densities.

Transport properties
After the convergence tests, we can obtain the transport properties including EOS, diffusion coefficients and viscosities for dense Fe. The summaries of the data are shown in Table. 1 and Table. 2. In the tables, pressures and diffusion coefficients from AAMD are also shown. When the solid phases are not broken, the diffusion coefficients are very small and the viscosities are very large (therefore they are not convergent in Green-Kubo equation).
For the density of 45 g/cm 3 , from Table. 1, the diffusion coefficients change rapidly between the temperatures of 40000 K (3.447 eV) and 5 eV, and reach the orders of liquids diffusion at 5 eV. Also, the viscosities above 5 eV go to convergence when we integrate the pressure autocorrelations in Eq. 3, indicating the existence of liquid phases. However, AAMD gives the liquid Fe at 4 eV (not shown here), showing its limitations from pair potentials in the very strongly coupling regimes. Furthermore, we can see that the diffusion coefficients from Eq. 1 and Eq. 2 are equal, which shows the validation and convergence of our simulations. For the pressure, AAMD gives much larger pressures than those of QLMD simulations, which is caused by the electronic structure calculations from Thomas-Fermi methods. It is very interesting that when the temperatures are above 5 eV, the diffusion coefficients from QLMD and AAMD are very close. One reason might be that the diffusion coefficients are strongly dependent on the nearest neighbors distributions. Pair potentials can deal with the nearest neighbors interactions between ions with good accuracy. This can be verified using the comparisons of RDF from QLMD and AAMD, as shown in Fig. 6, where the positions of the first peaks of RDF are almost overlapped. Therefore, the nearest neighbors interactions should be reasonable in AAMD. However, AAMD can not deal with many-body interactions, inducing the shear viscosities from AAMD are incorrect here, since dense Fe liquid holds a lot of medium ordered structures, and therefore the melting behavior is not correct within AAMD framework. For the density of 33.385 g/cm 3 , the behaviors are almost the same. The melting temperature locates below 3 eV, according to the RDF and diffusion coefficients. It is noticed that the pressures from AAMD are much closer to those from QLMD compared with the results of 45 g/cm 3 . This means that the electronic structures from AA model are reasonable here but give more free electrons at 45 g/cm 3 . Furthermore, the melting temperature from AAMD is also lower than that from QLMD, and the diffusion coefficients are in the same order when Fe melts above 3 eV. Table 1. Summary of the results of dense Fe at different temperatures at 45 g/cm 3 . D R and D V are respectively defined in Eq. 1 and 2. D AAMD and µ AAMD are the diffusion coefficient and viscosity from AAMD calculations, respectively.    The small sizes in the calculations can introduce statistical errors. Especially for the viscosities, the errors are within 10%. Therefore, calculations with more than 10000 atoms at least and accurate many-body potentials are needed basically. We need to construct more accurate semiclassical potentials for obtaining more accurate diffusion coefficients and viscosities.

SE relation
With respect to the simulation time and system sizes, we find that the diffusion coefficients can reach convergence much easier than the viscosities. Therefore, if we can Table 3. Summary of the effective diameters and its corresponding viscosities for dense Fe at different temperatures at 45 g/cm 3 . Z is the average ionization degree from AAMD calculations. µ i , µ f , µ E are the viscosities corresponding to different definitions of effective atomic diameters of ions r i , r f and r E , respectively. µ Y V M is the viscosity from the model of Yukawa viscosity model (YVM) (see details in Ref. [14].) obtain the viscosities from the calculated diffusions, a lot of computational resources can be saved. Therefore, the validation of SE relation is important and meaningful. In order to test the validation, we calculate viscosities from different definitions of effective atomic diameters: averaged ionic distance at a specific density r i = V (1/3) (where V is the averaged volume of one atom at a specific density), the position of the first peaks in RDF (r f ), and the effective radius (r E ) from effective coordination number (ECN) definition. ECN can reveal the topology of the structures partly, and also give the local information around one atom. The definition of ECN is as following where d ij is the distance between atoms i and j; N is the total number of atoms in the system; d i av and its average value d av are defined as With this definition, the effective diameters r E = d av , which includes the temperature and density effects.
The results of different effective atomic diameters and the calculated viscosities from SE relation are shown in Table. 3 and Table. 4. Comparing the results with the calculated viscosities from Green-Kubo equation, we find that the viscosities are very different at low temperatures. With the increasing temperature, SE relation seems to be more and more validated. Moreover, the choice of the effective atomic diameters plays an important role in the validation of the SE relation, where definitions of the positions of the first peaks of RDF and diameters from ECN calculations are more appropriated. Considering the calculation uncertainties and statistical errors, the SE relations can be considered validated only when the choice of the effective diameters are reasonable.  In order to avoid the usage of SE relation, some models have been constructed. For example, Murrilo reported the construction of the Yukawa viscosity model (YVM) based on molecular dynamics simulations [14], which can be used in warm dense regime. We verify the validation of YVM model here. First of all, we give the average ionization degree (Z) from AA model with Electronic Energy-level Broadening [22], and then show the viscosities (µ Y V M ) from YVM model, as shown in Table. 3 and Table 4. It can be shown that YVM model improve the accuracy of SE relation, especially for the cases at high temperatures. This suggests that the Yukawa potential might be reasonable for the states at so high density at relatively high temperatures.
What's the possible reason for the debatable validation of SE relation in this regime? The basic physics is whether the ions exhibit as Brownian particles. For the specific case in this work, heavy Fe ions are moving in the liquid-like electron sea, randomly collided with many free-like electrons. However, the other dominant factor for the ionic motions are the interactions between ions, i.e., the strong coupling of ions. Every ion moves in a specific force fields determined by both electrons and ions, and this force field is not constant. The ions in dense matter have memory effect, i.e., the ionic positions at this time plays important role in the ionic movements at next time. In other words, if the correlation time of ions is larger than the observable time, the motion of the ions should not be Morkov process. In order to analyze the correlation time of the systems, the velocity autocorrelation functions (VAF) of ions at 5 eV and 10 eV are shown in Fig. 7. Here, we adopt the simplest definition of the correlation time (τ c ) to compare the correlation time at different temperatures. Here, τ c is defined as the time for which the VAF(t) decreases to 1/e [48]. With this definition, τ c is about 4.1 fs and 2.7 fs for the temperatures of 5 eV and 10 eV, respectively. Therefore, with the increasing temperature, the behaviors of ions are more and more Brownian-like motion. In fact, if the system is ideal gas, the correlation time will be zero, indicating uncorrelated behaviors between two near time steps. That is to say, when temperature is high enough, the random collisions between particles are dominant, and the statistical behaviors can work well. If the correlation time of the system is zero, i.e, we can recover SE relation directly from Eq. 2 and Eq. 3 [49]. We can then conclude that the SE relation is strongly dependent on the coupling parameters of ions Γ = Z * 2 /(k B T a), with T the system temperature, k B the Boltzmann constant, a the mean ionic sphere radius defined as a = (3/(4πn i )) 1/3 , Z * the average ionization degree, n i the ionic number density. If Γ is small enough, SE relation can be valid. Furthermore, one possible reason that QMD or QLMD can obtain the reasonable transport properties within relatively short time (within 10 ps) is the short correlation time of particles as shown in Fig. 7.

Conclusion
In conclusion, the dynamical structures and transport properties including EOS, diffusion and viscosities are calculated using QLMD method within the framework of first principles. The dynamical structures show the information of RDF, melting behaviors, diffusions and viscosities. The validation of SE relation is also discussed, showing the importance of the methods of choosing the effective atomic diameters. This work studies the dense liquid Fe existed in giant planets such as Jupiter, exoplanets, indicating the necessity of first principles calculations or constructing accurate manybody interactions. Furthermore, the results are crucial for understanding the laser-shock induced compression experiments, reminding us of the complexity of dense matters. In this field, a lot of new physics needs to be deeper studied.