Unified first-principles equations of state of deuterium-tritium mixtures in the global inertial confinement fusion region

Inertial confinement fusion (ICF) is one of the most promising approaches to achieving an unlimited supply of clean energy. In conventional central ignition designs, cryogenic deuterium–tritium (DT) fuel is compressed to a state of high density and high temperature by an imploding ablator driven by strong sources such as intense laser pulses, x-rays generated by laser ablation, and Z-pinches.^1,2 In the process of ICF, the imploding DT fuel goes from a cold condensed state to one of warm dense matter, finally reaching the hot dense plasma regime, where the density covers a wide range from 0.1 g/cm³ to 1000 g/cm³ and the temperature varies from several hundred kelvin to a few thousand electronvolts.³ Accurate knowledge of the thermodynamic properties of DT fuel such as its equation of state (EOS) and its transport coefficient in these wide density and temperature ranges is essential for ICF designs using hydrodynamic simulations.^4–13 

To obtain an accurate EOS of DT mixtures, much effort has been devoted to measuring the Hugoniot and related properties of hydrogen and deuterium under shock compressions driven by gas guns,¹⁴ converging shocks,¹⁵ high-power lasers,^16–18 and magnetically driven fliers.^19,20 However, owing to the complex nature of condensed hydrogen and the challenges faced in the development of suitable diagnostic techniques, there is still a scarcity of experimental data. To date, the pressures reached in the laboratory are limited to several megabars, and the accuracy of shock-compression data is not yet sufficient to establish the reliability of various theoretical models. Consequently, theoretical calculations have become the most important approach to obtaining EOS data over the wide ranges of density and temperature relevant to ICF.¹² 

A number of theoretical methods have been developed for calculating the EOS of matter under extreme conditions. EOS data for hydrogen and its isotopes generated by chemical models^21–24 and the SESAME EOS tables^25,26 are widely used in radiation hydrodynamic simulations because these models are computationally efficient. To accurately describe the electronic structure of hot dense plasmas, the average atom (AA) model^27,28 was developed under the assumption of a single-particle spherically symmetric ionic potential. In this model, a pseudo-atom with an average fractional occupation number for electron orbitals is used to approximately describe the ions in the plasma environment. It should be noted that at high densities, ion–ion interactions could break the spherical symmetry of the ionic potential, and ionic correlation effects must therefore be taken into account. A method that combines the AA model and molecular dynamics simulations (AAMD)²⁹ has been proposed to treat ionic correlation effects at the level of pair correlations for calculating EOS data for hot dense plasmas. Although these methods have been extensively applied to the high-energy-density plasma regime, they cannot provide a satisfactory description of the strong-coupling state of DT fuel that exists in ICF processes.³⁰ 

At present, two ab initio approaches, namely, density functional theory^31–33 (DFT)-based quantum molecular dynamics (QMD)^34–36 and quantum Monte Carlo (QMC) methods,^37–40 are most widely used to calculate the EOS of materials at high densities and temperatures. One of the most promising QMC methods is path-integral Monte Carlo (PIMC),^37,38 which treats ions and electrons quantum-mechanically on the same footing. Another QMC algorithm is coupled electron–ion Monte Carlo (CEIMC),^39,40 which uses the conventional QMC method to obtain the potential energy surface directly. Although both of these have been employed to investigate the thermodynamic properties of hydrogen and its isotopes, the CEIMC predictions disagree with the experimental principal Hugoniot of deuterium,⁴¹ whereas PIMC becomes computationally prohibitive at low temperatures. On the other hand, within the Kohn–Sham–Mermin DFT framework, the ionic strong-coupling effect, which is significant in the high-density region of the DT phase diagram, can be included naturally in QMD, and the key approximation in principle is the exchange-correlation functional. To overcome the prohibitive computational cost of QMD at extremely high temperatures, orbital-free molecular dynamics (OFMD),^42–44 which constructs the approximate noninteracting free energy functional without the assistance of single-electron wavefunctions, and an extended QMD with a plane-wave approximation at high energy have been proposed.⁴⁵ Although OFMD is highly efficient, it is inaccurate at low temperatures because it lacks electron orbital information. Some EOS data for hydrogen and deuterium have been obtained by combining the two methods, using QMD at low temperature and a semiclassical method such as OFMD at high temperatures.⁶ In this case, the location of the boundary and the transition between the two methods is an essential but challenging task. The quantum Langevin molecular dynamics (QLMD) method was developed for a unified description of matter over a wide range from the cold condensed state to the ideal plasma gas.⁴⁶ It not only considers the electron–ion collision effects at high temperatures, which is usually neglected in QMD simulations, but also has a lower computational cost than conventional QMD.⁴⁷ QLMD has been successfully applied to calculate the wide-range EOS of hydrogen, hydrogen–helium mixtures, and iron.^48–52 

In this work, we perform extensive simulations to calculate the pressure and internal energy of DT mixtures over wide ranges of density and temperature using ab initio QLMD simulations. In contrast to EOS tables obtained by combining various theoretical methods each of which is suitable for different density and temperature ranges, the EOS data presented in this work, which can be used for hydrodynamic simulations of ICF implosions, are obtained for the first time within a unified ab initio framework. We also compare the EOS data obtained from QLMD simulations with those from AAMD and OFMD simulations to assess the accuracy of these data and the regions of validity of these methods.

In this section, we briefly introduce the QLMD simulation method. In conventional QMD simulations, ions move on the smooth potential surface obtained from Kohn–Sham DFT (KSDFT) calculations of electronic structure within the framework of the Born–Oppenheimer approximation. The physical quantities are averaged over all the configurations along the MD trajectories after a thermalization process. QMD has been extensively applied in a variety of fields from cold condensed matter to warm dense matter. We should note that when matter is in a warm or even a hot dense state, a large number of electrons are excited or ionized. These nearly free electrons form a sea of electrons in warm or hot dense matter, and there are frequent elastic or inelastic electron–ion collisions. There is an analogy between ions in warm or hot dense matter and heavy particles in Brownian motion. In the warm or hot dense regime, ions move in the electron sea as heavy Brownian particles, and electron–ion collisions occur frequently. The effects of these collisions, which are not included in conventional adiabatic QMD simulations, play important roles in determining the structures and thermodynamic properties of warm dense matter. We introduce such electron–ion collision-induced friction (EI-CIF) into the ion dynamics within an adiabatic framework, and we describe the ion motion using the Langevin equation, which takes the form

where M_I is the ion mass, R_I is the ion position, F is the force calculated from DFT, γ is the friction coefficient, and N_I is a Gaussian random force.

As a key parameter, the friction coefficient γ plays a central role in QLMD. There are three contributions to γ: γ = γ_B + γ_f + γ_a. The most important of these, γ_B, represents electron–ion collisions and is derived according to the assumptions of the Rayleigh model,⁵³ 

where m_e is the electron mass, n_i is the ion number density, and Z^* is the average degree of ionization, which is obtained by another approach such as the average atom model. The second contribution, γ_f, arises from force errors, i.e., it is the difference between the force obtained with insufficient convergence and that obtained with sufficient convergence in otherwise identical self-consistent-field calculations.⁴⁷ The Gaussian distribution of the force errors makes it possible to accelerate the QMD simulation with the Langevin equation. The third contribution, γ_a, is generally used as a conventional Langevin thermostat parameter to maintain a constant temperature. In the warm or hot dense regime, in particular at high temperatures, the electron–ion friction coefficient γ_B makes the dominant contribution to ion motion, and thus γ_f and γ_a can be neglected at relatively high temperatures.

It should be noted that the friction coefficient in Eq. (2) depends on an ion charge Z^* that is determined by an average atom model. This introduces a decidedly non-ab initio element to QLMD simulations. In fact, the DT mixture is fully ionized in most of the regime shown in Fig. 1 below. Even if it is only partially ionized, the error introduced in Z^* can be neglected because the friction coefficient is assumed to be in an appropriate range such that the ion–electron collisions and the dynamical behavior are described correctly.⁴⁶ 

The friction coefficient and the Gaussian random force are connected by the fluctuation–dissipation theorem

where dt is the molecular dynamics time step and the random forces are generated as $〈NI2〉=6γMIkBT/dt$⁠.

We use the Verlet algorithm to integrate the Langevin equation (1),

and the ion velocity at time t + dt is

We compare the EOS obtained from QLMD with the AAMD and OFMD simulations at typical density–temperature state points. In AAMD, the AA model is used to solve for the electron density, and then the modified temperature- and density-dependent Gordon–Kim (GK) theory is employed to obtain the ion–ion pair potential based on the electron density.^29,54 Finally, classical MD simulations are carried out for the ion motions. Specifically, we obtain the electron density by using a modified AA model to include the temperature and density effects on the electron distributions in a statistical manner.²⁷ The influence of the plasma environment on the atom is assumed to have spherical symmetry on average, and the occupation number of electrons on the orbitals of such a pseudo-atom is averaged over the possible ionic charge states. The electron orbitals of the ions are solved via the fully relativistic self-consistent-field Dirac equation

where P_nκ(r) and Q_nκ(r) are respectively the large and small components of the wave function. V(r) is the self-consistent potential, which consists of static, exchange, and correlation potentials. Because the thermal fluctuations of ions in a plasma produce dynamic energy level broadening of the ions, Gaussian functions centered at the corresponding electron orbital energies are introduced into the Fermi–Dirac distribution of electrons. With this approach, the instability of the pressure-induced electronic ionization with density can be avoided in a natural manner.²⁷ 

In the GK theory,^29,54 the total energy of a system includes the electrostatic Coulomb potential energy, the exchange potential energy, the correlation energy, and the kinetic energy. We construct a two-atom system, and the ion–ion pair interaction potential is then given by the difference in total energy between the coupled two-atom system and the two isolated-atom systems. It should be noted that at high temperatures, ionic many-body correlation effects will be very weak and can be neglected, and thus a pair potential is accurate enough to describe the ion correlations. However, at low temperatures, many-body correlations become important, and the pair-potential-based AAMD will deviate from ab initio methods.

The main difference between traditional QMD and OFMD is that the driving forces of the ions are obtained from two different DFT approaches: KSDFT and orbital-free DFT (OFDFT), respectively. In the framework of finite-temperature DFT,³³ the electron free energy is obtained by minimizing the grand canonical potential with respect to the electron density n(r). The grand canonical potential has the form⁴³ 

where v(r) is the external potential acting on the electrons corresponding to the density n, and μ is the chemical potential. The free energy functional F[n] is composed of the noninteracting free energy F_s[n], the classical Coulomb repulsion energy (i.e., Hartree energy) F_H[n], and the exchange-correlation free energy F_xc[n],

In conventional KSDFT, a sophisticated scheme exploits the one-electron orbitals of the noninteracting system to construct the electron density of the real system and thereby the total free energy. The advantage of KSDFT is that the noninteracting free energy functionals F_s[n] can be constructed exactly from the one-electron orbitals and electron Fermi–Dirac occupations, thereby giving an explicit Euler equation once a suitable approximate F_xc is provided.

In contrast to conventional KSDFT, in OFDFT, the noninteracting functionals T_s[n] and S_s[n] are formulated directly in terms of the electron density rather than the KS orbitals. Minimization of the grand canonical potential in Eq. (7) with respect to the electron density n(r) then gives the Euler–Lagrange equation

The computational cost of solving this equation scales linearly with the system size and is essentially independent of temperature. The accuracy of OFMD is largely determined by the quality of the noninteracting free energy functional.

We performed extensive calculations for the EOS of a DT mixture over wide ranges of temperature and density for ICF applications. The density–temperature state points chosen in this work are shown in Fig. 1. The density ranges from 0.1 g/cm³ to 2000 g/cm³ and the temperature from 500 K to 2000 eV. Regarding density and temperature, we can use two parameters, namely, the ion coupling parameter Γ = Z^*2/(k_BTa) and the electron degeneracy parameter θ = T/T_F, to define states of matter,⁵⁵ where Z^* is the average degree of ionization, T is the temperature, k_B is Boltzmann’s constant, $a=[3/(4πni)]1/3$ is the mean ion sphere radius, $TF=(3π2ne)2/3/2$ is the Fermi temperature, n_e is the electron number density, and n_i is the ion number density. When the values of θ and Γ are close to 1, matter is in a partially degenerate and moderately coupled state. When θ ≪ 1, matter is strongly degenerate; conversely, it is weakly degenerate when θ ≫ 1. Matter is strongly or weakly coupled when the coupling parameter Γ ≫ 1 or Γ ≪ 1, respectively. As shown in Fig. 1, the coupling parameter Γ corresponding to the state points in this work is greater than 0.1, and most of the state points have Γ ∼ 1 or ≫1. Thus, the DT ions in the ICF process are in states ranging from moderately to strongly coupled. The degeneracy parameter θ corresponding to most of the state points is between 0.01 and 1. Thus, the electrons of DT are partially in strongly degenerate states.

We performed QLMD simulations using our locally modified version of the Quantum-ESPRESSO package.⁵⁶ The generalized gradient approximation in the Perdew–Burke–Ernzerhof parametrization⁵⁷ was used to treat the electron exchange-correlation functional. A norm-conserving pseudopotential was used in low-density conditions and a Coulomb pseudopotential with a cutoff radius of 0.005 a.u. was used in high-density conditions. The plane-wave cutoff energy was from 100 Ry to 1000 Ry, depending on the temperature and density. In the finite-temperature DFT framework, electrons occupy orbitals according to the Fermi–Dirac distribution. We included sufficient energy bands to ensure that the highest occupied band energy was higher than the chemical potential by at least 10k_BT. Owing to the high computational efficiency of QLMD simulations resulting from the large self-consistent field tolerance, we can extend the QLMD simulation to extremely high temperatures at affordable computational cost. We used a supercell including 128–432 atoms, depending on the density. The mixing ratio of D and T atoms was 1:1. The Γ point was used to sample the Brillouin zone in the MD simulations. Convergence test calculations with denser k-point grids and larger supercells did not show any significant variations in the EOS data. The relative error was ∼1% for pressure and 5 meV/atom for energy. For thermodynamic statistics, 2000–5000 steps after thermalization, with time steps of 0.02 fs–1 fs, were used.

In the OFMD calculations, the finite-temperature Thomas–Fermi noninteracting free energy functional with the von Weizsäcker density gradient correction (TFVW)^58,59 was used. The Perdew–Burke–Ernzerhof parameterized generalized gradient approximation functional was adopted for the electronic exchange-correlation interaction. A local pseudopotential was employed in all OFMD calculations.⁴³ The numerical grid for real-space integrations was set to 96 × 96 × 96 to ensure convergence of the free energy and pressure. All the OFMD calculations were performed with our locally modified version of PROFESS.⁶⁰ 

A well-known ab initio EOS table of deuterium for ICF applications was derived by Hu et al.⁶ using PIMC simulations. In Table I, we compare the QLMD results with the PIMC data. We select ten density–temperature points for comparison. It can be seen that the QLMD results are in good agreement with the PIMC data. The pressures obtained from QLMD simulations are slightly lower than those from PIMC simulations at low temperatures. At high temperatures, however, the opposite trend is seen. It was recently demonstrated that the remarkable agreement of QMD simulations with the experimental first-shock Hugoniot of deuterium arises from a cancellation of errors in the DFT model, whereas many-body methods like CEIMC can introduce non-negligible and additive errors into the evaluation of the Hugoniot curve.⁴¹ This reminds us that in order to obtain more accurate EOS data, it is necessary to continuously reduce the size of the approximations in either the DFT model or the PIMC simulations.

We also make a direct comparison of pressures obtained from these methods to show the applicability of AAMD and OFMD. We can see from Fig. 2 that at densities of 10 g/cm³ and 100 g/cm³, the pressures obtained from AAMD and OFMD simulations are in good agreement with those from QLMD in the temperature range from 1 eV to 300 eV, although the AAMD result is slightly higher than that from the QLMD calculation at 1 eV. In contrast to what is found at the densities of 10 g/cm³ and 100 g/cm³, when the density is as low as 4.3 g/cm³ and 1 g/cm³, the pressures exhibit large differences with decreasing temperature. The pressures from AAMD are remarkably larger than those from QLMD at temperatures below 10 eV. We note that the calculation of the electronic structure in AAMD employs a statistical single-atom model for ions in the plasma environment, and it neglects the correlation effects of ions, which play significant roles in the warm dense regime.³⁰ Therefore, AAMD gives inaccurate EOS data in the strong-coupling regime. From the comparisons between QLMD and OFMD results, we can see that OFMD performs better than AAMD. However, the pressures obtained from OFMD at 1 g/cm³ and 4.3 g/cm³ are higher than those from QLMD, especially at low temperatures (<1 eV). This is a consequence of the inability of OFMD to provide satisfactory descriptions of the shell structure of bound electrons and the chemical bond, and therefore OFMD becomes invalid at relatively low temperatures. Moreover, it should be noted that the pressure from OFMD is strongly dependent on the choice of noninteracting free energy functional.^61,62 The accuracy of the noninteracting free energy functional in OFMD simulations plays a critical role in calculations of the thermodynamic properties of matter under extreme conditions.

To gain further insight into the differences between AAMD, OFMD, and QLMD, we compare the radial distribution function (RDF) between these simulations for a DT mixture at 1 g/cm³. Temperatures of 1000 K, 1 eV, and 5 eV are chosen for comparison. The results are presented in Fig. 3. We can see that at both 1000 K and 1 eV, the RDFs of DD, DT, and TT exhibit significant differences between AAMD, OFMD, and QLMD. When the temperature is as low as 1000 K, the RDFs obtained from QLMD have distinct peaks at 0.740 Å, 0.735 Å, and 0.725 Å for DD, TT, and DT, respectively. This means that there are a large number of molecules in QLMD simulations. By contrast, there are no molecular peaks of the RDFs in the OFMD simulations, which indicates that the DT mixture is in a dissociated atomic state in these simulations. The RDFs obtained from AAMD are similar to those from OFMD, although the first peaks are more structured in the case of AAMD. When the temperature is increased to 1 eV, there are still remarkable shoulders at about 0.9 Å in the QLMD simulations, indicating that there exist somewhat softened molecular structures. At 1 eV, the RDFs obtained both from AAMD and from OFMD show the distinct characteristics of atomic states, which is due to the fact that OFMD and AAMD cannot satisfactorily describe the bond formation in hydrogen molecules at low temperatures. When the temperature is increased to 5 eV, the RDFs obtained from three methods exhibit similar behavior, although there are still slight differences at the rising edge. At 5 eV, molecules have dissociated completely, and thus the spatial distributions of ions are in better agreement. From these comparisons, we can conclude that AAMD is accurate for hot dense DT mixtures and OFMD can be applied for calculating EOS data of DT mixtures over wider ranges of density and temperature than AAMD. However, neither of these two methods can provide accurate bonding information for DT mixtures at relatively low temperatures. In comparison with AAMD and OFMD, QLMD which is based on an accurate description of electronic structure and an efficient MD algorithm, can provide accurate EOS data for DT mixtures from low-temperature condensed states to the hot dense regime.

Table II shows the pressure and internal energy of DT mixtures obtained from QLMD simulations. The density ranges from 0.1 g/cm³ to 2000 g/cm³ and the temperature from 500 K to 2000 eV. Here, the total internal energy is obtained from

where F is the free energy of the simulation system obtained from the finite-temperature DFT self-consistent-field iterations, S is the entropy, T is the temperature, and E_kin is the kinetic energy of ions. The pressure is calculated from

where n is the ion number density, nk_BT is the ideal kinetic contribution of ions, and P_DFT is the interaction contribution calculated from DFT self-consistent-field iterations, which includes contributions from the kinetic energy of electrons, the ion–electron interaction, the Hartree interaction of electrons, and the electronic exchange-correlation interaction.

We should stress that the KSDFT calculations are still computationally expensive at high temperatures, even though QLMD lowers the convergence criteria of the self-consistent-field calculations. Therefore, 248 density–temperature state points are adopted to calculate the EOS of a DT mixture, as shown in Table II. These EOS data should be interpolated with a denser density–temperature grid set when applied to ICF hydrodynamic simulations. The interpolation approach can affect the accuracy of interpolated EOS data,⁶³ and therefore a proper interpolation method needs to be employed. However, this is outside the scope of this work. In addition, when the density of matter is extremely high and the temperature is relatively low, nuclear quantum effects become significant for atomic structures, transport properties, and thermodynamic properties, especially in the case of light elements such as hydrogen.^64,65 We can use a parameter α, which is defined as the ratio of the ionic thermal de Broglie wavelength to the mean distance between ions, to measure the degree to which the ions exhibit a quantum nature. For all the densities and temperatures considered here, α is less than 0.3. In terms of our previous investigations,⁶⁶ nuclear quantum effects play notable roles in static structures and thermodynamic properties only if α > 0.3. Therefore, the quantum nature of the DT ions can be neglected in this work.

We have performed extensive QLMD simulations to obtain EOS data for DT mixtures over wide ranges of density (from 0.1 g/cm³ to 2000 g/cm³) and temperature (from 500 K to 2000 eV) relevant to ICF implosions. Comparisons with AAMD and OFMD simulations reveal significant discrepancies at relatively low temperatures, where the strong ionic coupling plays a remarkable role in determining the EOS of DT mixtures. The DFT-based simulation methods provide more reliable EOS data than previous semiclassical methods. In the future, we should pay special attention to basic physical issues such as electronic many-body effects and nonlocal interactions⁶⁵ to meet the requirements for higher-precision EOS data not only for ICF applications, but also for planetary science and astrophysics.

This work was supported by the Science Challenge Project under Grant No. TZ2016001, the National Key R&D Program of China under Grant No. 2017YFA0403200, the National Natural Science Foundation of China under Grant Nos. 11874424 and 11774429, and the NSAF under Grant No. U1830206. All calculations were carried out at the Research Center of Supercomputing Applications at NUDT.