Lawrence Livermore National Laboratory

R. E. Rudd (15-ERD-052)


We developed a new capability to calculate species-transport properties of asymmetric plasma mixtures, i.e., plasmas with two or more components differing significantly in charge and mass in the presence of strong gradients. Inertial confinement fusion with ablator materials mixing into the deuterium-tritium plasma provides a good example of such gradients. We built the capability on the molecular-dynamics plasma-simulation code ddcMD (developed in the Cimarron project for prediction and measurement of properties of dense, highly collisional, and strongly-coupled plasmas) to deliver a high-fidelity description of plasmas, including those with strong ion–ion coupling for at least one component. We developed the multiscale orbital-free, density functional theory molecular dynamics (MOD-MD) technique that enables simulation of ion dynamics in non-uniform plasmas (such as interfaces) without explicitly calculating the electron dynamics. MOD-MD was implemented in the ddcMD code. We used this capability to simulate diffusive mixing at a fuel-ablator interface in inertial confinement fusion, super-diffusive transport at shock fronts, and mixing at a diamond–metal interface under the conditions of recent warm-dense-matter experiments. We quantified species separation at the inertial confinement fusion fuel-ablator interface and the shocked hydrogenic plasma. This species separation is described to some degree by barodiffusion, but in the shock case, we have shown that the flux is larger than predicted by barodiffusion and non-Fickian in nature. We used MOD-MD to calculate electrical fields at the interfaces that arise due to change in electron density across the interface, and we have shown how these electric fields decay in time as the plasmas mix. We have also developed a fast multipole capability to reduce the computational cost for weakly screened plasma. We calculated the viscosity of plasma mixtures and have shown that it is well described by a kinetics–molecular-dynamics viscosity model suitable for use in hydrodynamic simulations.

Background and Research Objectives

Critical aspects of Lawrence Livermore National Laboratory’s mission (such as inertial confinement fusion, ICF) involve plasma physics and warm dense matter. Plasma is formed when matter is heated sufficiently to be ionized (Spitzer 1962). Warm dense matter is also partially ionized, but at a somewhat lower temperature, so that the transport properties are affected by the potential energy of the ions and the quantum-mechanical degeneracy of the electrons (Whitley et al. 2015). Our research project developed techniques to calculate the transport properties of plasma that are valid in the presence of strong gradients and in the warm-dense-matter regime as the plasma becomes moderately to strongly coupled and is degenerated.

Conventional modes of plasma transport rely on kinetic theory, using techniques that were derived for dilute gases (Chapman and Cowling 1939). They rely on various approximations (such as binary collision integrals and small-scattering-angle approximations) that become problematic as plasma coupling increases and the spatial correlations of the ions matter. We used a different approach based on molecular dynamics (MD) (Frenkel and Smit 2001). Molecular dynamics simulates the motion of a large collection of ions, calculating the forces the ions exert on each other and integrating the equations of motion to determine the trajectories of the ions. The locations of the ions are known at every instant, so it is not necessary to assume that only two ions interact during a collision (the binary-collision approximation). Because the technique works even for strong coupling, MD is well suited for studying transport processes in the warm-dense-matter regime (Sanbonmatsu and Murillo 2001; Graziani et al. 2012). Bonds form at sufficiently cold temperatures and we assume that the systems of interest are ionized and that chemistry can be neglected. There are rare strong collisions that require a time step at sufficiently high temperatures that is much shorter than the mean collision time, and in which molecular dynamics becomes inefficient. There is a broad temperature range between these regimes over which MD with generalized screened Coulomb interactions provides a powerful description of dense plasma and warm dense matter.

One approach in MD is to simulate the motion of electrons, making use of quantum statistical potentials to account for quantum mechanical-diffraction effects (Graziani et al. 2012). Because the electrons are much lighter than the ions, simulation of the electrons requires a very short time step, compared to ion-transport times. We avoided the computational cost of a short time step by treating the electrons with a technique developed at the Laboratory called multiscale orbital-free, density functional theory (DFT) molecular dynamics (MOD-MD) (Stanton et al.). MOD-MD only tracks the motion of the ions, calculating their level of organization and the resulting electron density in the Born approximation that the electrons relax on a shorter time scale than that of the ion motion. The calculations account for electrical fields due to the electrons, including their ability to screen the ions, as described below. The result is a technique that advances with a time step set by ion motion rather than electron motion, so it is possible to simulate significantly longer periods of time.

We applied the MOD-MD technique to several problems in plasma transport: (1) mixing at the fuel-ablator interface in ICF capsules; (2) species separation in shocked plasma; and (3) mixing at the carbon–metal interface in warm dense matter experiments. These topics fall into the category of processes that may surpass conventional hydrodynamics and linear response transport. We also studied related topics such as fast ion dynamics at the interface related to ionic albedo and plasma viscosity, using uniform screened Coulomb MD for the calculations.

Scientific Approach and Accomplishments

The MOD-MD Technique

The MOD-MD model calculates plasma properties at length scales ranging from atomic to inter-ionic to continuum (Stanton et al.). It accounts for spatial and temporal variations of the electron density, resulting in a formulation similar to the screened Coulomb (Yukawa) model with a screening length that depends on the position of the ion. The ion motion is simulated in MD by integrating force (i.e., f = ma) in time (Frenkel and Smit 2001). The forces are determined by the spatially varying screened Coulomb interaction, including a dipole term that accounts for polarized screening clouds that are not spherically symmetric about the ions. The ion charge and the screening length are calculated using a Laplace–Thomas–Fermi solution on a coarse-grained mesh for which the ion charge density and temperature are the inputs. This calculation also gives the coarse-grained electron density and electric field. This framework was systematically derived from the underlying orbital-free DFT. MOD-MD was implemented in the code ddcMD (Glosli et al. 2010) in an algorithm that scales well on the Laboratory's supercomputers. It has been run on up to 65,536 processors (262,144 MPI tasks) with good performance. The configuration of a MOD-MD simulation is shown in Figure 1, in this case simulating mixing at the interface between a plastic ablator and deuterium-tritium (DT) ice fuel.

Figure 1.
Figure 1. Example of multiscale orbital-free density functional theory (DFT) molecular dynamics (MOD-MD). The lower image shows roughly half of the 57.6-million-ion simulation box in the long direction. The two insets show successively enlarged images of the interface. The colors indicate the ion species, as shown in the legend (bottom right). The greenish region is the plastic, the purplish region is deuterium-tritium fuel. The profiles for this snapshot 300 fs into the simulation are plotted in Figure 2 (Stanton et al.).

Diffusive Mixing at the Ablator–Fuel Interface

MOD-MD in the ddcMD code has been used to simulate ion dynamics at a heated interface characteristic of an ICF capsule, as shown in Figure 1 (Stanton et al.). In this series of simulations, the capsule consists of a plastic ablator and a DT ice fuel layer. The heating and resulting diffusive mixing near the interface were simulated on the scales of microns and tens of picoseconds (ps). The electrons in the capsule were heated from initially cryogenic temperatures to 50 eV, either instantaneously or ramped up for 10 ps. The ions also heated at a rate controlled by the electron–ion coupling. The simulations showed separation of the momentum fields among the species and strong hydrogen jetting from the plastic into the fuel region. The interface profiles after 0.3 ps and 10 ps are shown in Figure 2. The electron density was initially very different in the fuel and the ablator, as the higher Z elements in the ablator ionize more strongly than the hydrogenic species. The electric field created at the interface was extremely strong initially but decayed as the interface broadened and pressure equilibrated. Screening cloud polarization was found to reduce the electric field somewhat. In this case, there is little species separation in the DT fuel, so barodiffusion did not play an important role in the fuel. There was some C-H species separation in the plastic. These results were presented at the Laboratory's kinetics workshop and used for continuum code validation (Murillo 2016).

Figure 2.
Figure 2. Mixing at an inertial confinement fusion (ICF) capsule fuel–ablator interface. The concentrations of the ion species (H, D, T, C, and O) are plotted as the electron temperature is ramped to 50 eV over 10 ps. The left panel shows the profiles at 0.3 ps; the right panel, at 10 ps. The gray vertical line indicates the original position of the interface, with plastic at negative positions and deuterium-tritium (DT) fuel at positive positions. The oxygen density has been multiplied by 50 for visibility. Note that the D and T densities remain locked together, whereas the C and H densities separate. There is strong interpenetration of the H into the DT fuel region (Stanton et al.).

Species Separation in Shock Waves

We used MOD-MD in ddcMD to simulate the propagation of a shock wave through a binary plasma mixture of hydrogen isotopes to determine whether the shock propagation obeys the laws of hydrodynamics or if additional physical effects are important (Mackay et al.). This study was motivated by recent work that proposed kinetic effects are important in ICF (Amendt et al. 2010; Bellei et al. 2013; Rinderknecht et al. 2014). The kinetic effects arise when length scales in the ICF implosion are comparable to, or smaller than, the mean free path of the ions or when ion velocities deviate from a Maxwellian distribution. We simulated the shock-wave propagation with MD and a multispecies, multicomponent continuum hydrodynamics model (Figure 3). The continuum model served as a standard for hydrodynamics. It includes diffusion, viscosity, and thermal conductivity, but not kinetic effects.

Figure 3.
Figure 3. Profiles of shock fronts in a H-T plasma mixture simulated with (a) continuum and (b) molecular dynamics (MD) techniques with the same initial conditions and shock strength (Mach 4.91), showing density (ρ), mole fractions (xH and xT), temperature (T), and pressure (p), normalized to their maxima. The shock wave propagates in the direction of the arrow. The length scales are significantly different for continuum and MD simulations, owing to kinetic effects (Mackay et al. 2017).

Many aspects of the shock propagation in the two simulations were found to be in quantitative agreement, including shock velocity and post-shock temperature and density. We observed separation of the two species at the shock front in both simulations, owing to their mass difference. The leading edge of the shock was found to be more abrupt in the continuum simulation, as shown in Figure 3. This difference might have been due to the viscosity, but we performed separate equilibrium MD simulations that verified the accuracy of the recently proposed kinetic molecular dynamics (KMD) viscosity model, even at the relatively cold pre-shock temperature of 1 eV. The KMD model was used in the continuum simulations. The MD simulations predicted a significantly larger thickness of the leading light-ion shock front, with some associated differences in the species separation. We also found anisotropic non-Maxwellian ion-velocity distributions in MD that were not enabled in the continuum simulations.

Mixing at the Carbon–Metal Interface in Warm Dense Matter Experiments

Experiments conducted at the Laboratory employed a laser-generated ion beam to rapidly heat a diamond–metal interface to warm-dense-matter conditions and used side-on radiography to measure the mixing at the interface (Bang et al. 2015, 2016). We used MOD-MD to simulate the mixing process in a carbon–silver system model also at solid density and rapidly heated to a few electron volts (Haxhimali et al. 2017). This system differed from the ablator–fuel interface (described above) in several important ways. The silver has a significantly larger atomic number than any of the components of the plastic-fuel system. During the heating, it was ionized to Z*~8, producing a high electron density and strong screening. The screening was sufficiently strong that the first-order screening polarization correction induced a layering instability (Haxhimali 2016). This issue was resolved and simulations of the diffusive mixing at the interface were conducted (Figure 4).

Figure 4.
Figure 4. Simulation of diffusive mixing at a C–Ag interface at warm-dense-matter conditions (solid density, T = 10 eV). Left: Simulation configuration showing initially sharp interfaces. Right: Number density as a function of position near the interface at three different times. The carbon has a higher number density at these constant pressure conditions. Also shown is the electric field as a function of position near the interface at three different times (Haxhimali et al. 2017).

Mixing at the C–Ag interface was found to occur by Fickian diffusion enhanced by the electrical field. The number density profiles and the electrical fields from a 10-eV simulation are shown in Figure 4. The electrical field was strong (approximately 1 GV/m), but confined to a small region in space, decaying with time. There was no indication of super-diffusive or ballistic transport that lasted more than a few ion-collision times, much shorter than the earliest (33 ps) curve plotted here.

Viscosity of Plasma Mixtures

We used MD to calculate the shear viscosity for asymmetric mixed plasma for thermodynamic conditions relevant to astrophysical and inertial confinement fusion plasmas (Haxhimali et al. 2015). Specifically, we considered D–Ar mixtures at temperatures of 100–500 eV and a number density of 1,025 ions/cc. The ions were taken to interact via the Yukawa (screened Coulomb) potential, for which the screening length was uniform (i.e., not varying with position). The electrical field of the electrons was included in this effective interaction. Shear viscosity was calculated using the Green-Kubo approach with an integral of the shear-stress autocorrelation function, a quantity calculated in the equilibrium MD simulations. We studied the effect of varying the fraction of the minority high-Z element (Ar) in the D–Ar plasma mixture. In more weakly coupled plasmas, at 500 eV and low Ar fractions, results from MD compared very well with the Chapman–Enskog kinetic results. In more strongly coupled plasmas, the kinetic theory did not agree well with the MD results, but a hybrid KMD viscosity model did, as shown in Figure 5. This hybrid interpolates between classical kinetic theories at weak coupling and the Murillo Yukawa viscosity model at higher coupling, a model based on single-species MD viscosity. We presented these results at the Charged Particle Transport Workshop at Sandia National Laboratories in October 2016 (Rudd et al. 2016).

Figure 5.
Figure 5. Plot (left) of the viscosity of a deuterium–argon plasma mixture with ion density 1025 ion/cc at three temperatures as a function of the mole fraction of argon. The points with error bars are MD results. The curves denote the kinetics molecular dynamics (KMD) viscosity model, a hybrid of viscosity values from kinetic theory and single-species MD, which is in good agreement (within 16 percent) with the MD viscosity for mixtures. The stress autocorrelation function as a function of time (right) is used to calculate viscosity. The same results are plotted in a log scale (inset) (Haxhimali et al. 2015).

Accelerating MOD-MD Simulations with the Fast-Multipole Method

We developed fast multipole algorithms and code designed for plasma-transport simulations. The initial implementation of MOD-MD computed all of the forces on an ion directly from the neighboring ions within a specified cutoff distance. As the temperature increases and screening becomes weaker, the cutoff must be increased, which requires more calculations at a greater computational cost. To mitigate this cost, there are several efficient methods for screened potentials, such as the Ewald summation, particle–particle-particle–mesh (PPPM), and the fast multipole method (FMM). So far, none of these methods can deal with a spatially varying screening length, as in MOD-MD. In analyzing these methods, only the FMM appeared to be a viable candidate for varying screening. We developed an FMM algorithm suitable for MOD-MD calculations that only increases the cost by O(log(N)), compared to traditional FMM, where N is the number of ions. We examined its implementation for non-cubic simulation boxes with periodic boundary conditions. Further developments are needed for full parallel FMM-MOD-MD simulations on supercomputers.

Impact on Mission

Plasma-transport processes are relevant to the Laboratory's mission for inertial confinement fusion and stockpile stewardship. We developed a new capability using MD techniques for transport calculations for asymmetric plasma mixtures, even in the presence of strong gradients such as those found at sharp interfaces and shock waves. Because this capability is based on molecular dynamics, it is valid even at moderately to strongly coupled conditions. The capability was delivered to the program to conduct plasma simulations with strong gradients. We demonstrated the performance of the code in massively parallel simulations on the Livermore Computing supercomputers, making use of Computing Grand Challenge allocations. This demonstration is relevant to the question of what computational capabilities are feasible for next-generation computing and influences the Laboratory's mission for predictive computation.


We developed a new capability to calculate transport properties of asymmetric plasma mixtures, i.e., plasmas with two or more components differing significantly in charge and mass, in the presence of strong gradients. We built this capability on the MD plasma-simulation code ddcMD developed for the Cimarron project to deliver a high-fidelity description of plasmas. We developed the multiscale orbital-free DFT molecular dynamics (MOD-MD) technique to enable simulation of ion dynamics in nonuniform plasmas (such as interfaces) without explicitly calculating the electron dynamics. MOD-MD was implemented in the ddcMD code. We used this capability to simulate diffusive mixing at a fuel–ablator interface in ICF, super-diffusive transport at shock fronts, and mixing at a diamond–metal interface under conditions similar to recent warm-dense-matter experiments. This work will continue through the Center of Excellence effort for next-generation code development (a Tier-1 Computing Grand Challenge project) and collaborations through the Laboratory's High Energy Density Science Center. The goal of this work is to understand other aspects of diffusive mixing processes at plasma interfaces using extremely large-scale MD simulations. We also submitted a proposal to use this new capability for various applications. The proposal is currently under review.


Amendt, P., et al. 2010. "Plasma Barodiffusion in Inertial-Confinement-Fusion Implosions: Application to Observed Yield Anomalies in Thermonuclear Fuel Mixtures." Physical Review Letters 105 (11):4. doi: 10.1103/PhysRevLett.105.115005.

Bang, W., et al. 2015. "Uniform Heating of Materials into the Warm Dense Matter Regime with Laser-Driven Quasimonoenergetic Ion Beams." Physical Review E 92 (6):10. doi: 10.1103/PhysRevE.92.063101.

——— 2016. "Linear Dependence of Surface Expansion Speed on Initial Plasma Temperature in Warm Dense Matter." Scientific Reports 6:8. doi: 10.1038/srep29441.

Bellei, C., et al. 2013. "Species Separation in Inertial Confinement Fusion Fuels." Physics of Plasmas 20 (1):8. doi: 10.1063/1.4773291.

Chapman, S. and T. Cowling. 1939. "The Mathematical Theory of Non-Uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction, and Diffusion in Gases." London, Cambridge University Press.

Frenkel, D. and B. Smit. 2001. "Understanding Molecular Simulation: From Algorithms to Applications, Computational Science." New York, Elsevier Science.

Glosli, J. N., et al. 2010. "Extending Stability Beyond CPU Millennium: A Micron-Scale Atomistic Simulation of Kelvin-Helmholtz Instability." Proceedings of the ACM/IEEE Conference on High Performance Networking and Computing, Reno, Nevada, USA, 10–16 November 2007.

Graziani, F. R., et al. 2012. "Large-Scale Molecular Dynamics Simulations of Dense Plasmas: The Cimarron Project." High Energy Density Physics 8 (1):105–131. doi: 10.1016/j.hedp.2011.06.010.

Haxhimali, T. 2016. "Atomistic Study of Mixing at High Z/ Low Z Interfaces at Warm Dense Matter Conditions." 59th Annual Meeting of the APS Division of Plasma Physics Meeting, San Jose, CA.

Haxhimali, T., et al. 2015. "Shear Viscosity for Dense Plasmas by Equilibrium Molecular Dynamics in Asymmetric Yukawa Ionic Mixtures." Physical Review E 92 (5):14. doi: 10.1103/PhysRevE.92.053110.

Mackay, K. K., et al. (Forthcoming.) "Species Separation and Other Kinetic Effects at Shock Wave Fronts in Plasma Mixtures Assessed Using Molecular Dynamics." Journal of Chemical Physics.

Rinderknecht, H. G., et al. 2014. "First Observations of Nonhydrodynamic Mix at the Fuel-Shell Interface in Shock-Driven Inertial Confinement Implosions." Physical Review Letters 112 (13):5. doi: 10.1103/PhysRevLett.112.135001.

Rudd, R. E., et al. 2016. "Plasma Viscosity Calculated with KMD." Charged Particle Workshop, Sandia National Laboratories, Albuquerque, NM.

Sanbonmatsu, K. Y. and M. S. Murillo. 2001. "Shear Viscosity of Strongly Coupled Yukawa Systems on Finite Length Scales." Physical Review Letters 86 (7):1215—1218. doi: 10.1103/PhysRevLett.86.1215.

Spitzer, L. 1962. "Physics of Fully Ionized Gases." Geneva, Interscience Publishers.

Stanton, L. G., et al. (Forthcoming.) "A Multiscale Molecular Dynamics Model for Heterogeneous Charged Systems." Physical Review.

Whitley, H. D., et al. 2015. "Molecular Dynamics Simulations of Warm Dense Carbon." Contributions to Plasma Physics 55 (5):390—398. doi: 10.1002/ctpp.201400101.

Publications and Presentations

Glosli, J. N. 2015. “The Role of High-Performance Computing in Coulomb Systems.” AFOSR Workshop, East Lansing, MI. LLNL-PRES-691210.

Haxhimali, T. 2015. “Shear Viscosity for Dense Plasmas by Equilibrium Molecular Dynamics in Asymmetric Yukawa Ionic Mixtures.” 58th Annual Meeting of the American Physical Society Division of Plasma Physics, 17 November 2015. LLNL-PRES-679234.

Haxhimali, T., et al. 2015. “Shear Viscosity for Dense Plasmas by Equilibrium Molecular Dynamics in Asymmetric Yukawa Ionic Mixtures.” Physical Review E 92 (5):053110. LLNL-JRNL-674655.

Haxhimali, T. 2016. “Atomistic Study of Mixing at High Z/ Low Z interfaces at Warm Dense Matter Conditions.” Annual Meeting of the American Physical Society Division of Plasma Physics, 3 November 2016. LLNL-PRES-707978.

——— 2017. “An Atomistic Study of Transport and Mixing in High Energy Density Plasmas.” National Tsing Hua University, Hsinchu City, Taiwan, 6 September 2017. LLNL-PRES-737900.

Mackay, K. K., et al. (Forthcoming.) “Species Separation and Other Kinetic Effects at Shock Wave Fronts in Plasma Mixtures Assessed Using Molecular Dynamics.” Journal of Chemical Physics. LLNL-JRNL-733037.

Nicholas, M., et al. 2017. “Atomistic Study of Mixture at the Silver-Carbon Plasma Interface at Uniform High Temperatures.” North Carolina Central University, August 2017. LLNL-POST-736563.

Noble, B. A., et al., 2016. “Shear Viscosity of Asymmetric Strongly Coupled Dense Plasmas.” Lawrence Livermore National Laboratory Summer Student Symposium, August 2016. LLNL-POST-699161.

Rudd, R. E. 2016. “Species Diffusion in Plasma Mixtures.” Annual Meeting of the American Physical Society Division of Plasma Physics, San Jose, CA, 1 November 2016. LLNL-PRES-707402.

Rudd, R. E. and T. Haxhimali. 2016. “Plasma Viscosity Calculated with KMD.” Charged Particle Workshop, Sandia National Laboratories, Albuquerque, NM, October 25, 2016. LLNL-PRES-706641.

Stanton, L. G. 2014. “A Multi-Scale Molecular Dynamics Model for Heterogeneous Charged Systems.” IPAM Workshop, Lake Arrowhead, CA, 10 December 2014. LLNL-PRES-663470.

Stanton, L. G., et al. (Forthcoming.) “A Multiscale Molecular Dynamics Model for Heterogeneous Charged Systems.” Physical Review. LLNL-JRNL-732633.