{'fraction_non_alphanumeric': 0.04001978320900136, 'fraction_numerical': 0.02213246507027161, 'mean_word_length': 4.766159695817491, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
First-Principles Prediction of the Softening of the Silicon Shock Hugoniot Curve

Shock compression of silicon (Si) under extremely high pressures (>100 Mbar) was investigated by using two first-principles methods of orbital-free molecular dynamics (OFMD) and path integral Monte Carlo (PIMC). While pressures from the two methods agree very well, PIMC predicts a second compression maximum because of 1s electron ionization that is absent in OFMD calculations since Thomas-Fermi-based theories lack shell structure. The Kohn-Sham density functional theory is used to calculate the equation of state (EOS) of warm dense silicon for low-pressure loadings (P<100 Mbar). Combining these first-principles EOS results, the principal shock Hugoniot curve of silicon for pressures varying from 1 Mbar to above 10 Gbar was derived. We find that silicon is 20% or more softer than what was predicted by widely-used EOS models. Existing high-pressure experimental data (P = 1 - 2 Mbar) seem to indicate this softening behavior of Si, which calls for future strong-shock experiments (P>10 Mbar) to benchmark our results.


INTRODUCTION
Silicon, an important element widely used in the semiconductor industry [1], is closely related to many other scientific fields such as geophysics [2], planetary science [3], and astrophysics [4]. Moreover, silicon has also been used as a dopant to the ablator material for target designs in inertial confinement fusion (ICF) [5][6][7], as well as for mitigating laser-imprint [8,9] and two-plasmon-decay effects [10] in direct-drive ICF.
Although it is crucial to many scientific fields, ranging from understanding the geophysics of the Earth's outer core [2] to its application to high-energy-density physics (HEDP) and ICF, the properties of silicon at high pressures (>10 Mbar) have not yet been thoroughly studied. In the late 1960s, explosively driven shock experiments [11] only reached the highest pressure of ~2 Mbar. Laser-shock experiments in 1990s, with relatively high pressures but still below ~10 Mbar-revealed unusual behaviors of silicon under stress. For example, laser-shocked silicon exhibited a persistent thermal nonequilibrium between electrons and ions in the shock front even at ~6 Mbar pressures [12][13][14], while single crystals of silicon compressed by uniaxial shocks did not show the 3 normal hydrostatic-like compression [15]. Particularly, the reduction of lattice spacing occurs only along the shock-propagation direction. These abnormalities have called for more studies on the behavior of Si under the extreme pressures of P > 10 Mbar.
On the theoretical side, several classical simulations [16][17][18][19] have been devoted to study the shock-wave propagation in Si. Construction of a thermodynamically complete multi-phase equation of state, in order to accurately capture the many solid-solid and solid-liquid phase transitions exhibited by silicon, provides many challenges. Classical molecular dynamics simulations [18] indicated significant deviation in shock speed from experimental measurements for P < 1 Mbar. Using density functional theory (DFT), Swift et al. [20] investigated the equation of state (EOS) of Si through the various solidsolid phase changes and then into the liquid up to 0.5 Mbar. More recently, Strickson and Artacho [21] employed a DFT/constrained molecular dynamics approach (using a "Hugoniostat") and also obtained the Si Hugoniot up to 0.7 Mbar in good agreement with the results of Ref. [20]. Additional theoretical studies of the dynamic properties of shockcompressed Si have also been performed in the pressure regime of P < 10 Mbar [13,22,23]. Most recently, Militzer and Driver [24] have extended their path integral Monte Carlo (PIMC) method to study the EOS of hot dense Si over the P > 100 Mbar pressure range using novel Hartree-Fock nodes. There exists another first-principles technique for hot dense matter simulations-the orbital-free molecular-dynamics (OFMD) method [25]-based on the original idea of density functional theory [26]. The OFMD method has been extensively used to simulate a variety of hot dense materials from hydrogen/deuterium [27], helium-iron mixture [28], polystyrene [29,30] to plutonium [31]; however, it has never been tested against PIMC calculations for heavy elements.

4
In this paper, both the OFMD and PIMC calculations of the shock Hugoniot of silicon under extreme pressures (P > 100 Mbar) have been conducted. First we determined that the two first-principles calculations, PIMC and OFMD, agree well for hot dense material simulations. Next the PIMC/OFMD results were combined with the orbital-based Kohn-Sham DFT molecular-dynamics (KSMD) calculations to derive a global Hugoniot curve of Si in the fluid and plasma phases. Finally, we benchmark these first-principles calculations with the existing experimental data and EOS models to identify the shock behavior in Si. Our results from both KSMD + OFMD and KSMD + PIMC calculations indicate that silicon is much softer than the prediction from the extensively used SESAME and quotidian equation-of-state (QEOS) models. The stiff behavior predicted by these EOS models was caused by the overestimation of both pressure and internal energy. The QMD-predicted softening of silicon under extreme pressure might have implications in geophysics simulations and HEDP/ICF applications.
We hope these results will facilitate future high-pressure shock experiments of silicon.
The paper is organized as follows: a brief description of the three first-principles methods is given in Sec. II. Next, the derived Hugoniot curve of Si from the combined OFMD/PIMC and KSMD calculations is presented. We illustrate the softening of Si under high pressures by comparing with the widely used SESAME-EOS and QEOS models, as well as the existing experimental data in Sec. III. The heat capacity of shocked Si along its principal Hugoniot is also discussed in this section. Finally, the conclusions are presented in Sec. IV.

5
First-principles methods, such as DFT-based quantum molecular dynamics (QMD) [32][33][34][35], ground-state quantum Monte Carlo (QMC) [36,37], and finitetemperature PIMC [38], were developed to calculate the properties of materials under extreme conditions. The QMD method has two different implementations: (1) the orbitalbased Kohn-Sham formalism [39] with finite-temperature density functional theory [40] in conjunction with the molecular-dynamics method for ion motion (denoted here as "KSMD"); and (2) OFMD [25], which is based on the original DFT idea that the free energy of a many-electron system can be written as a function depending solely on the electron density. For most cases, the KSMD method has been proven to be an accurate and efficient way to calculate material properties under high compression at temperatures generally below the electron Fermi temperature T F . It becomes impractical for hightemperature (T > T F ) simulations because of the rapid increase in the number of orbitals required for convergence. The OFMD method is a natural extension of the KSMD method for high-T materials simulations. Using the density-matrix description of a manybody quantum system, the PIMC method is also currently applicable to simulate hot dense materials of elements up to Ar. Since the details of these first-principles methods have been documented elsewhere, only a brief description is presented here.

A. KSMD
The KSMD method has been successfully implemented in the Vienna ab initio simulation package (VASP) [41][42][43], in which electrons are treated quantum mechanically with a plane-wave, finite-temperature DFT description. The electrons and ions of the material are in thermodynamic equilibrium with equal temperature (T e = T i ).

6
The interaction between valence electrons and their parent ion is represented by a projector-augmented-wave (PAW) pseudopotential with "frozen" core electrons. The electron exchange-correlation potential is described by the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) functional [44]. Under were included such that the electron population of the highest band was less than 10 -5 .
The time step varied from δt = 1.25 fs to δt = 0.2 fs, respectively, for the lowest to highest compressions; good convergence was obtained for these parameters.

B. OFMD
The OFMD method [25] originates directly from the Hohenberg-Kohn theorem [26], with the free energy of an electron-ion system at any ion configuration can be 7 written as a function of the electron density. In the OFMD, the electronic density is found from a finite-temperature density functional theory [25] treatment in the Thomas-Fermi-Dirac mode with the kinetic-entropic form of Perrot [45]. The electron-ion interaction is obtained from a regularization prescription [25] with a cutoff radius selected at less than 10% of the Wigner-Seitz radius to avoid the overlap of the regularized spheres. Finally, the exchange-correlation takes a local-density Perdew-Zunger form [46].
At each time step of OFMD simulations, the electron free energy for an ionic configuration is first minimized in terms of local electron density. Then the classical ions are moved by the combined electronic and ionic forces as they were in the KSMD procedure. In the OFMD simulations of Si to determine its Hugoniot, 128 atoms in a cubic cell with a periodic boundary condition were used. The small time step, varying from δt = 2.4 × 10 -2 fs to δt = 4.8 × 10 -3 fs, is determined by the Si density and temperature considered. Finally, the thermodynamic quantities are averaged over the MD propagation of the system (5,000 to 10,000 steps).

C. PIMC
Path integral Monte Carlo (PIMC) is another first-principles method for hot dense material simulations, which uses the density matrix to describe a quantum many-body system in thermodynamic equilibrium. The density matrix has a convolution property; i.e., the density matrix at T can be expressed as a convolution of density matrices at very high temperatures (M × T), where the correlation effects between particles are small and a very good approximation for the density matrix exists [47][48][49]. The low-T density matrix can be found by carrying out the multidimensional integration along the imaginary time 8 path, which is needed to recover the full correlation effects at low temperatures. The Monte Carlo method is an effective way to perform such multidimensional integrations.
The electron-exchange effect is naturally taken into account through the inclusion of path permutations because the fermionic character of electrons requires an antisymmetric density matrix. The inclusion of permutation space results in inefficient sampling for lower temperatures since the positive and negative contributions to the integration are nearly cancelled out. This fermion sign problem [50] was overcome by restricting the paths using free-particle nodes or variational density-matrix nodes [51], which has succeeded in making PIMC simulations feasible to reach the lower temperatures of have produced accurate EOS tables of light elements such as hydrogen/deuterium [52][53][54][55], helium [56], carbon [57], water [58], neon [59], nitrogen [60], and oxygen [61].
For heavy elements such as Si, the challenge for PIMC simulations was to incorporate the effects of atomic bound states into the nodal structure. This was recently accomplished by constructing a thermal density matrix from Hartree-Fock (HF) orbitals [24].

III. RESULTS AND DISCUSSIONS
The KSMD method is employed for pressures less than ~100 Mbar, whereas both OFMD and PIMC methods are used for high-temperature points with P > 100 Mbar. Before comparing the high-pressure Hugoniot derived from the two calculations, two tests are performed: (1) a single Si atom in a cubic box of L = 5 bohr, and (2) EOS calculations for Si of different densities at T = 16,167,663 K. Figure 1 compares the predicted pressures and energies as a function of temperature from both PIMC and OFMD calculations for the first case. The temperature varies from T = 5 × 10 5 K to T = 4 × 10 7 K. The pressure comparison in Fig. 1(a) indicates that both methods give perfect agreement for all temperatures except the lowest temperature point at T = 5 × 10 5 Κ. At this low-T point, the PIMC simulations predict ~1% higher pressure than the OFMD result. Figure 1(b) shows that the predicted internal energies also compare well with each other except for the temperatures ~T = 10 6 to 10 7 K, for which the two results differ by a maximum of ~6% to 10%. It is noted that because of the different potential energy references used in these methods, a constant energy shift (+92.715 Ha/atom) has been applied to all OFMD-calculated energies. The zero energy was set to the completely ionized state, which is standard in quantum-chemistry calculations. Figure 1(c) shows the relative difference of internal energy calculated by KSMD, OFMD, and PIMC. It indicates that the OFMD calculation predicts ~10% higher energy than PIMC at T ~ 4 × 10 6 Κ, while the OFMD energies are ~6% lower than the PIMC ones around T = 10 7 K. This discrepancy may originate from the fact that the PIMC method handles the 1s core-electron ionization explicitly through the HF-nodal implementation while the OFMD method treats atoms without shell structures. Therefore, these two treatments of electron ionization can have slightly different features in the predicted Hugoniot curve at temperatures of ~T = 10 7 K (discussed below). Overall, the two methods provide reliable predictions of pressure and energy for such a simple situation.
Results from the second test case (using many particles) are presented in Fig. 2, where the quantity of PV and energy are plotted as functions of Si density for T = 16,167,663 K. Again, the pressures predicted by both OFMD and PIMC agree well overall within ~1%. The internal energies approach each other at the low density of ~0.1 g/cm 3 , but their difference increases gradually at high densities [see Fig. 2 Nevertheless, the predicted internal energy by OFMD and PIMC are still within ~3% at this temperature even for the highest explored density of ρ = 46.58 g/cm 3 . The Debye-Hückel model (the green dashed line in Fig. 2), gives excellent agreement with both results for low densities at this high temperature but deviates from the two first-principles calculations at higher densities, especially in pressures. These test results give us more confidence in pursuing the Hugoniot comparison.
To derive the principal Hugoniot, KSMD simulations were carried out for pressures below ~100 Mbar and OFMD/PIMC calculations for extremely high pressures (>100 Mbar). The initial state chosen was solid silicon (ρ 0 = 2.329 g/cm 3 ) in its diamond phase at an ambient pressure (P 0 = 1 bar   Figure 3 shows some small differences in the pressure range of several Gbars: the PIMC simulations predict a second maximum compression resulting from the 1s core-electron ionization, whereas such a feature is absent in the OFMD calculations. Namely, OFMD predicts the ionization to occur gradually along the Hugoniot, while PIMC shows distinct increases in shock compression because of the shell structure of the silicon atom. These differences can be attributed to the internal energy differences seen in Fig. 1(c) around T = 0.6 × 10 7 K -3 × 10 7 K, where the 1s electron ionization occurs.
To further examine the physics behind the 1s-electron ionization induced maximum compressions, we have calculated the heat capacity C v of silicon along its principal Hugoniot. As C v is a measure of the energy change rate with respect to temperature at fixed volume, one can perform two EOS calculations for each Hugoniot density (ρ) point with two temperatures of T+ΔT and T-ΔT slightly different from the Hugoniot temperature (T). The resulting internal energy difference ΔE from the two calculations can be used to compute the heat capacity (C v ≈ ΔE/2ΔT) for the Hugoniot point (ρ,T). The resulting C v is plotted in Fig. 4 as a function of the Hugoniot pressure for both KSMD+PIMC and KSMD+OFMD calculations. As the 1s-electron ionization process acts like a "heat-sink" for the system, one expects the heat capacity should dramatically increase during the ionization of inner-shell electrons. This is exactly what we see in Fig. 4 that the KSMD+PIMC calculation (blue dash-dotted line) gives a second peak of C v at P ≈ 3 Gbar, corresponding to the second compression peak seen in Fig. 3.
After the 1s-electron ionization completes, the heat capacity approaches the ideal-gas limit (black dashed line), as fully ionized Si plasma is formed. Instead of giving double peaks, the KSMD+OFMD calculation predicts a broad single peak of C v . Again, this originates from the lack of shell structure in Thomas-Fermi like theories. Whether or not the ionization-induced second maximum compression can be seen in experiments remains a challenge for future extremely high pressure measurements.
Returning to Fig. 3, the first-principles results are compared with the widely used SESAME-EOS [63] and the QEOS models [64]. Both models are based on the chemical picture of matter, meaning that the total free energy can be decomposed into the cold curve, the ionic excitation, and the electron thermal excitation. For example, Sesame EOS models are typically constructed (constrained) by the best available experimental data (typically limited). Specifically, for Sesame 3810 (Si) constructed in 1997, the EOS below the solid-liquid phase transition was based on experimental Hugoniot data [11,65,66]. For conditions above the liquid phase transition, the EOS was constructed such that the shock Hugoniot was "similar" to Germanium (Sesame 3950) up to 4.4 Mbar. The ion thermal contribution is based on a Debye model with a correction for the liquid specific heat beyond the melt temperature [67]. The correction also ensures that in the high temperature limit the proper model (ideal gas) is recovered that gives a shock Hugoniot compression ratio ρ/ρ 0 =4. This comparison in Fig. 3 indicates that under shock compression silicon is much softer than predicted by the traditional chemical-picture understanding of materials. For a shock density of ρ = 5 to 9 g/cm 3 , the pressures predicted from our first-principles calculations are smaller than the corresponding SESAME and QEOS pressures by a factor of ~2 to 3. In other words, for the same shock pressure in the range of 5 to 100 Mbar, the first-principles-predicted shock density is ~20% higher than the SESAME and QEOS prediction. The maximum compression (ρ/ρ 0 ) 14 changes from the model-predicted value of ~4.0 (QEOS) and ~4.6 (SESAME) to ~5.0 in our first-principles calculations. Perhaps in the case of SESAME this amount of difference is not surprising based on the way that SESAME 3810 was constructed. For QEOS, it appears that the EOS is constructed so to follow the Hugoniot until the ideal gas limit (four fold compression), at which point the Hugoniot climbs nearly vertically in pressures. Looking at a typical density/temperature condition along the principal Hugoniot (e.g., ρ=6.5 g/cm 3 and T=62500 K), we find the SESAME 3810 predicted pressure and internal energy (P≈11.45 Mbar and E≈34.5 eV/atom) are about ~45% higher than our first-principles calculations (P≈7.76 Mbar and E≈24.3 eV/atom). The overestimation of pressures and energies in SESAME 3810 is the cause of the stiff behavior we saw in Fig.3. Finally, the earlier experimental data [11] and other Hugoniot measurements in 1970s and 1980s in the liquid phase [65,66] are put into the same figure for comparison, respectively represented by the open triangles, diamonds, and squares in Fig. 3. These Hugoniot data were obtained from explosively driven shock experiments.
To the best of our knowledge, no published data exists for Hugoniot measurements in pressures above 10 Mbar. The opaqueness of Si for most VISAR (velocity interferometer system for any reflector) laser wavelengths [68] is one of the hurdles for accurate shock measurements in silicon. Since some of these explosively driven shock data at lowpressures (P<1 Mbar) were used to constrain the construction of SESAME 3810, the agreement there was guaranteed. The recent KSMD data by Strickson and Artacho [21] also agree with the data at these low pressures. However, for pressures between 1 and 2 Mbar our KSMD calculations are in better agreement with the available high-pressure data [11]; while the SESAME-EOS, guided by the Ge-based model, results in twice the 15 pressure compared to the KSMD result. It was noted that at the measured highest shock density of ρ = 4.6 g/cm 3 , the SESAME and QEOS predicted pressures are about twice higher than both our KSMD results and the experimental value of P  1 to 2 Mbar. It is expected that future experimental data in the 10-to 100-Mbar pressure regime will unambiguously benchmark the predicted softening of Si.
To further understand the shock-induced structure change, we turn to the pair correlation function along the predicted principal Hugoniot of Si (Fig. 5). A strong peak at the Si-Si inter-particle distance of ~2.25 bohr exists for compression at T = 5000 K and ρ = 4.38 g/cm 3 . As the shock compression increases to ρ = 5.50 g/cm 3 and T = 31,250 K, this peak moves to a smaller ionic distance, and the amplitude decreases.
Increasing the shock strength further, this peak becomes much less pronounced at ρ = 7.32 g/cm 3 and T = 125,000 K, which manifests the status of the structureless fluid phase.
Finally, the first-principles predicted shock temperatures are compared with the SESAME and QEOS models in Fig. 6. For shock pressures exceeding ~3 Mbar, the first-principlespredicted shock temperatures are ~30% higher than the SESAME EOS values and much higher (~50% to 200%) than the QEOS predictions. At a fixed pressure, the predicted shock from first-principles calculations, reaching both higher densities and temperatures, can have implications in HEDP and ICF applications. This is because the shock velocity from first-principles calculations can be significantly lower than the widely used EOS models, while the FP-predicted particle velocity is higher than SESAME and QEOS models.

IV. SUMMARY
Two first-principles methods of PIMC and OFMD have been used to calculate EOS of silicon at extremely high pressures. Combining these high-temperature calculations with the KSMD simulation of low-temperature EOS, the principal Hugoniot curve of silicon for pressures varying from ~1 Mbar to above 10 Gbar was derived.
Overall, the two results agree very well with each other, although a small difference exists at pressures around several Gbar. In this high-pressure regime, the 1s core-electron ionization of Si can cause a maximum compression, since the ionization process acts like a heat sink. This is evidenced by the peaked heat capacity in this pressure range of several Gbar, in which C v is larger than the ideal-gas limit of fully-ionized Si plasma.
When the first-principles predicted Hugoniot of Si is compared with the widely used SESAME and QEOS models, we find that silicon under pressure is ~20% or more softer than what was generally believed in the chemical picture of matter. In future work, the broad range of first-principles data will be incorporated into an improved Sesame model,  calculations compared to the two different models (SESAME-EOS and QEOS) and the existing experimental data by Pavlovskii [11], by Gust and Royce [65], and by Goto et al [66]. Also, the recent KSMD results by Strickson and Artacho [21] at lower pressures are also shown. The diamond-phase solid silicon was chosen as the initial state having a density of ρ 0 = 2.329 g/cm 3 . calculations compared to the two different models (SESAME-EOS and QEOS) and the existing experimental data by Pavlovskii [11], by Gust and Royce [65], and by Goto et al [66]. Also, the recent KSMD results by Strickson and Artacho [21] at lower pressures are also shown. The diamond-phase solid silicon was chosen as the initial state having a density of ρ 0 = 2.329 g/cm 3 . Hugoniot between our first-principles calculations and the two different SESAME-EOS model (SESAME and QEOS) predictions.