High temperature ion-thermal behavior from average-atom calculations

TL;DR

This paper develops a method to calculate ion-thermal behavior at high temperatures using average-atom models, reducing reliance on extensive molecular dynamics simulations, and applies it to carbon to analyze its equation of state.

Contribution

It introduces a new functional form for free energy correction based on the Maxwell-Boltzmann distribution, with a single parameter calibrated by molecular dynamics.

Findings

Ion displacement increases with temperature and compression.

The effective density of potential modes is approximately 0.2 per atom.

The method accurately predicts ion-thermal contributions across a wide range of conditions.

Abstract

Atom-in-jellium calculations of the Einstein frequency were used to calculate the mean displacement of an ion over a wide range of compression and temperature. Expressed as a fraction of the Wigner-Seitz radius, the displacement is a measure of the asymptotic freedom of the ion at high temperature, and thus of the change in heat capacity from 6 to 3 quadratic degrees of freedom per atom. A functional form for free energy was proposed based on the Maxwell-Boltzmann distribution as a correction to the Debye free energy, with a single free parameter representing the effective density of potential modes to be saturated. This parameter was investigated using molecular dynamics simulations, and found to be ~0.2 per atom. In this way, the ion-thermal contribution can be calculated for a wide-range equation of state (EOS) without requiring a large number of molecular dynamics simulations. Example calculations were performed for carbon, including the sensitivity of key EOS loci to ionic freedom.

Equations

$$e_i\frac{m_a}{3k_{\text{B}}T}=x_{\text{E}}\left[\frac{1}{\exp (x_{\text{E}})+1}+\frac{1}{2}\right]=D_3(x_{\text{D}})+\frac{3}{8}{x}_{\text{D}}$$

$$\overline{u^2}\frac{m_a}{3{\hslash }^2}=\frac{1}{{\theta }_{\text{E}}}\left[\frac{1}{\exp (x_{\text{E}})-1}+\frac{1}{2}\right]=\frac{1}{{\theta }_{\text{D}}}\left[\frac{D_1(x_{\text{D}})}{x_{\text{D}}}+\frac{1}{4}\right]$$

$$D_i(x)\equiv \frac{i}{x^i}{\int }_0^x\frac{x^{i}dx}{e^x-1}.$$

$$f_i=k_{\text{B}}T\left[3\ln \left(1-e^{-{\theta }_{\text{D}}/T}\right)+\frac{9{\theta }_{\text{D}}}{8T}-D_3({\theta }_{\text{D}}/T)\right]$$

$$c_v=\frac{3}{2}{k}_{\text{B}}\left\{1+\text{mbox}min\left[1,{\left(\frac{T}{T_m(\rho )}\right)}^{-1/3}\right]\right\}$$

$$f_i=f_b-k_{\text{B}}T\ln \left[\text{mbox}erf\left(\sqrt{\frac{T_r}{T}}\right)-\sqrt{\frac{4T_r}{\pi T}}{e}^{-T_r/T}\right]$$

$$f(v)={\left(\frac{m_a}{2\pi k_{\text{B}}T}\right)}^{3/2}4\pi v^2{e}^{-m_a{v}^2/2k_{\text{B}}T}$$

$$Pr(<v)=\text{mbox}erf\left(\sqrt{\frac{m_a{v}^2}{4k_{\text{B}}T}}\right)-\sqrt{\frac{2m_a{v}^2}{\pi k_{\text{B}}T}}{e}^{-m_a{v}^2/2k_{\text{B}}T}.$$

$$\frac{1}{2}{m}_a{v}^2\simeq \frac{N}{2}{k}_{\text{B}}{T}_b,$$

$$T_b\equiv T:u=r_{\text{WS}}.$$

$$f_i=f_b-k_{\text{B}}T\ln \left[\text{mbox}erf\left(\sqrt{\frac{N{T}_b}{T}}\right)-\sqrt{\frac{4N{T}_b}{\pi T}}{e}^{-N{T}_b/T}\right].$$

$$c_v\equiv {\frac{\partial e}{\partial T}|}_v=T{\frac{\partial s}{\partial T}|}_v,$$

Full text

equation of state

I Introduction

Accurate equations of state (EOS) are essential to understand stellar and planetary formation and evolution, astrophysical impacts, and engineering challenges such as the development of thermonuclear energy sources. However, the behavior of the ionic heat capacity of condensed matter at temperatures between melting and the formation of an ideal plasma is poorly understood, limiting the insight and accuracy of theoretical EOS. Wide-ranging EOS sesleos are almost invariably constructed using empirical models originally derived as approximate representations of observations of the variation of the heat capacity in the liquid of metals of low melting point at one atmosphere Grover1971 , and assumed to apply at arbitrary compression cowan ; Johnson1991 . Experiments to test or improve on this assumption are challenging: where even attempted, uncertainties on measurements of the temperature of warm dense matter are typically greater than 10% (see for instance Kritcher2011 ; Saunders2016 ), which is not adequate to distinguish the details of the ion-thermal heat capacity.

The most rigorous theoretical techniques applicable are path-integral Monte-Carlo (PIMC) and quantum molecular dynamics (QMD), in which the kinetic motion of an ensemble of atoms is simulated, where the distribution of the electrons is found with respect to the changing location of the ions using quantum mechanics pimc ; qmd . The energy of the ensemble is determined from an average over a sufficient time interval, and the heat capacity can be found from the variation of energy with temperature. This procedure is computationally expensive, requiring $o(10^{16})$ or more floating-point operations per state to determine the ionic heat capacity using QMD, equivalent to thousands of CPU-hours per state; PIMC requires roughly an order of magitude more. It is typically deemed impractical to perform these simulations for matter around or below ambient density and above a few tens of electron volts using QMD, limiting the regions of state space over which the EOS can be calibrated in this way.

Recent PIMC and QMD results have indicated that the simpler approach of calculating the electron states for a single atom in a spherical cavity within a uniform charge density of ions and electrons, representing the surrounding atoms, reproduces the electronic component of their more rigorous EOS Benedict2014 ; Driver2017 . This atom-in-jellium approach Liberman1979 has been used previously to predict the electron-thermal energy of matter at high temperatures and compressions, as an advance over Thomas-Fermi and related approaches tf . A development of atom-in-jellium was used to estimate ion-thermal properties using the Debye model Liberman1990 , and we have found that it can be used to construct the complete EOS Swift2018 . However, the model as originally implemented did not account for the decrease in ionic heat capacity at high temperatures as the ions cease to be caged between their neighbors, losing the contribution from potential energy.

In the work reported here, we extend the atom-in-jellium ion-thermal model developed previously to estimate the asymptotic freedom of ions at high temperature, and hence predict the form of the heat capacity of matter in the fluid-plasma regime. Because of the efficiency of average-atom calculations, this approach should for the first time enable EOS to be constructed with consistently accurate electronic contributions and the correct ion-thermal behavior over a wide range of states from expanded to compressed matter and over the full range of temperatures relevant to stellar and planetary interiors, and thermonuclear fusion science.

II Previous ion-thermal models

As condensed matter is heated from absolute zero, the heat capacity of the ions rises from zero as vibrational modes are excited. If all modes are excited before any dissociation occurs, the ionic heat capacity for quadratic degrees of freesom reaches $3k_{\text{B}}$ per atom, where $k_{\text{B}}$ is the Boltzmann constant, representing 3 kinetic and 3 potential modes Waldram1985 . The attractive potential between atoms has a finite depth, and once an ion has more energy, it behaves as a free particle with only kinetic modes available to it, and thus contributes $3k_{\text{B}}/2$ to the heat capacity.

The instantaneous separation and velocity of the ions are described by a distribution, and the most appropriate description of the ion energies is also by a distribution. Even at low temperatures some ions are free, and even at high temperatures some ions are bound. The detailed distribution of ion energies depends on the shape of the interatomic potential as well as the temperature, which complicates the analysis.

EOS have been constructed using very simple ion-thermal models, such as a constant specific heat capacity which may be the value at STP, $3k_{\text{B}}$ per atom, or $3k_{\text{B}}/2$ per atom simplecv . A generalization has been to use the Debye model debye , which assumes a simple form for the phonon density of states (PDOS), proportional to $\hbar\omega^{2}$ for phonon frequencies $0\leq\omega\leq k_{\text{B}}\theta_{\text{D}}$, and 0 for higher $\omega$. This model captures the rise of ion-thermal heat capacity from zero to $3k_{\text{B}}$ per atom as modes are excited, but does not capture the detailed heat capacity arising from the actual PDOS. EOS can be constructed using more accurate PDOS Swift2001 , though, in practice, because integrations are performed over the PDOS, the EOS is not sensitive to the full detail of the spectrum, and high-fidelity EOS have been constructed using a combination of a few Debye frequencies instead multidebye . More importantly for the present study, the Debye model ignores the asymptotic freedom of the ions at high temperature.

In the ion thermal model developed for use with atom-in-jellium calculations Liberman1990 , perturbation theory was used to calculate the Hellmann-Feynman force on the ion when displaced from the center of the cavity in the jellium. Given the force constant $k=-\partial f/\partial r$, the Einstein vibration frequency $\nu_{e}=\sqrt{k/m_{a}}$ was determined, where $m_{a}$ is the atomic mass, and hence the Einstein temperature $\theta_{\text{E}}=h\nu_{e}/k_{\text{B}}$. The Debye temperature $\theta_{\text{D}}$ was inferred from $\theta_{\text{E}}$, by equating either the ion-thermal energy $e_{i}$ or the mean square displacement $\bar{u^{2}}$, giving slightly different results. These calculations are, respectively,

[TABLE]

and

[TABLE]

where $x_{\text{E}}=\theta_{\text{E}}/T$, $x_{\text{D}}=\theta_{\text{D}}/T$, and $D_{i}$ is the Debye integral

[TABLE]

Below, we use the displacement calculation, for consistency. The ion-thermal free energy was then calculated from

[TABLE]

where $\frac{9}{8}k_{\text{B}}\theta_{\text{D}}$ is the zero-point energy. Unusually, $\theta_{\text{D}}$ is a function of temperature as well as compression, effectively because changes in ionization can result in a change to the stiffness and hence the vibration frequency.

An approach used widely in constructing EOS for fluid-plasma applications is the Cowan model cowan , in which the heat capacity is assumed to vary as

[TABLE]

where $T_{m}(\rho)$ is the melting temperature as a function of mass density $\rho$. The effect of the term in brackets is for the potential modes to fall as $(T/T_{m})^{-1/3}$ for $T>T_{m}$, and for $c_{v}$ to be held at $3k_{\text{B}}$ for $T\leq T_{m}$. Variants of the Cowan model have been used with more general dependence on $T/T_{m}$, also treating the latent heat of melting in ad hoc fashion as an increased ionic heat capacity over a finite temperature range Johnson1991 .

QMD studies of carbon Benedict2014 found that the heat capacity dropped more abruptly than in the Cowan model, and were represented better as a free energy of the form

[TABLE]

where $f_{b}$ is the free energy of bound ions and $T_{r}(\rho)$ is a reference temperature curve determined from QMD. This approach appears to be as accurate as QMD can achieve, but is limited by the restricted range of states accessible in practice to QMD.

III Asymptotic freedom of ions in jellium

In the method developed for calculating the vibrations of ions in jellium Liberman1990 , the Debye temperature $\theta_{\text{D}}$ was inferred from the Einstein temperature $\theta_{\text{E}}$ by equating the energy or displacement, as discussed above. The mean displacement $u$ can be used as a measure of ionic freedom: when it exceeds the Wigner-Seitz radius, $r_{\text{WS}}=\left({3m_{a}}/{4\pi\rho}\right)^{1/3}$, the ion is effectively free.

The Inferno program was modified to calculate and output the mean fractional displacement, $u/r_{\text{WS}}(\rho,T)$ (Fig. 1).

An average atom would be bound for $u<r_{\text{WS}}$, with ionic heat capacity $3k_{\text{B}}$, and then free for $u\geq r_{\text{WS}}$, with heat capacity $3k_{\text{B}}/2$. Given an energy distribution for the atoms, the change in heat capacity can be represented more accurately as a continuous variation with temperature. An accurate distribution could be calculated from the set of available ion-thermal energy levels populated using Boltzmann factors, but the average-atom-in-jellium method gives only a rough approximation to the states and hence to the energy levels. A simpler estimate can be made using the Maxwell-Boltzmann distribution for the velocity $v$ of free ions,

[TABLE]

Waldram1985 . The cumulative probability distribution, giving the fraction of ions whose velocity is less than $v$, is

[TABLE]

The critical displacement $u=r_{\text{WS}}$ can be equated to a critical velocity $v$, and hence to a characteristic binding temperature,

[TABLE]

where $N$ is the number of modes that must be saturated for the ion to become free, i.e.

[TABLE]

We have not so far found an integral of the Maxwell-Boltzmann derived heat capacity to represent the free energy in closed form, but we can generalize the similar functional form used previously Benedict2014 and derived from the partition function of a particle in a harmonic potential of finite volume Correa2018 , which we have shown exhibits the desired temperature dependence for the heat capacity, giving a modification to the Debye free energy (or any other free energy $f_{b}$ representing bound ions) so that the heat capacity falls at high temperature,

[TABLE]

This approach is similar, except that $T/T_{b}(\rho,T)$ is determined directly from the average-atom calculations, and the mode number $N$ is likely to be a constant $o(1)$ with $\rho$, $T$, and atom type. In contrast, the approach used in the previous study Benedict2014 requires QMD simulations to be performed for at least a few temperatures and over the full range of densities of interest, for each substance.

IV Quantum molecular dynamics simulations

QMD simulations were performed to predict the variation of ionic heat capacity $c_{vi}$ directly. Such simulations treat the motion of the ions as classical, with 3 kinetic modes. The total potential energy is calculated from the electron states with respect to the instantaneous configuration of the ions at the temperature of interest. The potential contribution to the $c_{vi}$ must thus be inferred from the total $c_{v}$, by calculating and subtracting the electronic heat capacity $c_{ve}$ Whitley2015 . Along each isochore, $c_{vi}$ was found to fall just as $c_{ve}$ started to rise, so the deduced variation in $c_{vi}(T)$ was sensitive to the treatment of $c_{ve}$. In addition, the drop in $c_{vi}$ started to become appreciable a little above the melt locus, so its precise variation depended on discriminating the latent heat of melting, which may be spread out in temperature in a relatively small ensemble of atoms.

QMD simulations were performed using the electronic structure program Vasp vasp . The projector augmented wave method paw was used, with carbon ions represented with a pseudopotential subsuming the inner two electrons. Electron wavefunctions were represented with a plane wave basis set cut off at 1000 eV, at the Baldereschi mean-value point in reciprocal space Baldereschi1972 . Density functional theory in the local density approximation Hohenberg1964 ; Kohn1965 ; Perdew1992 ; Perdew1994 was used for the exchange-correlation energy; some states were recalculated with the Perdew-Burke-Ernzerhof functional Perdew1996 for comparison. The simulations were in the NVT ensemble, using the Nosé-Hoover thermostat Hoover1996 , with periodic boundary conditions. For each state of density and temperature, the motion of 64 atoms was integrated for 20000-50000 steps of 0.05-0.5 fs. Convergence with respect to plane-wave cutoff energy was tested at 2000 K and 5 g/cm3 with calculations up to 3000 eV, and found to be converged to $<$0.06% in pressure and $<$4 meV/atom. Overall, given the finite cell size and use of the single mean-value $k$-point, the pressure is likely to be converged to $<$2% and the energy to $<$10 meV/atom. Simulations were performed along isochores at 5 and 10 g/cm3. Each simulation yielded a value for the total energy $e$ and the electronic entropy $s_{e}$, the latter from the population of electron states. Using

[TABLE]

the total heat capacity $c_{v}$ was deduced from $e(T)$, and the electronic heat capacity $c_{ve}$ from $s_{e}(T)$, by fitting polynomials to sections of each isotherm, and differentiating. Given the numerical uncertainties in convergence, differentiation, and subtraction, the uncertainty in ionic heat capacity was around 0.25 $k_{\text{B}}$/atom. (Fig. 2.)

To interpret the QMD results, the temperature along each isochore was expressed as $u^{2}/r_{\text{WS}}^{2}(\rho,T)\equiv T/T_{b}$ from the atom-in-jellium calculations. The best fit of the hypothesized free energy function, Eq. 11, was found for $N=0.2$ (Fig. 3.)

V Equation of state for carbon

Several wide-range EOS have been constructed for C, mostly using Thomas-Fermi theory for high compressions and temperatures. We compare against one such EOS employing a Cowan-like ion-thermal model, sesame EOS 7834 Crockett2006 , which was constrained by QMD calculations and Hugoniot measurements up to shock melting, but does not include the melting transition explicitly. In contrast, the 5-phase EOS Benedict2014 used atom-in-jellium calculations for the electron-thermal contribution only, and was constructed to include an explicit melting transition.

The principal shock Hugoniot deduced by numerical solution Swift2008 ; Swift2018a for each EOS. The treatment of ionic freedom caused a variation of up to 20% in pressure between 6 and 12 g/cm3 (1 and 10 TPa), large enough to investigate experimentally. (Figs. 4 and 5.)

VI Conclusions

Atom-in-jellium calculations of thermal vibrations of the ion were used to construct a model of the reduction in ionic heat capacity at high temperatures as the ions become free. The model has a single free parameter, equivalent to the effective number of potential modes that must be saturated before the high-energy tail of the Maxwell-Boltzmann distribution describes free particles. This parameter was investigated using molecular dynamics simulations and found to be $0.2\pm 0.03$, with a weak dependence on compression.

Carbon was chosen as a prototype material as its atomic number is low, emphasizing changes in the ionic heat capacity compared with the electrons, and it is widely used as a sample and ablator in high pressure experiments. The principal shock Hugoniot and ambient isochore showed a modest sensitivity to the treatment of ionic heat capacity, which may be experimentally detectable.

Atom-in-jellium calculations of electronic states were previously shown to be as accurate for warm dense matter as the most rigorous approaches of path-integral Monte Carlo and quantum molecular dynamics. The work reported here makes the atom-in-jellium based treatment of ion-thermal energies as complete as these multi-atom calculations, and means it is possible for the first time to calculate all contributions to the equation of state self-consistently, and efficiently enough to construct an entire wide-range equation of state for an element in a few CPU-hours at most.

Acknowledgments

This work was performed under the auspices of the U.S. Department of Energy under contract DE-AC52-07NA27344.