{'fraction_non_alphanumeric': 0.03953127760221407, 'fraction_numerical': 0.016899915598563214, 'mean_word_length': 4.362947368421053, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A multiphase equation of state for carbon addressing high pressures and temperatures

We present a 5-phase equation of state (EOS) for elemental carbon. The phases considered are: diamond, BC8, simple-cubic, simple-hexagonal, and the liquid/plasma state. The solid phase free energies are constrained by density functional theory (DFT) calculations. Vibrational contributions to the free energy of each solid phase are treated within the quasiharmonic framework. The liquid free energy model is constrained by fitting to a combination of DFT molecular dynamics performed over the range 10,000 K100,000 K (both for {\rho} between 3 and 12 g/cc, with select higher-{\rho} DFT calculations as well). The liquid free energy model includes an atom-in-jellium approach to account for the effects of ionization due to temperature and pressure in the plasma state, and an ion-thermal model which includes the approach to the ideal gas limit. The precise manner in which the ideal gas limit is reached is greatly constrained by both the highest temperature DFT data and the path integral data, forcing us to discard an ion-thermal model we had used previously in favor of a new one. Predictions are made for the principal Hugoniot and the room-temperature isotherm, and comparisons are made to recent experimental results.


I. INTRODUCTION
The high pressure and temperature equation of state (EOS) and phase diagram of carbon have received considerable attention of late. Recent laser-shock and ramp-compression studies [1][2][3][4][5][6][7], together with shock measurements performed with magnetically-driven flyer plates [8], have produced data in the range from P = 0 − 50 Mbar. In addition, theoretical work on the EOS and phase diagram in this same range have yielded predictions which are largely (though not completely) in accord with these experimental data [8,9]. These studies were conducted with density functional theory (DFT) molecular dynamics (MD). Much of this recent focus on the carbon EOS, specifically in states of compression reached when starting in the diamond phase, has arisen from the interest of using high density carbon as an ablator material for capsules designed to achieve fusion at the National Ignition Facility [10]. While these experimental and theoretical studies have been useful in constraining the EOS of carbon for this and related applications, it is crucial to note that the states reached by the ablator in inertial confinement fusion (ICF) are expected to include temperatures in excess of tens of eV [10]. Such conditions have not yet been sufficiently characterized.
At such high temperatures, DFT-MD is extremely challenging to perform, due to the large number of high-lying, partially-occupied single-electron states that must be included for an accurate rendering of the electronic thermal excitations. Indeed, an orbital-free variant of DFT-MD has been developed and used to handle such high temperature applications [11], but as yet, a satisfactory treatment of atomic shell structure in this approach is lacking.
The approach as implemented thus far [11] is (somewhat) reminiscent of a Thomas-Fermi treatment for the electrons [12]. Another altogether different approach to describing the electronic structure of the high-T plasma state is path integral Monte Carlo (PIMC) [13]. In this approach, unlike in typical implementations of DFT, there is no mean-field assumption made for the many-electron problem, and the imaginary time treatment makes high-T simulations more efficient to perform than low-T simulations. Although assumptions regarding the nodal surface of the many electron density matrix necessarily introduce approximations for the treatment of atomic shells, it is very encouraging that recent work on the C plasma has demonstrated that EOS predictions can be made with this method which smoothly interpolate between DFT-MD results at the lower temperatures and the high-T Debye-Hückel limit [14].
The ultra-high pressure regime near T = 0 has been studied recently as well, using ab initio electronic structure techniques of the DFT variety together with a random search method to find the thermodynamically stable crystalline phases at pressures between 10 -1000 Mbar [15]. This work extended earlier quantitative studies of the C phase diagram that treated the diamond and BC8 phases up to P ∼ 20 Mbar [16]. The random search established the sequence of stable phases to be diamond → BC8 → (slight modification of) simple cubic (sc) → simple hexagonal (sh) → face-centered cubic (fcc) [15]. Though these authors did not attempt to predict the details of the C melt curve in these extreme conditions, their analysis using a quasiharmonic description of the phonons for these solid phases suggested the existence of a BC8-sc-sh triple point, as seen in Fig. 3 of their manuscript [15].
A few years ago, Correa et al. used the available theoretical understanding of the C phase diagram at the time [16] to produce a three-phase EOS model for C focusing specifically on the regime of P = 0 − 20 Mbar and T = 0 -20,000 K [9]. The phases considered were: diamond, BC8, and the liquid. The individual phase free energy models were fit to the results of DFT-MD calculations. Since the final multiphase EOS model was designed to be used in hydrocode simulations spanning a wider range of conditions, these researchers embedded their detailed three-phase EOS model into a more coarse-grained model [17] that respected the ideal gas limits at high-T and low-ρ and the Thomas-Fermi limit [12] at high-ρ.
While satisfactory in a broad sense, this older EOS model suffers from three main drawbacks: 1. The embedding of the DFT-based 3-phase model into the coarse-grained model created kinks in the thermodynamic functions, even though an attempt was made to smooth out such features by interpolation. 2. There are no solid phases beyond BC8; sc and sh phases, for instance, are not included. 3. The electronic excitations of the high-T liquid are treated with an average-atom Thomas-Fermi model, so atomic shell structure is not included in sufficient detail.
In this work, we remedy these deficiencies by constructing a 5-phase EOS for C, again based on ab initio calculations (for densities between ∼ 1 -25 g/cc), in which the free energy models for each phase are defined over wide ranges of ρ and T and exhibit sensible limiting behavior. No embedding into a simpler model [17] is performed. The two extra phases included are the sc and sh phases. While our phase diagram differs slightly from that of Ref. [15], the inclusion of these extra phases allows our EOS to be in accord with ab initio predictions for pressures up to well over 100 Mbar . The electronic excitations   in the liquid are treated with an average-atom model called PURGATORIO which, unlike Thomas-Fermi, includes the treatment of atomic shell structure within a Kohn-Sham DFT framework [18]. To constrain the ion-thermal model for the liquid, we perform DFT-MD in C up to 100,000 K and PIMC caculations [14] from ∼ 200,000 K (depending on the density) to over 10 8 K. We establish that the best fit to these EOS predictions requires a model in which the decay of the ion-thermal specific heat is much faster than previously expected.
Our newly-developed "cell model" [19], in conjunction with the PURGATORIO electronthermal term, provides an excellent fit throughout the range of both the DFT-MD and the PIMC data. In what follows, we describe the details of our DFT and PIMC calculations, present the models we use for the free energies of the individual phases, and compare various thermodynamic tracks through our EOS with the results of recent experiments performed on high-energy laser platforms [1,6,7].

II. THEORY AND SIMULATION
Our computational results fall into two categories: 1. Calculations of the EOS, by which we mean internal energy, E, and pressure, P , as functions of density, ρ, and temperature, T , and 2. Calculations of intermediate quantities we which use to build the EOS. These include cold curves (E and P as functions of ρ, for ions fixed in position [20]), phonon densities of states (PDOS), and electronic excitation contributions to the free energy. The first category of quantities we extract from DFT-MD (for solid and liquid phases) and PIMC (liquid phase). The second group of quantities we extract solely from DFT calculations for the solid phases. It may seem at first as if our multitude of results overdetermines the EOS, particularly in the liquid phase where both DFT-MD and PIMC results reside. However we find it crucial to make use of this full suite of results; the approximations inherent in our DFT and PIMC calculations are altogether different, and it is important to check whether certain features of our results are artifacts of the simulation method or are robust indicators of the true EOS of C. This is especially the case for the liquid at high temperatures, where we find (with both methods) results which contradict some of our previous assumptions (see below). Fig.1 shows the locations of our DFT-MD and PIMC simulation points in the (P, T )plane. Also shown for reference are our predicted C phase lines, presented and discussed in  approaches to obtaining these data.

A. DFT calculations
With DFT methods, we perform calculations of cold curves, phonon densities of states, electronic excitation contributions to the free energy, and the EOS itself (E and P as functions of ρ and T ). We also present a limited number of calculations of ionic diffusivity, from which we can infer melt behavior. For all of these, we use the VASP code [21], together with projector augmented wave (PAW) pseudopotentials [22]. We use a "hard" PAW with a core radius of 1.1 bohr and 4 valence electrons and the generalized gradient approximation (GGA) of DFT with the Perdew-Burke-Ernzerhof (PBE) exchange correlation functional [23]. The plane-wave cutoff is set to 1300 eV.
For the phonon calculations, we use the primitive cell with a MP k-point grid of 40x40x40 for the simple cubic (sc), the simple hexagonal (sh) phases and for fcc. The sh is done with a c/a ratio of 0.986 evaluated by performing cell optimization for several compressions. No appreciable variation in the c/a ratio with volume was observed. For the diamond and bc8 phases we used a 8 atom unit cell with a 20x20x20 k-points grid. For the sc phase, we compared with a 2 atom unit cell with a pmma spacegroup (tetragonal distorsion) [15] which was calculated with a 20x40x40 k-point grid.
For the molecular dynamics simulations, we use 64 atoms and perform Born-Oppenheimer MD (BOMD) within the NVT-ensemble with a Nosé-Hoover thermostat [24]. We use a time step of 0.75 fs in order to converge the internal energy and pressure to the desired accuracy.
The electronic density is constructed from single-particle wave functions by sampling at the ( 1 4 , 1 4 , 1 4 ) point of the Brillouin zone. The electron occupation numbers are taken to be a Fermi-Dirac distribution set at the average temperature of the ions. For the different densities and temperatures, we use a sufficient number of bands such that we have at least 40 bands with occupation numbers smaller than 0.00001. The MD is run for 10000 steps with the last 5000 steps used for averaging the internal energy and the pressure. respectively; see below). At higher temperatures still, the thermal excitations of the 1s electrons can no longer be neglected. This can be addressed by having all 6 electrons in the valence (for the current DFT-MD calculations, however, we find this not to be necessary for temperatures of 100,000 K or lower). At the highest compressions, the perturbation of the 1s orbital may be significant and can again be addressed by having all 6 electrons in the valence.
In order to assess the sensitivity of our results to the PAW potential chosen, we construct PAW potentials with different cutoff radii, which range from 1.1 to 0.8 bohr radii. The Vanderbilt projectors generation scheme is utilized [25]. As the cutoff radii decrease, in order to maintain and optimize the PAW performance, we add additional partial waves for both s and p angular momenta. Different sets of reference energies are chosen, which are determined not to affect our computational results. Moreover, when the carbon-carbon distance is sufficiently small, the 1s core state is also included in the PAW. For all cases tested, we find that a 5442 eV plane-wave cutoff energy is enough to reach convergence.
For the PAW test, we use fcc cells with lattice constants chosen so that the carbon-carbon nearest neighbor distances are representative of the distances observed during the highest-P cold curve and MD runs. These tests are performed with the Quantum-Espresso (QE) package [26]. Figure 3 shows a comparison of the pressure obtained for this fcc cold curve with the different PAWs (for both VASP and QE). We only observe appreciable deviations between the different PAWs for densities above 60 g/cc, corresponding to carbon-carbon distances of ∼ 1.3 bohr (0.7Å). The difference in the calculated pressure between the VASP R cut =1.1 bohr PAW potential and the 6 electron PAW is only 1% up to densities of 100 g/cc. This corresponds to carbon-carbon distances of 1.1 bohr which is essentially the smallest carbon-carbon distances sampled during the highest temperature and highest density MD runs. However, this is no guarantee that the dynamics are not affected by the overlapping PAW spheres. In this sense, we submit that our DFT-MD EOS data for the very highest densities and temperatures reported here may be of lower accuracy than those for the more moderate conditions. Nevertheless, we stress that the bulk of our conclusions below are unaffected by this fact. PIMC is the most accurate and efficient first-principles simulation technique to study the equilibrium properties of quantum systems in high temperature plasma states. Fermionic PIMC simulations have been applied to study hydrogen [27][28][29][30][31][32][33][34][35], helium [36,37], hydrogenhelium mixtures [38], one component plasmas [39,40] and most recently to simulate carbon and water plasmas [14]. In PIMC simulations, electrons and nuclei are treated equally as Feynman paths in a stochastic framework for solving the full, finite-temperature, quantum many-body problem. The natural operator to work with in this context is the thermal density matrix represented as a path integral in real space, where H is a many-body Hamiltonian, β = 1/k B T assumes the role of imaginary time, R = (r 1 , r 2 , · · · r N ), r i is the position of the ith particle, and ∆τ = β/M is the step size in imaginary time intervals. S[R(τ )] is the action, which determines the weight of every path.
The kinetic energy operator controls the diffusion of the paths in imaginary time and keeps bead R i at adjacent time slices close together. A thermodynamic function corresponding to an operator,Ô, can be derived from Ô = T r(ρÔ)/Z c , where Z c = T r(ρ) is the canonical partition function.
PIMC explicitly addresses all the physics of high-T plasmas including effects of bonding, ionization, exchange-correlation, and quantum degeneracy on an equal footing, and thereby circumvents both the need to occupy single-particle states, and the need to employ exchangecorrelation approximations inherent in DFT. The Coulomb interaction between electrons and nuclei is introduced using pair density matrices derived using the eigenstates of the two-body Coulomb problem [41]. The periodic images are treated using an optimized Ewald break-up [27,42] applied to the pair action [13].
Challenges with the PIMC formalism arise in fermionic simulations because only antisymmetric eigenstates contribute to the partition function. Those can be projected out by introducing an additional sum over all permutations, ρ F (R, R , β) = 1/N ! P (−1) P ρ(R, PR , β). Straightforward integration methods all lead to unstable algorithms due to the fermion sign problem, which is the result of the near complete cancellation of positive and negative contributions to the many-body density matrix.
A fixed-node approximation [43,44] was initially introduced to solve the sign problem for ground-state quantum Monte Carlo calculations. For the finite-temperature PIMC simulation, a restricted-path method was developed [45], where one only integrates over all paths that do not cross the nodes of a trial density matrix, ρ T (R, R ; β), which must be available in analytic form. The fermionic version of Eq. 1 then reads In addition to weighting all paths according to their action and summing over permutations, one must also check whether the sign of the trial density matrix has remained positive during every Monte Carlo move. If a sign change occurs anywhere along the new path, the proposed move is rejected. The R argument in the nodal restriction plays a special role because the sign of the trial density matrix at all other time slices depends on it.
The restricted-path method provides the exact answer if the nodes from the exact manybody density matrix are used. Since exact nodes are not known for interacting systems, one must work with approximate nodes. However, within any given set of nodes, the restricted path method will obtain the best possible solution that includes all interaction and correlation effects. The PIMC method used here employs a free-particle nodal structure, which has been shown to be sufficient for carbon at temperatures where atoms still have occupied 1s states, but partially occupied 2s states (T ≥ 2.5 × 10 5 K) [14].
Using the restricted-path, free-particle nodal framework, we use all-electron PIMC to compute equations of state for carbon at densities of 0.1, 3.18, 8.5, and 11.18 g cm −3 over a temperature range of 10 5 -10 9 K (though we present data pertaining to all four isochores in the Supplementary Material, we stress that 0.1 g/cc is below the low-density limit of applicability of the EOS model described below) [46]. A sufficiently small time step associated with PIMC path slices is determined by converging total energy as a function of time step until it changes by less than 0.2%. We use a time step of 0.0078125 Ha −1 for temperatures below 5×10 6 K and, for higher temperatures, the time step decreases as 1/T while keeping at least four time slices in the path integral. In order to minimize finite size errors, the total energy was converged to better than 0.2% for a 24 atom cubic cell.

III. CONSTRUCTION AND DISCUSSION OF EOS MODEL
We make the fundamental assumption that the Helmholtz free enegy, F = E − T S, of each phase can be decomposed into the following terms: where F ion and F elec represent the contributions from ionic and electronic excitations, respectively. While the identification of a temperature-independent piece (E cold ) is always possible by virtue of the fact that it is merely definitional, the separation of the thermal part into decoupled ionic and electronic pieces is a major assumption. In a practical sense, this assumption is at least a reasonable starting point since for solids at lower-T , it has been shown to be accurate [9,47], and for the liquid at very high temperature (ideal gas limit) it is also trivially satisfied because interactions between particles are of no importance. It is not necessarily justified for intermediate temperatures, but we shall invoke the approximation at any rate, as has been done in the past [9,19,48]and take comfort in the fact that the ab initio EOS data to which we fit in the liquid (obtained from DFT-MD and PIMC) does not make this assumption.

A. solid phases
For the solid phases we consider diamond, BC8, simple cubic (sc), and simple hexagonal (sh), we take E cold (V ) to be the internal energy of the perfect crystalline lattice (with motionless ions) at T = 0 [20] with the electrons in their ground state. We fit this term to the DFT data for the internal energy of the perfect lattice using the Vinet equation of state [49], modulus at V 0 , and φ 0 is the internal energy at V 0 . Fig. 4 shows the computed cold curves of the various phases, together with our fits to them, using Eq.4. The cold curve parameters for each phase (φ 0 , B 0 , B , and V 0 ) are given in Table 1, along with the other phase-dependent parameters we discuss below. Our (PBE) GGA-DFT results for the solid cold curves are very similar to those of Ref. [15], but there are differences which result presumably from the use of different pseudopotentials and associated choices for the plane wave energy cutoff in the determination of the internal energy. Though small, these differences give rise to rather pronounced changes in the predicted phase transition pressures (BC8 → sc → sh) when compared to the phase diagram presented in Ref. [15]. This will be discussed more below, after the thermal components of the free energies are considered and our phase diagram is presented. At this point, we note that we have checked that our inferred cold transition pressures, determined from the two-phase Maxwell construction [47], should be inaccurate by no more than a few percent due to our use of the Vinet fitting form, Eq.4 (though we in no way imply that PBE GGA-DFT is necessarily this accurate for predicting these transition pressures at such high compressions).
For F ion (V, T ), we use the ion excitation free energies of the various lattices computed for cases in which the electrons are in their instantaneous ground state (for each ionic configuration). In other words, we invoke the Born-Oppenheimer approximation. In addition, the quasiharmonic approximation is assumed for each phase. We investigated anharmonicity for diamond and BC8 phases in Ref. [9] and showed its effects to be small; the effects of anharmonicity in sc and sh phases await further study. Given this quasiharmonic description, we use double-Debye models [9,19] for the diamond, BC8, and sc phases, and a regular single-Debye model for the sh phase. All are fit to the V -dependent phonon densities of states (PDOSs) for these various phases computed with DFT using a linear response approach. In the double-Debye prescription, the PDOS for a given volume is presumed to consist of two overlapping features, each of which is of the form: aε 2 for ε ≤ θ and 0 for ε > θ, where the two unequal θ(V ) functions are related to various moments of the PDOS [9]. Since the double-Debye model allows for two different values of θ (say, θ A and θ B ) for each V , it is useful in cases where the high-T phonon moment, k B θ 0 (V ) = e 1/3 exp ∞ 0 dε ln(ε)D V (ε) , is rather different from the low-T zero-point moment, k B θ 1 (V ) = 4/3 ∞ 0 dεεD V (ε). This is certainly the case for the diamond and BC8 phases, especially at high compressions, as discussed in detail in Ref. [9]. Our calculations for the sh phase shows that θ 0 and θ 1 are very nearly equal throughout the predicted range of its stability. Thus, we use a single-Debye free energy for the sh ion-thermal term to reduce the total number of fitting parameters. Fig. 5 shows the three phonon moments, θ 0 , θ 1 , and θ 2 (k B θ 2 (V ) = 5 3 ε 2 D V (ε)dε), as functions of V for BC8, sc, and sh phases as determined from our computed PDOSs. Note that the θ 0 (V ) for the different phases nearly coincide in the neighborhood of the volumes where their cold curves cross. For this reason, the T vs. P phase lines are rather vertical in nature (see below). It is from these moments that we derive (for the diamond, BC8, and sc phases where the double-Debye model is used) the two Debye temperatures, θ A and θ B , by solving the system of equations given as Eqs.9 -12 of Ref. [9]. For all θ A,B (V ) functions for the various phases, we assume the V -dependence: where θ 0 is the value at a reference volume V p . This arises from the assumption that the ion-thermal Grüneisen parameter for θ(V ) is equal to γ = AV + B [9,48]. The relevant ion-thermal parameters (θ 0 , A, B, and V p ) which determine the various Debye temperatures for each phase are presented in Table 1. In terms of θ A , θ B , and θ 1 , the double-Debye free energy is written as where F A and F B are the single-Debye free energies for the individual A and B peaks associated with a given PDOS.
The electronic excitation term for each solid phase, F elec (V, T ), is taken be of the form, which is motivated by a performing a Sommerfeld expansion of the electronic excitation free energy, assuming T /T Fermi to be a small parameter [47]. This contribution is extracted from PBE GGA-DFT calculations in which electronic occupancies are constrained by a Fermi-Dirac distribution with temperature T , but for ions fixed in their lattice positions.
As discussed in Ref. [9], this term is essentially zero for the diamond phase throughout its stability field, and is exceedingly small (though non-zero) for the BC8 phase as well [50]. For this reason, electronic excitation contributions were intentionally neglected by those authors in their solid-phase EOS models. Here we choose to include these terms; their values for the sc and sh phases are somewhat larger than those of BC8, because the electronic density of states near E Fermi generally increases with compression. Still, their inclusion has a rather small effect on EOS and the resulting phase lines. We fit these small DFT-based terms by assuming α(V ) = α 0 (V /V e ) κ [9,47,48]; our choices for α 0 , κ, and V e are presented in Table 1 for each phase. Note that there are nonzero values for α 0 which appear for the diamond and BC8 phases in this table; these result from earlier extractions of anharmonic ion-thermal free energy contributions [9] (which happen also to have a T 2 dependence), and are not due to electronic excitations [51].
In order to check the validity of our solid-phase EOS model for C, we compare to our DFT-MD results for pressure and internal energy in the diamond phase. These results exist throughout a range where diamond is predicted to be the stable phase (3.0Å 3 /atom < V < This is likely to be a very small contributor to the EOS differences, however, for the relevant moments of the PDOSs at various volumes are strikingly similar to those resulting from the computational framework adopted in this work.

B. liquid phase
For liquid C, we assume a model much like that of our solid phases, consisting of separate cold, ion-thermal, and electron-thermal terms as shown in Eq.3. Here, however, there is no natural way to determine them all independently, because the liquid cannot be described as a perturbation away from a fixed configuration of ions. Nevertheless, we assume that Eq.3 holds for the liquid and base our model on the liquid free energy assumed in Ref. [9], though with some important differences. In this approach, the cold piece is once again assumed to be of the Vinet form [49] (Eq.4), though augmented with a few so-called break-points which introduce bends in E cold (V ) over localized regions in V [9,52]; additional bends at small V are included to improve agreement with our liquid DFT-MD data performed at densities far higher than those considered in the work of Ref. [9]. The ion-thermal term at lower-T (up to ∼ twice T melt ) is treated in the scheme of Chisolm and Wallace [53], which is an effective Mie-Grüneisen model for the liquid with a V -dependent characteristic temperature, θ(V ), which we assume to have the form of Eq.5. This model was used before for liquid C [9], and is constrained to have an ion-thermal specific heat, C ion V = 3k B /atom; this is reasonable considering the (T < 20,000 K) DFT-MD data that was provided in order to inform the liquid C EOS model in that earlier work [9]. The contribution to the liquid C EOS from electronic excitations is predicted to be quite small for the lower temperatures, even though liquid C is predicted to be metallic [9]. Nevertheless, we use the average-atom-in-jellium DFT model, PURGATORIO [18], for the electron-thermal model of the liquid throughout the entire field of stability of the phase. As expected, its effects are pronounced at high-T ; we discuss them below. The comparison of the resulting liquid C EOS model to our DFT-MD  EOS data taken over a wide range of compression is shown below, after the high-T liquid is addressed.

The high-T liquid: approach to the ideal gas
Since we now have first-principles electronic structure data (DFT-MD and PIMC) for liquid C at far higher T than we had previously [9], we are able to test various assumptions regarding the manner in which C ion V evolves from its Dulong-Petit value (3k B per atom) to the ideal gas limit (3k B /2 per atom). To frame this discussion, we use two different high-T additions to the liquid ion-thermal free energy, the predictions of which we will compare to these data: 1. A variant of the the Cowan model [54,55], and 2. The cell model [19]. Both models essentially satisfy the limit: for sufficiently high T , where T ref (V ) is a V -dependent reference temperature. For the Cowan model, the exponent ν= 1/3; for the cell model, ν= 2. Each have C ion V = 3k B for T well below T M and C ion V = 3k B /2 as T → ∞, but the rate at which the high-T limit is approached is very different. In our version of the Cowan model, we choose the specific form of C ion V to be where θ(V ) = e −1/3 θ 0 (V ) is the characteristic temperature of the liquid model in Ref. [9].
This expression is used for T > T M ; for T ≤ T M , the Mie-Grüneisen form is used: In the cell model, the decay to the ideal gas form is affected by an additive (rather than piecewise) correction, derived by appealing to the classical approximation to the partition function for an ensemble of Einstein oscillators with each particle strictly confined [19] to a radius R. The resulting free energy, giving rise to C ion V approximately (for T T M ) equal to that of Eq.8 with ν = 2, is where here, k B T M (V ) = mk 2 B θ 2 R 2 /2 2 [56]. In the simplest prescription, the confining radius is related to the volume per atom by 4πR 3 /3 = V . This gives rise to a slightly incorrect value for the entropy at T = ∞, due to the explicit neglect of permutations between particles which are confined to individual cells. The correct high-T limit of the entropy can be recovered by instead choosing the relationship, 4πR 3 /3 = V e (where ln(e) = 1) [19]. However, we find below that for fitting to E and P within the range of (ρ, T ) considered here, the choice of 4πR 3 /3 = V is preferred. In our practical implementation of the cell model we replace the first term in Eq.11 (Mie-Grüneisen) by a Debye model free energy with the identical characteristic temperature, θ, to capture the thermodynamics of quantum ions at the lower temperatures (however subtle that distinction may be for liquid C). Fig. 7 shows internal energy, E, versus T for liquid C on a densely-spaced grid of isochores for densities: 3 g/cc < ρ < 12 g/cc. Red points are the results of DFT-MD extending up to T = 100,000 K as described in Section II A. Black dashed-dotted lines are the predictions of a liquid C EOS model for these same (ρ, T ) states, which has been constructed exactly along the lines of Ref. [9], but with two important additions: 1. F ion includes the Cowan model correction as described above (Eq.10) and in Ref. [55]; 2. The PURGATORIO model for C is used for F elec . Note that while the agreement is good for T < 20,000 K, as ensured by the fitting scheme of Ref. [9], the slope of E vs. T for higher T is much larger in the model, indicating that the model's C V is too large above 20,000 K. Also shown in Fig. 7 is the results of an otherwise identical model (blue lines), but with F ion chosen to be the cell model; once again, PURGATORIO is used for F elec . Red points and blue lines in this plot are in far better accord. Assuming the electronic excitation contribution to be that of the PURGATORIO model, this comparison is showing that a faster decay of C ion V is greatly favored over the slower decay previously assumed [55]. This is at least vaguely troubling, since many EOS tables in the Livermore Equation of State (LEOS) database are constructed using the paradigm of QEOS [17], which assumes the Cowan model for F ion . Of course, our present comparison is just for carbon, and is also for a somewhat restricted range of ρ and T . Furthermore, it is entirely possible that 1. This conclusion depends on the electron-thermal model chosen, here PURGATORIO rather than the choice outlined in QEOS, average-atom Thomas-Fermi [12,17], and 2. The DFT-MD data to which we compare is somehow lacking in fidelity at the higher-T .
To explore both possibilities while investigating a far wider range of T , we compare to our PIMC results for E and P , as generated by the prescription outlined in Section II B. Fig. 8 shows the fractional differences, (X model − X PIMC )/X PIMC [58], between liquid C EOS model predictions and the results of PIMC for three isochores, ρ = 3.18 g/cc, 8.5 g/cc, and 11.18 g/cc. Both internal energy and pressure differences are displayed. For the liquid C EOS models, we consider four variants: 1. F ion = Cowan, F elec = PURGATORIO (green), 2. F ion = Cowan, F elec = Thomas-Fermi (magenta), 3. F ion = cell model, F elec = PURGATORIO (red), 4. F ion = cell model, and F elec = Thomas-Fermi (blue). Note the exceptionally wide range of T displayed in log-scale on the x-axis (5 × 10 5 K → 10 8 K). Each set of points (of a given color/symbol) includes data from three isochores with densities: ρ = 3.18 g/cc, 8.5 g/cc, and 11.18 g/cc. All model variants give results which are coincident with PIMC at sufficiently high-T . However, there are significant deviations in E and P for lower temperatures. Clearly, it is the red set of points, corresponding to F ion = cell model, F elec = PURGATORIO, which presents the smallest deviations from PIMC. Furthermore, while the preference for PURGATORIO over Thomas-Fermi when fixing the ion-thermal model is notable but modest, the preference for the cell model over the Cowan model is quite dramatic for both choices of F elec . Finally, it is interesting that there is a pronounced peak in these fractional differences at T ∼ 10 6 K for both E and P , and for all three isochores presented. It is possible that this arises from an ionization feature in the PIMC which is somehow manifested rather differently in the Thomas-Fermi and PURGATORIO models.
However, it is also curious that these deviations attain their maxima at roughly the same T even for very different densities, since one might expect such ionization features to be rather dependent on density. The detailed shape of these curves awaits further analysis.
The conclusion of these comparisons is that both our DFT-MD and our PIMC results indicate that the rate of decay of C ion V to its ideal gas value is far too slow in the Cowan model, while the faster decay of the cell model provides far better agreement for liquid C in this range of densities (3 g/cc < ρ < 12 g/cc). We therefore use the cell model together with the PURGATORIO electron-thermal model to build our liquid C free energy.

Final liquid free energy model
The parameters of our liquid cold curve (including break-points) and ion-thermal characteristic temperature are listed in Table 2. Minimal changes have been made to the liquid free energy model of Ref. [9]. These changes include: 1. An altered cold curve at high compressions, 2. The use of the Purgatorio electron-thermal term, and 3. The use of the cell model addition to the Chisolm-Wallace Mie-Grüneisen term (see above). The lower pressure (< 25 GPa) cold curve and the Debye-like ion-thermal model itself is unchanged from that of Ref. [9]. Fig. 9 shows isotherms of internal energy for the liquid as computed by both DFT-MD and our liquid EOS model. The range of V is chosen to correspond to that for which diamond and BC8 phases are stable at lower temperatures. Note that while the lowest temperature (T = 10000 K) isotherm of the model and the DFT-MD are in reasonable accord, especially for smaller V , higher-T isotherms deviate in a manner which depends on V . This results from the Chisolm-Wallace + cell + Purgatorio models being less than perfect for C, a fact which is not surprising given their relative simplicity. Still, as we have discussed, the agreement is far better than that obtained by using the Cowan and/or Thomas-Fermi free energy contributions (see Figs.7 and 8). Pressure isotherms and isochores (Fig. 10) on the same grid of (V, T ) exhibit somewhat better overall agreement. The ultra-high compression pressure isotherms of the liquid are shown in Fig. 11; agreement between our model and our DFT-MD results is quite good but not perfect.

C. phase diagram
Given the free energy models for the individual solid (diamond, BC8, sc, sh) and liquid phases, we construct the phase diagram and multiphase EOS by invoking the Maxwell construction for each pair of phases (here denoted 1 and 2), and at each temperature, T , of interest [47]: where F 1,2 are total free energies for each phase, V 1,2 are the transition volumes, and P is the transition pressure. Since the individual solid phases are only meaningful (as in metastable) in the neighborhood of their fields of thermodynamic stability, we only consider each solid phase free energy function throughout a restricted range of V , so as to minimize the spurious parameter liquid effects of improper extrapolation. In particular, we consider diamond from 2.5Å 3 /atom → 7.0Å 3 /atom, BC8 from 1.3Å 3 /atom → 3.5Å 3 /atom, sc from 1.0Å 3 /atom → 7.0Å 3 /atom, and sh from 0.01Å 3 /atom (undoubtedly well beyond its actually range of stability) → 7.0Å 3 /atom. The liquid is considered throughout the wide range: 0.01Å 3 /atom → 7.0 0.75Å 3 /atom < V < 5.4Å 3 /atom (as we stated above, we choose not to require agreement with the 0.1 g/cc PIMC isochore in this work). Fig. 12a shows the resulting phase diagram in the regime of pressure and temperature addressed in the work of Ref. [9]. Our phase lines computed in this work differ only slightly from those in that earlier study. In particular, the diamond-BC8 transition pressure is slightly lower; this is due to our diamond and BC8 cold curves being slightly different from those of Ref. [9] (higher plane wave cutoffs and different pseudopotentials were used here). Also, the BC8-phase melt temperature is slightly lower, as evidenced by our comparison to the Kechin fits to the carbon melt curves of Ref. [16] to which the EOS model of Ref. [9] was fit. Nevertheless, our phase diagram is quite close to that of the older model in the range of its applicability. Also shown in Fig. 12  The principal Hugoniot of carbon computed with the model is shown as well in Fig. 12.
This is defined to be the locus of final states accessible via a planar one-dimensional shock, given a particular assumed initial state, here taken to be ρ = 3.52 g/cc and T = 300 K.
This curve has two branches: one for which the final state is the diamond phase, and one in which the liquid is the final state. They are separated by a flat region in T straddling the pressure at which the diamond → BC8 transition is predicted to occur. In this portion of the Hugoniot curve the final state is in a mixed-phase region, the size of which is directly related to the latent heat of melting, as discussed in Ref. [9]. We predict a principal Hugoniot which is similar to that predicted by the earlier carbon EOS model [9] (indeed, it is nearly identical in the diamond-phase portion). However, there is a pronounced increase (∼ 100 GPa) in the Hugoniot final state pressure of the liquid portion, which is explained by the added latent heat of melting resulting from the addition of the PURGATORIO electronthermal contribution. Our change in entropy from solid to liquid at constant V at the shock melt conditions is ∼ 3.2k B /atom, while in the model of Ref. [9] it is ∼ 2.9k B /atom. This in turn moves our liquid branch of the Hugoniot a bit closer to the measured T vs. P Hugoniot data of Eggert et al. [6] (the upper and lower bounds of their experimental error bars are displayed as well in Fig. 12a). Our Hugoniot still fails to fall within the bounds of their measurements, however.
The phase diagram of our multiphase C EOS model over a larger range of pressure and temperature is shown in Fig. 12b. The higher-P solid phases, sc and sh, appear as well, in addition to diamond, BC8, and liquid phases. Several important features are worth noting: 1. The BC8 → sc transition pressure is nearly independent of temperature, as is the case for the diamond → BC8 transition pressure. This is because the relevant moments, θ(V ), of the PDOSs of the phases on either side of each transition are very similar at the densities at which the transitions occur, as discussed above. 2. The sc → sh phase line is less vertical, and curves to lower-P for higher T , but we do not predict that this phase line intersects the BC8-sc phase line before intersecting the melt curve. This is in stark contrast to the prediction of Ref. [15], which shows a sc-sh phase line which intersects the BC8-sc line at T ∼ 5000 K. Though we use similar ab initio and EOS modeling methods to theirs, we suspect that their treatment of the T -dependence of the transition pressure in their Gibbs free energy matching prescription is sufficiently different from our approach to produce this discrepancy. We have checked that alternate ways of fitting our EOS data (cold curves, PDOSs and resulting ion-thermal free energy terms) still give the general phase diagram topology we present here, though we do not discount the possibility that the extreme sensitivity of the phase lines to small changes in the individual phase free energies can render our predictions somewhat inaccurate, especially at these higher pressures where the F (V ) curves are nearly parallel. 3. The value of the sc-sh transition pressure at low-T is ∼ 20% higher than that reported in Ref. [15], though the position of the BC8-sc transition line is in reasonable agreement between the two studies. This is likely a result of the use of different pseudopotentials and plane wave convergence criteria. Again, the nearly coincident free energy functions for different phases at high pressures makes the phase lines extremely sensitive to any otherwise subtle changes. 4. The melt line of the sh phase possesses a maximum at P ∼ 10, 000 GPa. This is striking, but is one of the least certain predictions we make in this work; we have seen that slight changes to our EOS model parameters which produce small changes to our liquid free energy (which are within both our assumed model and DFT-MD uncertainties) yield sc melt temperatures which exhibit markedly different behavior.
Our reason for favoring this particular EOS model parameterization over other nearly equivalent ones can be found in our computation of the ionic diffusivities using the DFT-MD. Fig. 13 shows two diffusivity isotherms, the lower one for T = 10,000 K and the upper one for T = 20,000 K. A very small (approaching zero) diffusivity is indicative of a solid phase. The ionic diffusivity at 20,000 K shows the expected decrease as V decreases, but is otherwise large as expected for a liquid. The diffusivity at 10,000 K shows a pronounced dip towards zero at a V ∼ 1.2Å 3 /atom, corresponding to a pressure of roughly 10,000 GPa (see Fig. 11). This is a clear indication of solidification, and indeed, the phase diagram of Fig.   12b shows a melt curve which attains its maximum (of just over 10,000 K) at P = 10, 000 GPa. Though this provides a nice consistency check on our DFT-MD-derived C EOS model, we caution that the diffusivity calculations were performed in a cell of 64 atoms (as were the DFT-MD calculations of EOS), and solidification can be biased in simulations with small cells of fixed shape [59].
The overarching points to be made about our prediction of the phase diagram of C are two-fold: First, our predictions of phase lines are necessarily less accurate the higher the pressure is, for the reasons mentioned directly above. These inaccuracies are the combined result of the inherent uncertainties in the underlying ab initio electronic structure data, and the deficiencies of EOS models which have a finite number of adjustable parameters.
Second, the predicted volume changes across such high-P phase lines are often exceedingly small (which is what makes the accurate determination of these transition lines difficult).
Thus, the least accurate phase lines have a minimal impact on the EOS itself (P and E at a given ρ and T ).

IV. COMPARISONS TO EXPERIMENTAL RESULTS
We have already discussed the comparison of the results of our C EOS model to principal Hugoniot data in the (P, T )-plane. These recent measurements made use of optical pyrometry techniques [6] which lie outside the typical purview of classic shock-compression experiments. Now we compare to a larger set of principal Hugoniot data for elemental carbon transformed more typically to the (ρ, P )-plane, as well as to room-temperature isotherm data from static compression experiments.
Fig. 14 shows the 300 K isotherm (red) and the principal Hugoniot (blue) in (ρ, P )-space as calculated from our multiphase C EOS model. The isotherm (which is very similar to the principal isentrope) shows flat regions where the phase transitions occur (BC8 → sc and sc → sh in this plot). The P vs. ρ Hugoniot is nearly indistinguishable from the isotherm below roughly 1000 GPa, but deviates dramatically from it at higher stresses. The point at which they diverge corresponds to the point where the Hugoniot final states begin to reside in the liquid phase (see Fig. 12). Also shown on this figure are data from static compression measurements of the room-temperature isotherm [60] (green circles at the lowest P ), magnetically-driven flyer plate studies [8] (dense set of magenta + symbols with ρ ∼ 6 -7 g/cc) and multiple sets of laser-shock data on the principal Hugoniot of diamond (symbols with error bars) [1][2][3][4][5]. Note that the much of the highest-P Hugoniot data [1] seems to straddle the 300 K isotherm, even well above P = 1000 GPa; indeed, the data of Brygoo et al. [4] even falls below our prediction of the room-T isotherm. This is very puzzling in light of our theoretical results which suggest that the Hugoniot should be much stiffer than the isotherm at these larger compressions. Our previous ab initio-based C EOS [9] was only fit to DFT-MD data up to T = 20,000 K, which translates to P = 1600 GPa on the principal Hugoniot. At this pressure, the Hugoniot and the 300 K isotherm are still rather close. The In order to check the robustness of our conclusions that the principal Hugoniot should indeed be far stiffer than the 300 K isotherm for P > 1000 GPa, we have verified that the different variants of liquid C EOS we discussed in Section III.B.1 (Cowan + Purgatorio, cell + Thomas-Fermi, etc.) all show this same basic relationship between the low-T isotherm and the Hugoniot. Looking back again to the experimental results pictured in Fig. 14, we note that there are data with P > 1000 GPa that seem to be closer to our Hugoniot predictions, in particular those of Nagao et al. [2] and Hicks et al. [3]. It must be mentioned that the Hicks et al. data shown here is not that as presented in the original Hicks et al.
reference, but is instead a data set which results from a reanalysis [7] of the quartz EOS [61] which was used as an impedance-matched standard in that work. The original uncorrected data set [3] lies much closer to our low-T isotherm than to our principal Hugoniot prediction in this higher-P region. These experimental improvements notwithstanding, we still submit that further work must be done to resolve the theory-experiment discrepancy at the very highest compressions shown here.

V. CONCLUSIONS
We have constructed a 5-phase EOS for elemental carbon based entirely on ab initio electronic structure calculations of the density functional theory and path integral Monte Carlo varieties. The PIMC and high-temperature DFT-MD data were particularly useful in helping us differentiate between different choices of ion-thermal (Cowan, cell) and electronthermal (Thomas-Fermi, PURGATORIO) models. We found that these data strongly favor the cell + PURGATORIO combination. The multiphase EOS is constructed entirely without patching different models together in an ad hoc fashion; rather, full thermodynamic consistency is maintained and the appropriate high-T limit is reached in a seamless manner.
Our EOS model includes, in addition to diamond, BC8, and liquid, two of the ultra-high pressure solid phases recently predicted by Martinez-Canales et al. [15]: simple-cubic, and simple-hexagonal. Though our prediction of the phase diagram of carbon in this high-P region is broadly similar to theirs, the detailed positions of phase lines and resulting triplepoints are somewhat different, owing to the sensitivity of the phase lines to the detailed prescriptions for obtaining the individual phase free energies and computing the transition pressures. Further work should be done to resolve these details, though we maintain that the EOS itself should be reasonably accurate at high compressions in spite of this.
Comparison to recent laser-shock compression data [1][2][3][4][5] shows notable disagreement with the subset of results [1,4] which suggest that the principal Hugoniot is nearly coincident with the low-T isotherm even at pressures in excess of 1000 GPa. It is not at all clear at the moment as to what is causing such a discrepancy, particularly since the general features of our predictions (seen in Fig. 14) are quite robust and seemingly independent of many of the choices we have made in modeling the EOS of liquid C. We are encouraged, however, that a more recent reanalysis [7] of the data of Hicks et al. [3] seems to bring their principal Hugoniot into better agreement with our prediction.
One of the central results of this work is that our high-T DFT-MD and PIMC data on carbon are found to be in excellent agreement with a new ion-thermal free energy model.
is increased which is more rapid than that of the older Cowan model. Although the results presented here are specific to carbon, existing EOS models for other materials that are based on the widely-used Cowan model may need to be re-examined in this light, as those materials may also possess a somewehat rapid decay of C ion V with T . Recent work with simpler classical inter-ionic potentials suggests that such rapid decays may be common [62]. FIG. 14: 300 K isotherm (red curve) and principal Hugoniot (blue curve) as computed by our EOS model. Green points indicate static compression data from Ref. [60], and all other symbols (+ error bars) indicate shock data from Refs. [1][2][3][4][5]8].