We study nonideal mixing effects in the regime of warm dense matter (WDM) by computing the shock Hugoniot curves of BN, MgO, and MgSiO3. First, we derive these curves from the equations of state (EOS) of the fully interacting systems, which were obtained using a combination of path integral Monte Carlo calculations at high temperature and density functional molecular dynamics simulations at lower temperatures. We then use the ideal mixing approximation at constant pressure and temperature to rederive these Hugoniot curves from the EOS tables of the individual elements. We find that the linear mixing approximation works remarkably well at temperatures above ∼2 × 105 K, where the shock compression ratio exceeds ∼3.2. The shape of the Hugoniot curve of each compound is well reproduced. Regions of increased shock compression, which emerge because of the ionization of L and K shell electrons, are well represented, and the maximum compression ratio of the Hugoniot curves is reproduced with high precision. Some deviations are seen near the onset of the L shell ionization regime, where ionization equilibrium in the fully interacting system cannot be well reproduced by the ideal mixing approximation. This approximation also breaks down at lower temperatures, where chemical bonds play an increasingly important role. However, the results imply that the equilibrium properties of binary and ternary mixtures in the regime of WDM can be derived from the EOS tables of the individual elements. This significantly simplifies the characterization of binary and ternary mixtures in the WDM and plasma phases, which otherwise requires large numbers of more computationally expensive first-principles computer simulations.
The physical properties of hot, dense plasmas have been studied with experimental and theoretical techniques for decades1 because their behavior is important for a number of energy technologies, including inertial confinement fusion (ICF).2–7 On the path to fusion, the sample material typically passes through the regime of warm dense matter (WDM), which encompasses matter at solid-state densities and elevated temperatures of 104 K–107 K. This regime is particularly difficult to describe with theoretical methods because the densities are too high and interaction effects are too strong for typical plasma theory to be applicable8,9 or for Saha ionization models to work properly.10 On the other hand, the temperatures are too high and the fraction of excited electrons too large for conventional condensed matter theory to be applicable. The temperature is also not high enough for screening effects to become the dominant type of interaction, and thus, Debye plasma models11 do not work well. All particles are strongly interacting, which renders the system nonideal. There is no small parameter that would allow for analytical descriptions to be appropriate. Chemical bonds still play a role, even though they are typically short-lived. The electrons may be highly excited and partially ionized. Pauli exclusion effects are relevant when the ionization equilibrium is established, which renders the system partially degenerate. A good fraction of the electrons occupy core states because density is orders of magnitude too low for them to form a rigid neutralizing background. In this regard, a one-component plasma model would be a poor description of WDM. Despite these challenges, the development of a rigorous and consistent theoretical framework to describe WDM remains of high importance.
Significant progress has been made with laboratory experiments and first-principles computer simulations. A number of different simulation methods have been advanced.12 These simulation methods enable us to compute the equation of state (EOS) of materials over a wide range of conditions that are also relevant in astrophysics. In the interiors of giant planets, for example, not only hydrogen–helium mixtures but also rocky materials are exposed to tens of megabars and ∼104 K.13–22 Accurate EOSs are needed to characterize their interior structure and evolution.16,23 The discovery of thousands of exoplanets with ground-based observations and space missions as well as the unexpected diversity in their masses and radii24,25 considerably broadened the range of conditions and materials of interest.26–28
In the interior of stars, matter is exposed to a wide range of temperatures ∼104 K−108 K. The most detailed constraints on the interior conditions come from the measurements of normal mode oscillations of our Sun.29–31 Such astero-seismological observations are also employed to constrain the interiors of distant stars.32 The interpretation of these observations would not be possible without a comparable development of laboratory experiments to probe such extreme conditions,33 which employ a variety of high-velocity impacts,34,35 lasers,3,4,7,36–40 and magnetic compression techniques.41,42 The goal of this article is to support these activities by providing theoretical methods to predict the state of WDM with computer simulations and to derive the EOS of materials over a wide range of temperature–pressure conditions. This will aid in the interpretation of current experiments or help with their design in the future. Therefore, any material and condition that can be probed with current experimental facilities is of potential interest.
The range of pressure–temperature conditions of interest is very large and so is the space of possible chemical compositions. The challenge of dealing with this huge number of materials and conditions is not unique to the field of WDM.43 For binary and ternary mixtures, the number of relevant conditions scales as and , where Nρ, NT, NE, and Nmix are the numbers of densities, temperatures, elements, and mixing ratios of interest, respectively. The resulting numbers are often so large that an exhaustive coverage is impractical not only for laboratory experiments but also for computer simulations, although exceptions exist.44 To simplify the computation of WDM, we investigate the validity of the linear mixing approximation at high pressures and temperatures in this article.
The linear mixing approximation is a widely used approach to obtain the equation of state (EOS) of materials from the individual components if the information about the fully interacting mixture is lacking. The ideal mixing rule is often used to study gaseous mixtures of simple elements, such as hydrogen, helium, carbon, and oxygen, to understand the atmospheres and interior of gas giant planets, where complex mixtures at different concentrations emerge and whose physical properties are unknown.45–47 The accuracy of this approach depends largely on the thermodynamic conditions under which it is applied.48 A wide variety of systems have been explored under the assumption of ideal mixing.49–56 The hypothesis works remarkably well for hydrocarbon mixtures in the warm dense matter regime,57 water–hydrogen mixtures,58 as well as in hydrogen–helium mixtures enriched in heavier elements.17,33,59,60 The range of validity of the linear mixing hypothesis has been explored in great detail for hydrogen–helium mixtures.61 As an example of a nonlinear mixing effect, Vorberger et al.62 showed that the presence of helium stabilizes the hydrogen molecules in the mixture and thus moves their molecular-to-metallic transition to higher pressures.
In this article, we employ two first-principles computer simulation methods, path integral Monte Carlo (PIMC) calculations and density functional molecular dynamics (DFT-MD) simulations, to study nonideal mixing effects in the regime of WDM. With this approach, we have been able to produce several EOSs in previous years, which have supported the laboratory experiments. In particular, our simulation results for warm dense carbon63 are currently being used to design NIF targets.
In this article, we will show with our first-principles computer simulations that the compositional dependence of the EOS is manageable in the regime of WDM. We will demonstrate that the shock Hugoniot curves of various mixtures can be derived with good accuracy for temperatures above ∼2 × 105 K by invoking the ideal mixing approximation at constant pressure and temperature. This means chemical bonds between species no longer play an important role at these temperatures. It is still surprising that the properties of hot, dense MgSiO3 can be derived from those of the pure substances because one essentially assumes that, e.g., a Mg ion in a dense MgSiO3 environment behaves very similarly to the one that is surrounded by other Mg ions under the same P–T conditions. Without verification, there is no guarantee that the degree of ionization will be similar in the two systems. The goal of this article is to investigate these similarities and to characterize the nonlinear mixing effects quantitatively for the three representative WDM materials BN, MgO, and MgSiO3.
II. METHODS AND ASSUMPTIONS
We derived the equation of state of every material under consideration, Mg, Si, O, B, N, MgO, MgSiO3, and BN, by performing a series of first-principles computer simulations that employed PIMC simulations at high temperature and standard Kohn–Sham DFT-MD calculations at low temperature. We describe these methods in Secs. II A and II B.
A. PIMC simulations
Path integral Monte Carlo (PIMC) methods have gained considerable interest as a state-of-the-art, stochastic first-principles technique to compute the properties of interacting quantum systems at finite temperature. This formalism results in a highly parallel implementation and an accurate description of the properties of materials at high temperature where the electrons are excited to a significant degree.63,68,70–72 The application of the PIMC method to light elements from hydrogen through neon71 has been possible due to the development of free-particle nodes,73,74 while simulations of heavier elements relied on the advancement of Hartree–Fock nodes.75 The latter approach enables one to efficiently incorporate localized electronic states into the nodal structure, which extends the applicability of the path integral formalism to heavier elements and lower temperatures.76,77 Furthermore, PIMC treats all electrons explicitly and avoids the use of pseudopotentials. The PIMC simulation time scales as 1/T, proportional to the length of the paths, which is efficient under high-temperature conditions, where most electrons including the K shell are excited. Predictions from PIMC simulations at intermediate temperatures have been shown to be in good agreement with predictions from density functional theory molecular dynamics (DFT-MD) simulations.14,44
The fundamental techniques for the PIMC simulations of bosonic systems were developed in Ref. 78 and reviewed in Ref. 73. Subsequently, the algorithm was extended to fermionic systems by introducing the restricted path approach.74,79,80 The first results of this simulation method were reported in the seminal work on liquid 3He80 and dense hydrogen.81 In subsequent articles, this method was applied to study hydrogen,82–87 helium,55,70,88 hydrogen–helium mixtures,48 and one-component plasmas.89–91 In recent years, the PIMC method was extended to simulate the plasmas of various first-row elements,40,44,57,63,92,93 and with the development of Hartree–Fock nodes, the simulations with heavier nuclei up to silicon became possible.72,75–77
The PIMC method is based on the thermal density matrix of a quantum system, , which is expressed as a product of higher-temperature matrices by means of the identity , where M is an integer and τ ≡ β/M represents the time step of a path integral in imaginary time. The path integral emerges when the operator is evaluated in real space,
The sum includes all permutations, , of N identical fermions in order to project out the antisymmetric states. For sufficiently small time steps, τ, all many-body correlation effects vanish, and the action, S[Rt], can be computed by solving a series of two-particle problems.78,94,95 The advantage of this approach is that all many-body quantum correlations are recovered through the integration over paths. The integration also enables one to compute the quantum mechanical expectation values of thermodynamic observables, such as the kinetic and potential energies, pressure, pair correlation functions, and the momentum distribution.73,96 Most practical implementations of the path integral techniques rely on Monte Carlo sampling techniques because the integral has D × N × M dimensions in addition to sum over permutations. (D is the number of spatial dimensions.) The method becomes increasingly efficient at high temperature because the length of the paths scales like 1/T. In the limit of low temperature, where few electronic excitations are present, the PIMC method becomes computationally demanding and the Monte Carlo sampling can become inefficient. Still, the PIMC method avoids any exchange–correlation approximation and the calculation of single-particle eigenstates, which are embedded in all standard Kohn–Sham DFT calculations.
The only uncontrolled approximation within fermionic PIMC calculations is the use of the fixed-node approximation, which restricts the paths in order to avoid the well-known fermion sign problem.74,79,80 Addressing this problem in PIMC is crucial, as it causes large fluctuations in computed averages due to the cancellation of positive and negative permutations in Eq. (1). We solve the sign problem approximately by restricting the paths to stay within nodes of a trial density matrix that we obtain from a Slater determinant of single-particle density matrices,
which combined free and bound electronic states,75,77
The first sum includes all plane waves, Ψk, while the second represents n bound states Ψs with energy Es that are localized around all atoms I. Predictions from various slightly differing forms of this approach have been compared in Ref. 76.
The PIMC simulations were performed with the CUPID code84 using periodic boundary conditions. For pure O, Mg, Si, and N systems, we considered simulation cells with 8 nuclei and 64, 96, 112, and 56 electrons, respectively. For the simulations of boron, slightly larger cells with 30 nuclei and 150 electrons are used. For MgO, we considered 80 electrons, 4 Mg, and 4 O nuclei. For BN, we considered 144 electrons, 12 B, and 12 N nuclei, while our MgSiO3 simulations consisted of 3 Mg, 3 Si, and 9 oxygen nuclei as well as 144 electrons. A finite size effect was discussed when we first reported these EOS calculations. A detailed finite-size study is provided in the supplementary material of Ref. 71 that showed that finite size effects are relatively small because the most important changes in the energy and pressure are caused by the ionization of various electronic states. While the ionization equilibrium depends on the thermodynamic conditions of the plasma, it does not require large simulation cells to capture these effects.
We enforced the nodal constraint in small steps of imaginary time of τ = 1/8192 Ha−1, while the pair density matrices95 were evaluated in steps of 1/1024 Ha−1. This results in using between 2560 and 5 time slices for the temperature range that was studied with PIMC simulations. These choices converged the internal energy per atom to better than 1%. We have shown that the associated error is small for relevant systems at sufficiently high temperatures.97
B. DFT-MD simulations
Kohn–Sham DFT98,99 is a first-principles simulation method that determines the ground state of quantum systems with high efficiency and reasonable accuracy, which has gained considerable use in computational materials science. The introduction of the Mermin scheme100 enabled the inclusion of excited electronic states, which extended the applicability range of the DFT method to higher temperatures. The combination of this method with molecular dynamics has been widely applied to compute the EOS of condensed matter, warm dense matter (WDM), and some dense plasmas.40,101–103 Unless the number of partially occupied orbitals is impractically large, DFT is typically the most suitable computational method to derive the EOS because it accounts for electronic shell and bonding effects. The main source of uncertainty in DFT is the use of an approximate exchange–correlation (XC) functional. The errors resulting from the XC functional often cancel between different thermodynamic conditions. Furthermore, this error may only be a small fraction of the internal energy, which besides pressure is the most relevant quantity for the EOS and the derivation of the shock Hugoniot curve.104 However, the range of validity of this assumption in the WDM regime remains to be verified for different classes of materials through the comparison with laboratory experiments and other computational technique such as PIMC simulations.
With the VASP code,105 we performed simulations from 104 up to 2 × 106 K. We employed a Nosé thermostat106 to keep the temperature constant. As illustrated in Fig. 1, we explored densities from 6.89 g cm−3 to 51.67 g cm−3 for Mg, 1 g cm−3 to 100 g cm−3 for O, 2.3 g cm−3 to 18.6 g cm−3 for Si, 0.35 g cm−3 to 71 g cm−3 for MgO, and 0.321 g cm−3 to 64.2 g cm−3 for MgSiO3. We used cubic simulation cells with periodic boundary conditions, and to improve efficiency, we used a smaller number of atoms at the highest temperatures. As shown in our previous work,57,77,107,108 this is not detrimental to the accuracy of the EOS data at high temperatures. More details of the simulations for BN, B, and N can be found in our previous publications.40,68,92 We employed projector augmented wave (PAW)109 pseudopotentials with a 1 s2 frozen core for all elements, and mostly used the Perdew–Burke–Ernzerhof (PBE) functional110 at the majority of conditions to describe the exchange–correlation effects. For some materials and a small number of conditions, we had to switch to the local density approximation (see details in Refs. 108 and 111) but have very good agreement between the obtained results when we switched functionals. The time step was adapted according to the corresponding temperature, and a large energy cutoff was used for the plane wave basis set.
C. Shock Hugoniot curves
The shock Hugoniot curves of many materials have been measured up to megabar, and in some cases gigabar, pressures.101,112–114 Even under extreme conditions,44,68,102,103,108,115 predictions from first-principles simulations and experiments have been shown to be in good agreement.
The EOS can be used to predict the thermodynamic conditions that are reached when a material is subjected to dynamical shock compression. Assuming that the materials reached thermodynamic equilibrium during the experiments, the measured shock and particle velocities can be converted into pressure, density, and energy through the Rankine–Hugoniot equations.116–118 The energy conservation equation
is particularly convenient to derive the shock Hugoniot curve with theoretical methods. Here, E0, V0, and P0 represent the initial conditions of energy, volume, and pressure, respectively. E, V, and P are the final conditions after the material behind the shock front has reached an equilibrium state. We solve Eq. (4) for T and V by interpolating E(V, T) and P(V, T) in our EOS tables. Most simply, one solves for V at given T because there is only one solution.
D. Linear mixing
The linear mixing approximation at constant pressure and temperature is the most common form in astrophysics,64 but it is often used also in plasma physics.119 Still, variations and alternate mixing rules have been invoked.120 In Ref. 59, the linear mixing approximation was employed to perturb the helium fraction in an interacting H–He EOS. When the mixture of carbon, oxygen, and neon nuclei are studied under conditions in White Dwarf interiors, one would want to mix the individual EOSs at constant temperature and nuclear charge density for the following reason. The density in white dwarfs is sufficiently high for the electrons to decouple from the motion of the nuclei and to form a rigid neutralizing background. That background, however, provide the dominant contribution to the pressure and that is a function of the electron density, which is equal to the nuclear charge density. Therefore, mixing at constant P translates into mixing at constant nuclear charge.
Plasmas have also been studied with two-temperature models119 that treat nuclei and electrons as two independent thermodynamic ensembles with differing temperatures because the inter-species thermalization is delayed by the difference in mass. A number of other approaches to study mixtures120–124 have been proposed with the goal of facilitating large scale hydrodynamic simulations. Some of these approaches have been verified by orbital-free molecular dynamics.125
In this article, however, we use the simplest form of the linear mixing approximation. For a mixture of species A and B, one neglects all inter-species interactions and, for given P and T, one assumes that the volume of the mixture is given by Vmix(P, T) = VA(P, T) + VB(P, T). The mass density, ρmix, is given
where xA and xB are the mass fractions of each species in the mixture. The internal energy is then given by
where all three energy terms are normalized per unit mass. When theoretical and computational results are employed, it is often more convenient to normalize all quantities by formula unit (FU). Let us assume that V1(P, T) and V2(P, T) are the volumes per formula unit of species 1 and 2. N1 and N2 specify how many formula units of species 1 and 2 are contained in one unit of the mixture. For given P and T, the volume, mass, and internal energy of one mixture unit are obtained from
The mass density of the mixture is given by ρmix = mmix/Vmix.
The linear mixing approximation only provides reasonable results for the mixture at elevated temperatures where chemical bonds do not affect the EOS significantly. For this reason, we always use the E0, P0, and V0 of the fully interacting system when we solve Hugoniot equation (4). E and V, however, can be approximated by Emix and Vmix. To solve Eq. (4) for a mixture, we assume a temperature and a value for Vmix. Then, we determine the pressure that matches Vmix and we derive the corresponding Emix. We iterate over Vmix to find a solution of Eq. (4).
III. RESULTS AND DISCUSSION
In Figs. 2–4, we compare the shock Hugoniot curves that were computed with the fully interacting internal energy and pressure for BN,68 MgO,108 and MgSiO3115,127 with those derived from the linear mixing approximation using the elemental EOS tables. The agreement between the pairs of curves is remarkably good in pressure–density and temperature–density spaces. For temperatures above ∼2 × 105 K corresponding to shock compression ratios above ∼3.2, the shape of the Hugoniot curve is very well reproduced by the linear mixing approximation. This includes the regimes of K and L shell ionization. The compression maximum is also well reproduced. The linear mixtures of three elements show the same level of agreement with the fully interacting Hugoniot curves as linear mixtures of two elements. We do see some deviations for MgSiO3 at 2 × 106 where the Hugoniot curves transition between the K and L shell ionization regimes. Under these conditions, the linear mixing approximation does not accurately capture the ionization equilibrium of the interacting system.
For comparison, we also show the Hugoniot curves for the individual elements in Figs. 2–4. The differences in the Hugoniots for the individual elements and that of the compounds are primarily due to differences in the initial densities. In Fig. 2, we also show the results from laser shock experiments68 that reached up to a pressure of 2643 GPa and a compression ratio of 2.66. This is not sufficiently high for the linear mixing approximation to work well. While in pressure-compression space, the linear mixing and the fully interacting Hugoniot curves both agree the experimental data, but the temperature-compression graph of Fig. 2 reveals that the shock temperatures are underestimated for compression ratios below 3.2 if the linear mixing approximation is invoked. We see the same trend in Fig. 3 where we compare our theoretical predictions for MgO with the experimental results from Ref. 67 that reached up pressures of 2303 GPa and compression ratios of 2.68. While in pressure-compression space, the predictions from the linear mixing approximations appear to be reasonable, the shock temperature is underestimated for compression ratios smaller than 3.2. Finally in Fig. 4, we compare with the shock experiments in Ref. 69 that reached up to 1426 GPa and a compression ratio of 2.26. While these results are in good agreement with fully interacting DFT-MD simulations,69 the shock temperatures are underestimated if the linear mixing approximation is employed. The reason for this discrepancy is that, for given pressure and temperature, the linear mixing approximation underestimates the density and the internal energy for T ≦ 2 × 105 K, as we confirm in the following analysis:
The shock Hugoniot curves can only be reproduced well by the linear mixing approximation as long as the volume and internal energy, which enter into Eq. (4), are reasonably accurate. In Fig. 5, we plot the deviation in ρmix and the interacting ρ for two temperatures as a function of density. The error in density is less than 1% for all three materials with the exception of a B + N mixture, in which case the error reaches 2% at low and high densities.
In Fig. 6, we plot the linear mixing error in the internal energy that we normalized by dividing by nuclear kinetic energy . The deviations are 0.1 or less for all three materials and conditions under consideration.
In Fig. 7, we plot the linear mixing errors as a function of temperature for three relevant density sets equal to 4.5 times the initial shock density, ρ0. The deviations in density are 2% or less except at 1.5 × 106 K where we switch between PIMC and DFT-MD EOS computations, in which case the deviations reach 4%. This, however, does not reflect any insufficiency in the linear mixing approximation, but the underlying EOS tables are imperfect. When we study the linear mixing error in the internal energy, we also find a discrepancy for MgSiO3 at 1.5 × 106 K.
For temperature below 2 × 105 K, the errors in density and internal energy of the linear mixing approximation increase with decreasing temperature because chemical bonds and interactions between different species play an increasingly important role. Chemical bonds lower the internal energy and pressure. Since bonding effects are absent from the linear mixing approximation, it overestimates the internal energy and density for given pressure and temperature, which explains the trends at lower temperature in Fig. 7. Still, already for T ≧ 2 × 105 K, the linear mixing approximation works well.
Finally, we performed three tests how much of a change in pressure and in internal energy is needed to shift the Hugoniot curve in temperature-compression and pressure-compression spaces. First, we determined what fractional change in pressure is needed to reduce the maximum compression ratio on the three principal Hugoniot curves by 10%. The results of our calculations in Table I show that a 3% increase in pressure would trigger such a shift. Alternatively, such a reduction in compression ratio can be introduced by lowering the internal energy by 0.21 ⋯ 0.34 . The magnitude of both corrections is much larger than the linear mixing error that we reported in Figs. 5–7, which explains why we were able to reproduce the compression maxima of the Hugoniot curves very well with the linear mixing approximation.
|Material .||ρ0 (g cm−3) .||.||.||.||.||.||.|
|Material .||ρ0 (g cm−3) .||.||.||.||.||.||.|
In Table I, we report the results from two more tests at lower temperatures and pressures. Starting the point at 3.5-fold compression on the Hugoniot curve, we asked what fractional pressure change and what energy correction in units of would be needed to move ρ/ρ0 = 3.5 point up in temperature by 5% or up in pressure by 5%. The required pressure and energy corrections are reported in columns 5–8 of Table I. Energy changes between −0.097 and are needed to change the temperature and pressure on the Hugoniot curve by 5%. These changes are comparable in magnitude to the linear mixing errors we have reported in Figs. 6 and 7. Hence, at 3.5-fold compression, the accuracy of the linear mixing approximation is about 5%.
In Table I, we also show that fraction pressure changes between 0.009 and 0.019 are needed to move ρ/ρ0 = 3.5 point by 5% in pressure and temperature. These values are larger than the linear mixing errors in density at 5 × 105 K that we show in Fig. 5 that amount to less than 1% at this temperature. This suggest that the changes in the internal energy are slightly more difficult to reproduce with the linear mixing approximation than the pressure.
We have validated the linear mixing approximation across a wide range of temperature and pressure conditions for MgO, MgSiO3, and BN plasmas. Under this approximation, accurate shock Hugoniot curves can be obtained for temperatures of T ≧ 2 × 105 K and compression ratios of ρ/ρ0 ≧ 3.2, correctly predicting the maximal compression ratio and the K- and L-shell ionization regimes. This will greatly simplify the computations for the regime of WDM and may even help reduce the number of experiments. This conclusion is further supported by the first-principles calculations for CH,57 which reported that the maximal compression ratio on the Hugoniot curve can be derived with an accuracy of 1% by combining the EOSs of elemental carbon and hydrogen. Similarly, we are able to reproduce the maximum compression ratio of B4C66 with an accuracy of 0.4% if we mix the EOS of boron and carbon. Soubiran and Militzer58 determined that the hydrogen and water form an ideal mixture under conditions of ice giant envelopes. On the other hand, mixtures of hydrogen, helium, and heavier elements in giant planet envelopes could only be accurately represented as a linear mixture after the volumes of the heavier species had been adjusted to match the results from fully interacting DFT-MD simulations.17 The validity of the linear mixing approximation depends on the pressure and temperature conditions. Magyar and Mattsson showed that errors of 10% can be expected for xenon–deuterium mixtures at megabar pressures and 10 000 K120 while we have shown here that at much higher temperatures, the linear mixing approximation works very well. Still, in this article we only investigate how well the volume and internal energies can be derived by combining the EOSs of the various elements at constant pressure and temperature. Future work should thus be directed on understanding to which degree the transport properties of WDM are affected by nonlinear mixing effects.
This work was, in part, supported by the National Science Foundation-Department of Energy (DOE) partnership for plasma science and engineering (Grant No. DE-SC0016248) and by the DOE-National Nuclear Security Administration (Grant No. DE-NA0003842). Part of this work was performed under the auspices of the U.S. DOE by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. Computational support was provided by the Blue Waters sustained-petascale computing project (No. NSF ACI 1640776) and the National Energy Research Scientific Computing Center (NERSC).
H.D.W., D.C.S., and M.M. state that their work was sponsored by an agency of the United States government. Neither the United States government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof, and shall not be used for advertising or product endorsement purposes.
The datasets of all figures are available from the corresponding author upon request.