High temperature ion-thermal behavior from average-atom calculations
Damian C. Swift, Mandy Bethkenhagen, Alfredo A. Correa, Thomas Lockard, Sebastien Hamel, Lorin X. Benedict, Philip A. Sterne, Bard I. Bennett
aa r X i v : . [ phy s i c s . c o m p - ph ] M a y High temperature ion-thermal behavior from average-atom calculations
Damian C. Swift, ∗ Mandy Bethkenhagen † , Alfredo A. Correa, Thomas Lockard, Sebastien Hamel, Lorin X. Benedict, Philip A. Sterne, and Bard I. Bennett Lawrence Livermore National Laboratory, 7000 East Avenue, Livermore, California 94551, USA Los Alamos National Laboratory, PO Box 1663, Los Alamos, New Mexico 87545, USA (Dated: November 1, 2018; revisions to March 13, 2019 – LLNL-JRNL-768634)Atom-in-jellium calculations of the Einstein frequency were used to calculate the mean displace-ment of an ion over a wide range of compression and temperature. Expressed as a fraction of theWigner-Seitz radius, the displacement is a measure of the asymptotic freedom of the ion at hightemperature, and thus of the change in heat capacity from 6 to 3 quadratic degrees of freedom peratom. A functional form for free energy was proposed based on the Maxwell-Boltzmann distributionas a correction to the Debye free energy, with a single free parameter representing the effective den-sity of potential modes to be saturated. This parameter was investigated using molecular dynamicssimulations, and found to be ∼ Keywords: equation of state
INTRODUCTION
Accurate equations of state (EOS) are essential to un-derstand stellar and planetary formation and evolution,astrophysical impacts, and engineering challenges such asthe development of thermonuclear energy sources. How-ever, the behavior of the ionic heat capacity of condensedmatter at temperatures between melting and the forma-tion of an ideal plasma is poorly understood, limiting theinsight and accuracy of theoretical EOS. Wide-rangingEOS [1] are almost invariably constructed using empiri-cal models originally derived as approximate representa-tions of observations of the variation of the heat capac-ity in the liquid of metals of low melting point at oneatmosphere [2], and assumed to apply at arbitrary com-pression [3, 4]. Experiments to test or improve on thisassumption are challenging: where even attempted, un-certainties on measurements of the temperature of warmdense matter are typically greater than 10% (see for in-stance [6, 7]), which is not adequate to distinguish thedetails of the ion-thermal heat capacity.The most rigorous theoretical techniques applicable arepath-integral Monte-Carlo (PIMC) and quantum molec-ular dynamics (QMD), in which the kinetic motion ofan ensemble of atoms is simulated, where the distribu-tion of the electrons is found with respect to the chang-ing location of the ions using quantum mechanics [8, 9].The energy of the ensemble is determined from an aver-age over a sufficient time interval, and the heat capac-ity can be found from the variation of energy with tem-perature. This procedure is computationally expensive,requiring o (10 ) or more floating-point operations per † Present affiliation: Universit¨at Rostock, 18051 Rostock, Germany state to determine the ionic heat capacity using QMD,equivalent to thousands of CPU-hours per state; PIMCrequires roughly an order of magitude more. It is typi-cally deemed impractical to perform these simulations formatter around or below ambient density and above a fewtens of electron volts using QMD, limiting the regions ofstate space over which the EOS can be calibrated in thisway.Recent PIMC and QMD results have indicated thatthe simpler approach of calculating the electron statesfor a single atom in a spherical cavity within a uniformcharge density of ions and electrons, representing the sur-rounding atoms, reproduces the electronic component oftheir more rigorous EOS [10, 11]. This atom-in-jelliumapproach [12] has been used previously to predict theelectron-thermal energy of matter at high temperaturesand compressions, as an advance over Thomas-Fermiand related approaches [13]. A development of atom-in-jellium was used to estimate ion-thermal properties usingthe Debye model [14], and we have found that it can beused to construct the complete EOS [15]. However, themodel as originally implemented did not account for thedecrease in ionic heat capacity at high temperatures asthe ions cease to be caged between their neighbors, losingthe contribution from potential energy.In the work reported here, we extend the atom-in-jellium ion-thermal model developed previously to esti-mate the asymptotic freedom of ions at high temperature,and hence predict the form of the heat capacity of mat-ter in the fluid-plasma regime. Because of the efficiencyof average-atom calculations, this approach should forthe first time enable EOS to be constructed with consis-tently accurate electronic contributions and the correction-thermal behavior over a wide range of states from ex-panded to compressed matter and over the full range oftemperatures relevant to stellar and planetary interiors,and thermonuclear fusion science.
PREVIOUS ION-THERMAL MODELS
As condensed matter is heated from absolute zero, theheat capacity of the ions rises from zero as vibrationalmodes are excited. If all modes are excited before anydissociation occurs, the ionic heat capacity for quadraticdegrees of freesom reaches 3 k B per atom, where k B is theBoltzmann constant, representing 3 kinetic and 3 poten-tial modes [16]. The attractive potential between atomshas a finite depth, and once an ion has more energy, it be-haves as a free particle with only kinetic modes availableto it, and thus contributes 3 k B / k B per atom, or3 k B / hω forphonon frequencies 0 ≤ ω ≤ k B θ D , and 0 for higher ω .This model captures the rise of ion-thermal heat capacityfrom zero to 3 k B per atom as modes are excited, but doesnot capture the detailed heat capacity arising from theactual PDOS. EOS can be constructed using more accu-rate PDOS [19], though, in practice, because integrationsare performed over the PDOS, the EOS is not sensitiveto the full detail of the spectrum, and high-fidelity EOShave been constructed using a combination of a few De-bye frequencies instead [20]. More importantly for thepresent study, the Debye model ignores the asymptoticfreedom of the ions at high temperature.In the ion thermal model developed for use with atom-in-jellium calculations [14], perturbation theory was usedto calculate the Hellmann-Feynman force on the ion whendisplaced from the center of the cavity in the jellium.Given the force constant k = − ∂f /∂r , the Einstein vibra-tion frequency ν e = p k/m a was determined, where m a is the atomic mass, and hence the Einstein temperature θ E = hν e /k B . The Debye temperature θ D was inferredfrom θ E , by equating either the ion-thermal energy e i orthe mean square displacement ¯ u , giving slightly differentresults. These calculations are, respectively, e i m a k B T = x E (cid:20) x E ) + 1 + 12 (cid:21) = D ( x D ) + 38 x D (1) and¯ u m a h = 1 θ E (cid:20) x E ) − (cid:21) = 1 θ D (cid:20) D ( x D ) x D + 14 (cid:21) (2)where x E = θ E /T , x D = θ D /T , and D i is the Debyeintegral D i ( x ) ≡ ix i Z x x i dxe x − . (3)Below, we use the displacement calculation, for consis-tency. The ion-thermal free energy was then calculatedfrom f i = k B T (cid:20) (cid:16) − e − θ D /T (cid:17) + 9 θ D T − D ( θ D /T ) (cid:21) (4)where k B θ D is the zero-point energy. Unusually, θ D isa function of temperature as well as compression, effec-tively because changes in ionization can result in a changeto the stiffness and hence the vibration frequency.An approach used widely in constructing EOS for fluid-plasma applications is the Cowan model [3], in which theheat capacity is assumed to vary as c v = 32 k B ( " , (cid:18) TT m ( ρ ) (cid:19) − / (5)where T m ( ρ ) is the melting temperature as a function ofmass density ρ . The effect of the term in brackets is forthe potential modes to fall as ( T /T m ) − / for T > T m ,and for c v to be held at 3 k B for T ≤ T m . Variants of theCowan model have been used with more general depen-dence on T /T m , also treating the latent heat of meltingin ad hoc fashion as an increased ionic heat capacity overa finite temperature range [4].QMD studies of carbon [10] found that the heat ca-pacity dropped more abruptly than in the Cowan model,and were represented better as a free energy of the form f i = f b − k B T ln " erf r T r T ! − r T r πT e − T r /T (6)where f b is the free energy of bound ions and T r ( ρ ) is areference temperature curve determined from QMD. Thisapproach appears to be as accurate as QMD can achieve,but is limited by the restricted range of states accessiblein practice to QMD. ASYMPTOTIC FREEDOM OF IONS IN JELLIUM
In the method developed for calculating the vibra-tions of ions in jellium [14], the Debye temperature θ D was inferred from the Einstein temperature θ E byequating the energy or displacement, as discussed above. )0.010.1110100100010000 t e m pe r a t u r e ( e V ) FIG. 1: Mean fractional displacement calculated for carbon.Contours shown are from 0.1 to 1.0 at intervals of 0.1, andthen 2.0 to 4.0 at intervals of 1.0. A fractional displacementof 1 would correspond to ionic freedom, ignoring the velocitydistribution of the ions.
The mean displacement u can be used as a measure ofionic freedom: when it exceeds the Wigner-Seitz radius, r WS = (3 m a / πρ ) / , the ion is effectively free.The Inferno program was modified to calculate andoutput the mean fractional displacement, u/r WS ( ρ, T )(Fig. 1).An average atom would be bound for u < r WS , withionic heat capacity 3 k B , and then free for u ≥ r WS , withheat capacity 3 k B /
2. Given an energy distribution forthe atoms, the change in heat capacity can be representedmore accurately as a continuous variation with tempera-ture. An accurate distribution could be calculated fromthe set of available ion-thermal energy levels populatedusing Boltzmann factors, but the average-atom-in-jelliummethod gives only a rough approximation to the statesand hence to the energy levels. A simpler estimate canbe made using the Maxwell-Boltzmann distribution forthe velocity v of free ions, f ( v ) = (cid:18) m a πk B T (cid:19) / πv e − m a v / k B T (7)[16]. The cumulative probability distribution, giving thefraction of ions whose velocity is less than v , is P r ( < v ) = erf s m a v k B T − s m a v πk B T e − m a v / k B T . (8)The critical displacement u = r WS can be equated to acritical velocity v , and hence to a characteristic bindingtemperature, 12 m a v ≃ N k B T b , (9)where N is the number of modes that must be saturated for the ion to become free, i.e. T b ≡ T : u = r WS . (10)We have not so far found an integral of the Maxwell-Boltzmann derived heat capacity to represent the freeenergy in closed form, but we can generalize the similarfunctional form used previously [10] and derived from thepartition function of a particle in a harmonic potentialof finite volume [5], which we have shown exhibits thedesired temperature dependence for the heat capacity,giving a modification to the Debye free energy (or anyother free energy f b representing bound ions) so that theheat capacity falls at high temperature, f i = f b − k B T ln " erf r N T b T ! − r N T b πT e − NT b /T . (11)This approach is similar, except that T /T b ( ρ, T ) is deter-mined directly from the average-atom calculations, andthe mode number N is likely to be a constant o (1) with ρ , T , and atom type. In contrast, the approach used inthe previous study [10] requires QMD simulations to beperformed for at least a few temperatures and over thefull range of densities of interest, for each substance. QUANTUM MOLECULAR DYNAMICSSIMULATIONS
QMD simulations were performed to predict the varia-tion of ionic heat capacity c vi directly. Such simulationstreat the motion of the ions as classical, with 3 kineticmodes. The total potential energy is calculated from theelectron states with respect to the instantaneous config-uration of the ions at the temperature of interest. Thepotential contribution to the c vi must thus be inferredfrom the total c v , by calculating and subtracting the elec-tronic heat capacity c ve [21]. Along each isochore, c vi wasfound to fall just as c ve started to rise, so the deducedvariation in c vi ( T ) was sensitive to the treatment of c ve .In addition, the drop in c vi started to become apprecia-ble a little above the melt locus, so its precise variationdepended on discriminating the latent heat of melting,which may be spread out in temperature in a relativelysmall ensemble of atoms.QMD simulations were performed using the electronicstructure program Vasp [22]. The projector augmentedwave method [23] was used, with carbon ions representedwith a pseudopotential subsuming the inner two elec-trons. Electron wavefunctions were represented with aplane wave basis set cut off at 1000 eV, at the Baldereschimean-value point in reciprocal space [24]. Density func-tional theory in the local density approximation [25–28]was used for the exchange-correlation energy; some stateswere recalculated with the Perdew-Burke-Ernzerhof func-tional [29] for comparison. The simulations were in the hea t c apa c i t y ( k B / a t o m ) temperature (kK)totalelectronicionic5 g/cm
10 g/cm FIG. 2: Variation of heat capacities along sample isochores.Fully activated vibrational modes give 3 k B /atom, and un-bound kinetic modes k B /atom. NVT ensemble, using the Nos´e-Hoover thermostat [30],with periodic boundary conditions. For each state ofdensity and temperature, the motion of 64 atoms wasintegrated for 20000-50000 steps of 0.05-0.5 fs. Conver-gence with respect to plane-wave cutoff energy was testedat 2000 K and 5 g/cm with calculations up to 3000 eV,and found to be converged to < < k -point, the pressure is likely tobe converged to <
2% and the energy to <
10 meV/atom.Simulations were performed along isochores at 5 and10 g/cm . Each simulation yielded a value for the totalenergy e and the electronic entropy s e , the latter fromthe population of electron states. Using c v ≡ ∂e∂T (cid:12)(cid:12)(cid:12)(cid:12) v = T ∂s∂T (cid:12)(cid:12)(cid:12)(cid:12) v , (12)the total heat capacity c v was deduced from e ( T ), andthe electronic heat capacity c ve from s e ( T ), by fittingpolynomials to sections of each isotherm, and differenti-ating. Given the numerical uncertainties in convergence,differentiation, and subtraction, the uncertainty in ionicheat capacity was around 0.25 k B /atom. (Fig. 2.)To interpret the QMD results, the temperature alongeach isochore was expressed as u /r ( ρ, T ) ≡ T /T b from the atom-in-jellium calculations. The best fit ofthe hypothesized free energy function, Eq. 11, was foundfor N = 0 . EQUATION OF STATE FOR CARBON
Several wide-range EOS have been constructed for C,mostly using Thomas-Fermi theory for high compressions hea t c apa c i t y ( k B / a t o m ) temperature (kK) N=0.2QMD FIG. 3: Variation of atom-in-jellium ion-thermal heat capac-ity along sample isochores. p r e ss u r e ( T P a ) mass density (g/cm )Cowan (SESAME 7834)Benedict 20140.10.20.5 FIG. 4: Principal shock Hugoniot for carbon calculated usingdifferent EOS, including atom-in-jellium with different valuesof the parameter N in Eq.11, and compared with experimentaldata [34]. and temperatures. We compare against one such EOSemploying a Cowan-like ion-thermal model, sesame EOS7834 [33], which was constrained by QMD calculationsand Hugoniot measurements up to shock melting, butdoes not include the melting transition explicitly. In con-trast, the 5-phase EOS [10] used atom-in-jellium calcula-tions for the electron-thermal contribution only, and wasconstructed to include an explicit melting transition.The principal shock Hugoniot deduced by numericalsolution [31, 32] for each EOS. The treatment of ionicfreedom caused a variation of up to 20% in pressure be-tween 6 and 12 g/cm (1 and 10 TPa), large enough toinvestigate experimentally. (Figs. 4 and 5.) p r e ss u r e ( T P a ) temperature (eV)Cowan (SESAME 7834)Benedict 20140.10.20.5 FIG. 5: Ambient isochore for carbon from different EOS, in-cluding atom-in-jellium with different values of the parameter N in Eq.11. CONCLUSIONS
Atom-in-jellium calculations of thermal vibrations ofthe ion were used to construct a model of the reduc-tion in ionic heat capacity at high temperatures as theions become free. The model has a single free parame-ter, equivalent to the effective number of potential modesthat must be saturated before the high-energy tail of theMaxwell-Boltzmann distribution describes free particles.This parameter was investigated using molecular dynam-ics simulations and found to be 0 . ± .
03, with a weakdependence on compression.Carbon was chosen as a prototype material as itsatomic number is low, emphasizing changes in the ionicheat capacity compared with the electrons, and it iswidely used as a sample and ablator in high pressureexperiments. The principal shock Hugoniot and ambientisochore showed a modest sensitivity to the treatmentof ionic heat capacity, which may be experimentally de-tectable.Atom-in-jellium calculations of electronic states werepreviously shown to be as accurate for warm dense matteras the most rigorous approaches of path-integral MonteCarlo and quantum molecular dynamics. The work re-ported here makes the atom-in-jellium based treatmentof ion-thermal energies as complete as these multi-atomcalculations, and means it is possible for the first time tocalculate all contributions to the equation of state self-consistently, and efficiently enough to construct an en-tire wide-range equation of state for an element in a fewCPU-hours at most.
Acknowledgments
This work was performed under the auspices of theU.S. Department of Energy under contract DE-AC52-07NA27344. ∗ Electronic address: [email protected][1] Prominent examples are the Los Alamos and LawrenceLivermore National Laboratories’ EOS libraries:S.P. Lyon and J.D. Johnson, Los Alamos NationalLaboratory report LA-UR-92-3407 (1992); R.M. More,K.H. Warren, D.A. Young and G.B. Zimmerman, Phys.Fluids , 3059 (1988); D.A. Young and E.M. Corey,J. Appl. Phys. , 3748 (1995).[2] R. Grover, J. Chem. Phys. , 3435 (1971).[3] C.W. Cranfill and R. More, Los Alamos Scientific Labo-ratory report LA-7313-MS (1978); R.M. More, K.H. War-ren, D.A. Young, and G.B. Zimmerman, Phys. Fluids ,3059 (1988).[4] J.D. Johnson, High Press. Res. (5), 277-285 (1991).[5] A.A. Correa, L.X. Benedict, M.A. Morales, P.A. Sterne,J. I. Castor, and E. Schwegler, arXiv:1806.01346 (2018).[6] A.L. Kritcher, T. D¨oppner, C. Fortmann, T. Ma,O.L. Landen, R. Wallace, and S.H. Glenzer, Phys. Rev.Lett. , 015002 (2011).[7] A.M. Saunders, A. Jenei, T. D¨oppner, R.W. Falcone,D. Kraus, A. Kritcher, O.L. Landen, J. Nilsen, andD. Swift, Rev. Sci. Instrum. , 11E724 (2016).[8] E.L. Pollock and D.M. Ceperley, Phys. Rev. B , 2555(1984).[9] For example, L. Collins, I. Kwon, J. Kress, N. Troullier,and D. Lynch, Phys. Rev. E , 6202 (1995).[10] L.X. Benedict, K.P. Driver, S. Hamel, B. Militzer, T. Qi,A.A. Correa, A. Saul, and E. Schwegler, Phys. Rev. B , 224109 (2014). Note: in Eq. 11 of this reference, theexponent is positive, which gives a heat capacity of 2 . k B as T →
0. The negative exponent gives 3 k B , as desired.[11] K.P. Driver and B. Militzer, Phys. Rev. E , 043205(2017).[12] D.A. Liberman, Phys. Rev. B , 12, 4981 (1979).[13] L.H. Thomas, Proc. Cambridge Phil. Soc. , 5, 542548(1927); E. Fermi, Rend. Accad. Naz. Lincei. , 602607(1927).[14] D.A. Liberman and B.I. Bennett, Phys. Rev. B , 2475(1990).[15] D.C. Swift, T. Lockard, M. Bethkenhagen, R.G. Kraus,L.X. Benedict, P. Sterne, M. Bethkenhagen, S. Hamel,and B.I. Bennett, submitted and arXiv:1903.00163 (2018).[16] J.R. Waldram, “The Theory of Thermodynamics”, Cam-bridge (1985).[17] For example, D.J. Steinberg, Lawrence Livermore Na-tional Laboratory report UCRL-MA-106439 change 1(1996).[18] P. Debye, Ann. Phys. , 4, 789839 (1912).[19] D.C. Swift, G.J. Ackland, A. Hauer, and G.A. Kyrala,Phys. Rev. B , 214107, (2001).[20] A.A. Correa, L.X. Benedict, D.A. Young, E. Schwegler,and S.A. Bonev, Phys. Rev. B , 024101 (2008). [21] H.D. Whitley, D.M. Sanchez, S. Hamel, A.A. Correa, andL.X. Benedict, Contrib. Plasma Phys. , 5, pp. 390-398(2015).[22] G. Kresse and J. Furthm¨uller, Phys. Rev. B , 11169(1996).[23] P.E. Bl¨ochl, Phys. Rev. B , 17953 (1994).[24] A. Baldereschi, Phys. Rev. B , 5212 (1972).[25] P. Hohenberg and W. Kohn, Phys. Rev. B , 3B,pp 864–871 (1964).[26] W. Kohn and L.J. Sham, Phys Rev , 4A, pp 1133–1138 (1965).[27] J. Perdew, Phys Rev B46 , 6671 (1992).[28] J. Perdew, Phys Rev
B50 , 4954 (1994).[29] J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996).[30] W.G. Hoover and B.L. Holian, Phys. Lett. A , 5,pp 253257 (1996). [31] D.C. Swift, J. Appl. Phys. , 073536 (2008).[32] D.C. Swift and M. Millot, in preparation.[33] S. Crockett, documentation for sesame EOS 7834, LosAlamos National Laboratory (2006).[34] M.N. Pavlovskii, Sov. Phys. Solid State , 741 (1971);K. Kondo and T.J. Ahrens, Geophys. Res. Lett. ,281 (1983); D.K. Bradley et al, Phys. Rev. Lett. ,19, 195506 (2004); H. Nagao et al., Phys. Plasmas ,052705 (2006); S. Brygoo et al., Nat. Mater. , 274(2007); D.G. Hicks et al, Phys. Rev. B , 174102 (2008);R.S. McWilliams et al, Phys. Rev. B , 014111 (2010);M. Knudson et al, Science , 1822 (2008); M.C. Gregoret al, Phys. Rev. B95