Multiscale Molecular Dynamics Model for Heterogeneous Charged Systems

Modeling matter across large length scales and timescales using molecular dynamics simulations poses significant challenges. These challenges are typically addressed through the use of precomputed pair potentials that depend on thermodynamic properties like temperature and density; however, many scenarios of interest involve spatiotemporal variations in these properties, and such variations can violate assumptions made in constructing these potentials, thus precluding their use. In particular, when a system is strongly heterogeneous, most of the usual simplifying assumptions (e.g., spherical potentials) do not apply. Here, we present a multiscale approach to orbital-free density functional theory molecular dynamics (OFDFT-MD) simulations that bridges atomic, interionic, and continuum length scales to allow for variations in hydrodynamic quantities in a consistent way. Our multiscale approach enables simulations on the order of micron length scales and 10’s of picosecond timescales, which exceeds current OFDFT-MD simulations by many orders of magnitude. This new capability is then used to study the heterogeneous, nonequilibrium dynamics of a heated interface characteristic of an inertial-confinement-fusion capsule containing a plastic ablator near a fuel layer composed of deuterium-tritium ice. At these scales, fundamental assumptions of continuum models are explored; features such as the separation of the momentum fields among the species and strong hydrogen jetting from the plastic into the fuel region are observed, which had previously not been seen in hydrodynamic simulations.


I. INTRODUCTION
For a wide range of problems in science and technology, essential descriptive information can be found over multiple timescales and/or length scales [1,2].In contrast, computational methods typically take advantage of limited ranges of scales; for example, molecular dynamics (MD) is used to study systems at the microscopic scale, and hydrodynamics, at the macroscopic scale.Information from other scales can be included in such methods in a variety of ways.Typical hydrodynamic codes incorporate subgrid microscopic information through precomputed, near-equilibrium equations of state and transport coefficients.When this is not possible, frameworks such as the heterogeneous multiscale method [3] couple scales by linking disparate models, usually with MD providing closure information to an incomplete macroscale model.Conversely, a brute-force strategy is to rely on the development of ever-faster algorithms and hardware that allow MD to natively reach relevant scales [4,5].However, obtaining MD forces through an on-the-fly electronic structure calculation remains relatively expensive [6,7], considerably limiting achievable timescales and length scales.
High-energy-density (HED) physics is an important example of a field of research for which computational tools are essential: HED experiments are both very expensive and relatively poorly diagnosed, so computational methods play a crucial role in our understanding of HED environments.For example, facilities such as the Linac Coherent Light Source [8], the National Ignition Facility [9,10], and the Omega Laser Facility [11] regularly produce highly transient and heterogeneous HED matter into which modeling codes provide critical insights.Traditionally, such modeling codes use hydrodynamic models [12] coupled to precomputed equations of state [13,14].
At the same time, the past decade has seen enormous progress in the development of MD methods for the study of HED physics.MD codes based on Kohn-Sham-Mermin density functional theory (DFT) can reach scales of about 10 Å, at the price of severely limiting the method to uniform properties of matter, and the method also becomes more expensive at higher temperatures.Thus, orbital-free DFT (OFDFT) methods have been developed [15][16][17][18][19][20][21][22][23] to overcome some of these limitations.For zero-temperature systems with fixed nuclei, finite element methods have been successful in scaling OFDFT calculations to very large systems [24]; however, coupling this approach to the MD equations of motion would require an adaptive mesh generation routine to be run at each time step for each configuration of nuclear coordinates, which could introduce a significant bottleneck to the simulation.OFDFT methods are particularly well suited to HED matter because higher temperatures generate enough disorder to decrease the importance of subtler quantum effects (e.g., band structure, bonding, etc.).Furthermore, many OFDFT models become increasingly accurate at the higher densities encountered in HED systems.However, even the simplest OFDFT models are expensive enough to prohibit reaching length scales and timescales relevant for important nonequilibrium processes in heterogeneous experiments.As a result, large-scale MD simulations of HED matter tend to employ simplified potentials [25].
Here, we are interested in new methods that permit MD simulations of nonequilibrium, heterogeneous HED physics environments using OFDFT-based techniques.We develop a multiscale model that computes interatomic forces on the fly, without assuming a precomputed pair potential (e.g., of the Yukawa form), while also including a fast electronicstructure calculation that allows us to explore hydrodynamic phenomena in moderately to strongly collisional systems.This paper is organized as follows.In Sec.II, we present equations of motion for a system of interacting ions and electrons.We then develop the multiscale model in its most general form in Sec.III, with some details given in the appendixes.A specific application of the model is presented in Sec.IV; atomic mixing at an interface typical of those found in the ablator-fuel region of an inertial-confinementfusion (ICF) capsule is used as the system of study.Finally, general conclusions and a discussion of the results are given in Sec.V. Throughout this work, the equations will be presented in Gaussian-cgs units (i.e., the Coulomb constant is unity) unless otherwise stated, with temperature being expressed in units of energy, where the Boltzmann constant k B has been absorbed into the temperature.

II. EQUATIONS OF MOTION
At the core of our model is a collection of nuclei, with coordinates fr i ðtÞg and velocities fv i ðtÞg, surrounded by a distribution of electrons n e ðr; tÞ, the details of which are discussed in the next section.The classical equation of motion for the ith nucleus of mass m i is given by where the force F i has been decomposed, without loss of generality, into its nuclear and electronic components, respectively.In general, these components will each still depend on both the nuclear and electronic degrees of freedom.Because of the large difference in mass between nuclei and electrons, the Born-Oppenheimer approximation, in which the electrons are assumed to evolve implicitly with the nuclei, is often invoked to drastically reduce the computational cost of simulations.We have partially relaxed this approximation [26] by allowing the forces resulting from the electrons to be decomposed into "fast" and "slow" components as The slow component can be calculated in the limit of m e =m i → 0, where m e and m i are the electron and ion masses, respectively.In this approximation, the electron density responds instantaneously to the motion of the nuclei; the calculation of this term will be discussed in detail in Sec.III.The fluctuations about this density are then captured in the remaining fast component, and we treat this term using a simple Langevin model: As with most Langevin models, this form is purely dissipative [27], meaning that it does not add or subtract any systematic forces.This can be seen by noting that if the ions were in a frozen configuration and the average force were computed, hF fast i i would be zero.It is in this sense that we place the electrostatic forces associated with the electronic structure into the slow portion of the force and assume that the two contributions are "orthogonal." The Langevin parameter γ ie is chosen subject to the constraint that a nucleus will slow down because of electron drag, in accordance with an appropriate theoretical model [28,29].The noise term in Eq. (3) has zero mean [hξðtÞi ¼ 0], and its correlation properties are determined by the fluctuation-dissipation theorem to ensure that the nuclei tend toward the electron-bath temperature T e , resulting in hξ i ðtÞξ j ðt 0 Þi ¼ ð2γ ei T e =m e Þδ ij δðt − t 0 Þ, where δ ij is the Kronecker delta and δðtÞ is the Dirac delta function.Dai and co-workers [30,31] have discussed Langevin models of this form; however, we use a different approach to determine the Langevin parameters.
To study the fast portion of the force and obtain the parameter γ ei , consider a single nucleus traveling through the electron gas with no nucleus-nucleus interactions.For an ensemble of such nuclei (averaged over the noise), we can find the mean velocity hv i i and the mean kinetic energy E i ¼ ðm i =2Þhv i i 2 .Using Eq. (3) in Eq. ( 1) and setting the slow component to zero, the rate of change of the kinetic energy can be found from the average of the solution of Eq. ( 1), to yield a stopping power (energy loss per unit distance) of the form thereby connecting the frequency γ ie to the low-velocity stopping power [32,33].While we could similarly apply a low-velocity stopping power formulation, such as the model of Nagy et al. [34], we desire a model that both applies to a wide range of plasma conditions and can be computed extremely rapidly.For this reason, we have employed the model of Skupsky [35].Comparing the Langevin prediction (4) with Skupsky's model [35], we find where η ¼ μ=T e , μ is the chemical potential, Z i is the charge of the ith nucleus, e is the elementary charge, ℏ is the reduced Planck constant, and the effective Coulomb logarithm is given by [36] ln The Coulomb logarithm ( 6) is naturally convergent as a result of long-wave (k → 0) screening characterized by the Thomas-Fermi (TF) length scale k −1 TF and of a purely quantum short-wave (k → ∞) cutoff on the scale of the thermal de Broglie wavelength λ T ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2πℏ 2 =ðm e T e Þ p .The frequency γ ie is valid for all degeneracies, and the values of η and ln Λ eff can be computed using the relations where the pth-order Fermi integral is defined here as We have also introduced the degeneracy parameter, which is defined in terms of the electron temperature and Fermi energy as . Accurate and efficient fits to these integrals and their inverses can be found in Refs.[37][38][39].To avoid the integration in Eq. ( 6), a simple interpolation formula is also presented in Ref. [35], which approximates the effective Coulomb logarithm as

: ð10Þ
Note that the variables n e ðr; tÞ, ηðr; tÞ, and T e ðtÞ vary throughout the simulation discussed in Sec.IV, so the nuclei experience energy transfer to and from the electrons consistent with local conditions.In other words, we convert all local, uniform models, such as the stopping-power model and the internal-energy model given below, to their local-density approximations; this is a potential source of error but one that is difficult to quantify because of a lack of stopping models that incorporate macroscopic gradients in input quantities such as the dielectric response function [40].
The Langevin forces on the nuclei are associated with momentum and energy transfer between the electronic and nuclear subsystems, subject to energy and momentum conservation [41][42][43].For example, fast nuclei will experience electronic stopping, which in turn causes electron heating.The electron-bath temperature therefore evolves as where Ω is the volume of the simulation cell, In the second line, the sum over each nucleus i is replaced with a sum over the species j.The electronic heat capacity C e ¼ ∂u e =∂T e can be computed using an equation of state given by an ideal Fermi gas with specific internal energy Finally, the term S e ðtÞ is the power delivered to the electrons by an external source (e.g., particles and/or radiation).In Sec.IV, we employ this model to examine the heating of an interface; however, we use a reduced form of the model in which we inject energy into the electrons only through a prescribed T e ðtÞ, and there is no backreaction from the nuclei.

III. MULTISCALE OFDFT MD METHOD
The principal challenge in evolving Eq. ( 1) is the evaluation of the forces F i on each nucleus.In particular, calculations of the screened forces arising from the combination of the nuclei and the "slow" component of the electrons, which we denote as , are the most computationally demanding because of the manybody and potentially long-range nature of these forces.In the absence of magnetic fields, these forces result from the electrostatic interactions between the charge distributions and can thus be expressed in terms of an electrostatic potential and the corresponding Poisson equation as The electronic density n e ðrÞ then remains to be determined.The choice of an electronic structure is needed before proceeding; we compromise between the computational cost and physics fidelity.For computational reasons, it is desirable to exploit the lower computational costs of OFDFTs.Within that class of methods, we must choose the forms of the kinetic-energy and exchange-correlation functionals.First, our model must accurately predict the ionization level for a high-fidelity separation of the microscale from the macroscale.Both Fromy et al. [44] and Murillo et al. [45] have shown that the ionization level predicted from various models is insensitive to the model choice above approximately 10 eV.Similarly, Fromy et al. [44] show that the pressure is also insensitive, which is important for obtaining accurate effective ion-ion forces.Yonei, Ozaki, and Tomishima [46] have found similar results for the exchange-correlation potential, showing that the exchange and correlation free energies decrease rapidly at higher temperatures and that the most important correction to the basic TF model is the gradient correction, consistent with the observations of Fromy et al. [44].In comparison with Kohn-Sham DFT, Danel, Kazangjian, and Zerah [47,48] have shown that the Thomas-Fermi-Weizsacker method converges to the Kohn-Sham result rapidly above approximately 10 g=cc for boron.For hot, dense plasmas, the functionally similar Thomas-Fermi-Kirzhnits (TFK) model has been found to be accurate [49,50], as the gradient expansion is derived in the longwavelength limit.Thus, from a physics-fidelity point of view, we begin with the TFK formulation, which is among the least computationally expensive approaches; more details are given below.While this choice may not be justifiable in other applications, such as those at lower temperatures [51], our numerical study will elucidate the relative costs of each portion of the multiscale method.
We now turn to a complete description of our OFDFT model, which is formulated from the variational minimization of an approximate free-energy functional.Consider the grand potential [52] for the electrons in the presence of an external potential produced by the nuclear cores: where F ½n e is the Helmholtz free energy, and μ is the chemical potential, which ensures conservation of particle number.Within the total free energy (17), F e ½n e is the free energy of noninteracting electrons, the second term is the Hartree electron-electron interaction, v ext ðrÞ is the external potential arising from the nuclear cores, and F xc ½n e is the remaining exchange-correlation contribution.The electronic structure is then obtained through the variational minimization δΩ δn e ½n e ¼ 0; ð18Þ which yields the Euler-Lagrange (EL) equation Once a particular form of the free energy ( 17) is prescribed, OFDFT can be used to determine the electronic density at a particular time step, together with standard MD (i.e., without the Langevin terms) to evolve the nuclear coordinates; this combined method is called the "OFDFT molecular dynamics" (OFDFT-MD, or sometimes just OF-MD) method [15,[53][54][55].The OFDFT-MD formulation is preferred here over the more commonly employed orbital-based approaches because of the computational advantage it offers with little loss of accuracy at high temperature; however, aspects of our model can be modified to include orbitals (e.g., in the average-atom calculation) for lower-temperature applications.Note that this notation is consistent with the commonly used convention of describing the electronic method before the dash and the ionic method after the dash.
To illustrate our multiscale method, we choose a specific macroscale model relevant to HED matter; however, the method is generic to an arbitrary choice in the functional.For HED environments, the random-phase approximation, in which exchange and correlation effects are ignored (F xc ≈ 0), can often be invoked, as these contributions become increasingly negligible in the high-temperature limit [46].A reasonable starting point to model the remaining free energy of the noninteracting electrons, F e , is the TF functional This functional is exact in the (nonrelativistic) dense and/or hot limit, where the electron distribution becomes increasingly uniform.Equation (19) thus reduces to Alternatively, Eq. ( 22) can be written in differential form as the Poisson-Thomas-Fermi equation where β ¼ 1=T e is the inverse thermal energy of the electrons.An advantage of the DFT formalism is that one can systematically improve upon the model by adding new contributions to Eq. ( 16).For example, as the TF functional assumes a uniform electron density, we can explore the effects of gradient corrections to the free energy, where the finite-temperature Kirzhnits correction [49] to the kinetic energy functional is given by The prime in Eq. ( 27) denotes differentiation with respect to η, and the parameter λ allows the model to span both the true gradient-corrected TF limit (λ ¼ 1=9) and the traditional Weizsäcker correction at T e ¼ 0 (λ ¼ 1).The function h λ ðηÞ, in turn, ranges from ℏ 2 λ=ð8m e Þ in the strongly degenerate limit (Θ ≪ 1) to 3ℏ 2 λ=ð8m e Þ in the weakly degenerate limit (Θ ≫ 1).An accurate fit to h λ ðηÞ for the λ ¼ 1=9 case (denoted simply as h) can be found in Ref. [49], and further gradient corrections can be found in Ref. [17].In differential form, the EL equation corresponding to the TFK model becomes Both the TF and TFK models offer multiple advantages; they become exact in the dense limit, are valid over all finite temperatures, and yield simpler calculations than their orbital-based counterparts.Furthermore, they are sufficiently accurate over large regions of the HED parameter space [50].
It is important to note that even with the simpler TF model, the evolution of the OFDFT-MD equations can still be computationally challenging, as the calculations of the electron density and resulting forces are typically nonlinear and nonlocal and involve a variety of length scales.Current attempts at a direct solution of the OFDFT-MD system, which require empirical modifications to the equations to circumvent numerical instabilities arising from singular electron densities near nuclear cores [56], are still limited to thousands of nuclei and hundreds of femtoseconds [57].We attempt to solve these same equations but mitigate some complexities of the calculations by separating the model into components occurring on three length scales: intraatomic, interatomic, and continuum.At the intra-atomic (micro) scale, atomic physics, in which complex interactions occur between bound (core) electronic states and free (valence) states, dominates.At the interatomic (meso) scale, the dynamics is governed by ionic collisions and many-body phenomena.Finally, at the continuum (macro) scale, large-scale variations in the bulk parameters can naturally be described.The remainder of this section is devoted to calculations at these three scales and ultimately to recombining the calculations.For brevity, we refer to this Multiscale OF-DFT MD approach as "MOD-MD" within the remainder of this paper.

A. Microscale: Ionization states
At the smallest length scale in our model, we can greatly simplify our calculations by treating bound (core) and free (valence) electrons of a given atom separately.Once determined, the free states of the electrons can be distributed throughout the system, and the remaining bound states are treated as point particles suppressed to the nucleus, giving the resulting ion an effective charge of Z Ã i e.As this ionization state might not be known a priori and could depend on the local density and temperature for the particle's position, we present a model to approximate Z Ã i e for each particle.It should be noted that we use the notation Z Ã to refer to the mean ionization state (MIS) in general, rather than a specific theoretical definition [45].
We begin by introducing a fixed, coarse-grained mesh with grid spacings that are large compared to the mean free paths of the ions, as illustrated in Fig. 1.It is then assumed that the electrons in each cell between mesh points are in local thermodynamic equilibrium, and thus the local temperature and mean ionic density within the cell are used to calculate electronic properties.A multispecies ionization model that can be used to rapidly compute the finitetemperature and finite-density atomic physics is required.If a precomputed table of values is not readily available, then several options for computing the ionization states on the fly are available.Saha models, which are accurate in the lowdensity limit, are derived from a chemical picture in which experimentally measured ionization energies are combined with Boltzmann transition probabilities and known degeneracies for each state [58,59].Density effects can be included in a limited sense through accurate continuumlowering models.Alternatively, average-atom (AA) models, which are accurate in the high-density limit, can be used to calculate the spherically symmetric electronic structure about a nucleus of charge Ze within a cavity embedded in a jellium background.The simplified geometries used in AA models allow the rapid calculation of various thermodynamic quantities, such as the ionization state.
For our own calculations, we use the well-studied TF average-atom (TF-AA) model [60], for which accurate fits are available [61]; this model has also often been shown to agree well with more sophisticated models, and it becomes increasingly accurate for hot, dense systems [45].In this model, the electrons are approximated by minimizing the TF functional (20) in a spherical volume chosen to have the ion-sphere radius a i ¼ ð4πn ion =3Þ −1=3 , where n ion is the mean ionic number density of the coarse-grained cell.On the sphere, the electric field is assumed to be zero, and the electrons are assumed to be free beyond the sphere.The approximate ionization state is thus calculated as The TF-AA model can easily be generalized to mixtures of multiple atomic species, as well, by iteratively converging on a self-consistent chemical potential within the numerical cell; the numerical scheme for this generalization is described in detail in Appendix A. Decomposing the density as a superposition of average atoms at the microscale step has the distinct advantage that the MOD-MD steps begin with a finite-temperature, finite-density, allelectron calculation that does not employ a pseudopotential nor require careful treatment of a grid near the nuclei.Alternatively, because we ultimately obtain a set of fZ Ã i g for the species from this step, the average-atom step can also be viewed as including the pseudopotential calculation in line, where it is computed on the fly as local conditions change in space and time.Once the ionization state is calculated for a given particle, the free and bound states are separated, and the bound states are taken to be point charges located at the nucleus.This allows ionization states to be assigned to ions and the total number of free electrons to be determined, as well.More importantly, the macroscale calculation need only be applied to the free electrons, as we discuss next in Sec.III B. Note that this will also alter the damping rate (5) in the Langevin process, so it will now be expressed in terms of the ionic charge number as

B. Macroscale: Coarse-grained fields
Once an ionic charge is assigned to each nucleus in the above calculation, a coarse-grained (CG) ionic charge density ρ cg ion ðrÞ can be constructed on the grid points associated with the CG mesh.Using this quantity, we can then calculate the corresponding CG free-electron density n cg e ðrÞ and electrostatic potential Φ cg ðrÞ selfconsistently using the corresponding OFDFT model [Eqs.( 16)- (19)] at the CG level.As with the full model, we use the grand potential for the CG free electrons in the presence of an external potential produced by the CG ionic charge density: Using this potential, we obtain the following EL and Poisson equations to be solved on the CG mesh: There are two advantages to this approach.First, the coarse-grained electronic density is smoother than the total density, greatly relaxing the conditions on the grid spacing relative to what would be needed for the bound electrons near the nuclei.Second, the macroscale density is allowed to have a different symmetry than the total density.For example, a planar interface can be approximated as one dimensional, as in the application discussed in Sec.IV.In that simulation, we construct the CG mesh by dividing the numerical domain at that scale into bins parallel to the interface, with enough resolution to capture the gradient of the free-electron density near the interface.Capitalizing on both the smoothness and the higher symmetry in the FIG. 1. Coarse-grained mesh.Illustration of point ions and coarse-grained mesh used to find the local ionization levels of the mixtures in each cell.Note that, in practice, the mesh and ion coordinates are in three dimensions, and the mesh spacings need not be equal in each direction.macroscale electronic structure calculation considerably reduces the computational cost.We would like to emphasize several points that highlight the multiscale aspect of the MOD-MD model.The grid is chosen to match the macroscale variations in the free-electron density, which allows the grid to be in fewer than three dimensions and to have a spacing that resolves only free-electron variations.For example, for problems such as planar interfaces and simple shocks, this allows the use of a 1D grid that is fine in a subregion near the largest inhomogeneities and very coarse in other regions.This approach stands in stark contrast to other DFT methods that enforce full self-consistency by using a grid that resolves all scales associated with all of the electrons.
Within each spatial zone, we begin each time step by computing the atomic-physics properties for the new composition, temperature, and number densities for that time step.These calculations are necessary when studying HED matter because the ionization of the time-evolving mixtures is not known a priori before the simulation.

C. Mesoscale: Effective pair interactions
Decomposing the electron density into free states distributed throughout the system and bound states suppressed to their respective nuclei, we can approximate Eq. ( 15) as where n f ðrÞ denotes the free electrons.By introducing the free-electron-density fluctuation Δn e ðrÞ ¼ n f ðrÞ − n ðeÞ 0 and electrostatic potential fluctuation ΔΦðrÞ ¼ ΦðrÞ − Φ 0 , where n ðeÞ 0 is the mean free-electron density, which acts as a uniform neutralizing background, and Φ 0 is the mean electrostatic potential, we can then rewrite the above expression in the form where Here, KðrÞ can be thought of as a generalized screening function; note that this function becomes constant to the value of the inverse screening length in the Debye-Hückel-Yukawa limit.We exploit this fact by assuming that KðrÞ is itself a slowly varying function in space, even if its constituent functions (Δn e and ΔΦ) are not; this approximation is the central assumption of this multiscale model.In this slowly varying limit, we thus approximate KðrÞ using the CG fields as Δn cg e ðrÞ ΔΦ cg ðrÞ : Once the function KðrÞ is known from the CG calculation, Eq. ( 37) is used to calculate the resulting electrostatic potential, and the particle positions can be numerically updated from Eq. (1); the entire process is then repeated.
Figure 2 shows a schematic of the concurrent steps used during each time step of MOD-MD to evolve the particle positions, and we show how the various length scales enter into solving Eq. (37).Deeper insight into approximation (38) is gained by examining the limits of the CG mesh size.In the extreme limit of an infinitely resolved 3D mesh, the CG ionic charge distribution simply recovers the bare charge distribution of discrete point nuclei, and thus the calculation of Eq. ( 37) recovers the full all-electron equations of OFDFT.Conversely, the opposite limit of an infinitely coarse mesh yields a constant function K calculated from the mean free electron density, which is identical to a system described by screened Coulomb (or Yukawa) interactions.Using an example system of mixing at an interface, which will be explored in great detail in Sec.IV, we can demonstrate the convergence of K on the CG mesh.In this example, the interface is planar, so a 1D mesh can be used to approximate the full system, where the direction normal to the interface is denoted by the z axis (note that this 1D approximation will break down as the mesh becomes infinitely resolved).Calculations of the screening function KðzÞ for various bin sizes are shown in Fig. 3 an interpolation scheme is used to make KðzÞ continuous throughout the system.As the number of bins within the fixed domain increases, the solution converges fairly rapidly.This convergence can be quantified using the relative error where K c is the converged function, which is plotted in the bottom row of Fig. 3.As expected, the smoother case at the later time converges more rapidly (by roughly an order of magnitude).
For the general problem, calculating KðrÞ separately at the macroscale simplifies Eq. (37) to a (quasi)linear equation at a particular time step, and this simplification allows us to decompose the many-body problem into a superposition of one-body problems as As solving Eq. (41) for each particle could still be computationally expensive, we also exploit the slow variations in the function KðrÞ and derive perturbative solutions in this limit; see Appendix B for details of these calculations.At the lowest-order approximation, where variations in KðrÞ can be considered locally negligible, the Oð1Þ asymptotic solution to the one-body problem in Eq. ( 41) is given by Note that this solution is equivalent to the standard Yukawa model with a locally defined screening length at each particle's position.We can easily improve upon Eq. ( 42) by including the next-order corrections to obtain We can now calculate the screened force (14) on the ith particle (i.e., the combination of the nuclear and slowelectronic contributions), modified by the macroscopic variations in the system.Using Eq. ( 42), the Oð1Þ force is asymptotic to where r ij ¼ r i − r j and r ij ¼ jr ij j.In Eq. ( 44), we find that the so-called "self-energy" terms vanish; however, at the next order, the anisotropy of the screening cloud due to the heterogeneous background induces a dipole contribution to the force.To understand the validity of expanding in terms of slow variations in the screening function KðrÞ, which is calculated using a CG mesh, we again employ the example of a planar interface used to calculate Eq. (39).For this expansion to be valid, the magnitude of the expansion terms must be successively smaller, and thus the first-order correction ðr=2Þr • ∇Kðr i Þ from Eq. ( 43) needs to be small compared to unity.On the 1D CG mesh under consideration, this correction will be largest along the z direction and must be small for a given bin width.Evaluating these correction magnitudes over the entire domain, we thus require where L is the bin size.As expected, these maximal corrections occur near the large gradients induced by the interface, and we show C max as a function of the number of bins N in Fig. 4. The two cases of a sharp and a diffuse interface are again examined, where it can be seen that over a thousand bins are needed to keep these corrections small for the initial interface, while only several hundred bins are needed for later times in the simulation.This correction appears to roughly scale as N −2 for large N, which is expected to break down once individual particles are being resolved and the 1D approximation is no longer valid.However, at that point, the CG calculation will begin to recover the full OFDFT solution, making the expansion obsolete.Using Eq. ( 43) instead, the next-order correction to the force yields the improved form Higher-order corrections to both the electrostatic potential and the corresponding force are presented in Appendix C.
To understand the source of the dipole term in Eq. ( 46), we examine a screening cloud of electrons around a given ion; this screening can be calculated from At Oð1Þ, we have the traditional Yukawa form n ðiÞ e ðrÞ ¼ Z Ã i κ 2 i e −κ i r =ð4πrÞ, which is clearly spherically symmetric.The dipole appears at the next order, where we have While n ðiÞ e ðrÞ can be negative in this approximation, recall that this quantity is the deviation in the density with respect to the mean.In general, the nth-order asymptotic expansions in Eqs.(B5) and (B6) will produce the screening cloud where ε is the expansion parameter associated with gradients in K 2 ðrÞ; see Appendix B. We present plots of these corrections to the electron density in Fig. 5. Here, we have used κ i ¼ 1, h 1 ¼ 1=2, and h 2 ¼ 1=4 in arbitrary units, where h 1;2 are coefficients defined in the expansion (B6), and we have assumed an electric field directed solely in the negative-z direction.The top row of the figure shows the firstorder correction corresponding to Eq. ( 49), and the bottom row shows the next-order correction corresponding to Eq. (C17).Within each row, the left panel shows rn e ðrÞ as a color map in the (x, z) plane, while the right panel shows a cross section of the same quantity along the z axis, with the unperturbed density plotted (dashed line) for comparison.
The electric field can be seen to cause the electron cloud to shift to the right, which in turn weakens the force on the ion, and the next-order correction appears to enhance this polarization.
In summary, the MOD-MD method separates the system into two timescales (fast and slow) and three length scales (atomic, interionic, and continuum) as shown in Fig. 2. The microscale is used to resolve atomic physics, while the macroscale captures variations in the bulk parameters.These scales are then combined to inform the screened interactions between ions at the mesoscale.Finally, the nuclear positions are updated using these screened forces, along with the electronic fluctuations associated with the FIG. 4. Maximal correction.Maximum magnitude of the firstorder correction in Eq. ( 43) evaluated across a bin given by the expression (45) as a function of number of bins.The two cases shown are for an initial sharp interface at t ¼ 0 ps (red curve, circles) and a later diffuse interface at t ¼ 15 ps (green curve, triangles).The curve N −2 is also plotted to show the rough scaling of the convergence.faster scale in accordance with Eq. ( 1).The models that encompass MOD-MD are summarized in Table I.While we have argued that good modeling choices that balance physics fidelity and computational cost have been made, we also show a column with obvious improvements that could be made.Which improvements are needed depends, of course, on the specific problem, and the computational costs are currently mostly unknown.However, the important aspect of this multiscale approach is that its main approximation is breaking with perfect self-consistency, which allows improvements to various model inputs without necessarily impacting the others.Next, we apply this model to a specific application for which we can examine the relative computational costs.

IV. INTERFACIAL MIXING IN HIGH ENERGY-DENSITY MATTER
We now illustrate the use of the MOD-MD method for studies of mixing at interfaces.Several specific aspects of such systems can be studied using this approach, including microscale mixing rates computed at MD fidelity (i.e., across coupling regimes), the nature of transport (e.g., the role of electric fields), and nonhydrodynamic features (e.g., anisotropic velocity distributions).While MD has been used previously to examine electric fields near metalaqueous interfaces [62], our interest is in HED environments.In particular, targets in ICF experiments typically have a central deuterium-tritium (DT) fuel surrounded by higher-Z elements that form pusher-ablator layers; mixing of the higher-Z elements into the fuel can produce an undesirable energy loss (through enhanced bremsstrahlung emission).In some cases, the fuel is initially separated into an interior gas volume surrounded by a solid (ice) fuel layer, both composed of a DT mixture.The targets are compressed by shocks produced in the ablator layers, and these shocks can potentially cause mixing (or unmixing) at the interface.We examine the preshock conditions at the ice-ablator interface, which is preheated before the arrival of the shock resulting from radiation and fast electrons.Of particular interest are the transport processes [63] at the heated interface as a function of the heating rate and how these processes impact the material properties and interface structure that the shock will experience.
The simulations were set up as follows.Because transport processes scale as T 3=2 in the high-temperature limit, FIG. 5. Polarized screening clouds.Comparisons between the perturbed and unperturbed electron densities.On the left, we plot a color map of rn e ðrÞ in the (x, z) plane, where Eq. ( 49) is used for the top row and Eq.(C17) is shown in the bottom row.On the right, we have plotted these same quantities as cross sections along the z axis, with the perturbed density shown as a dashed curve for comparison.In each case, we have set and h 2 ¼ 1=4 in arbitrary units, where the coefficients h 1;2 are defined in Eq. (B6).TABLE I.A summary of the various inputs to the MOD-MD model is given, along with improvements we feel are the next obvious step.Because of the multiscale nature of MOD-MD, most models can be improved independently; however, improvements are also possible by enforcing higher degrees of self-consistency among the models.their dominant contributions will occur at the highest temperature of the simulation (50 eV).As such, we focus on the relevant physics at that temperature scale and thus use the basic Thomas-Fermi formulation of OFDFT for both the microscale and macroscale components of the model as detailed in Table I [44][45][46][47][48]. Periodic boundary conditions were used in all Cartesian directions, which results in two interfaces; in most of the plots below, we average the two macroscopically equivalent portions of the simulation cell to improve statistics.To accommodate the spatial scale needed to examine micron-scale mixing at the many-picosecond timescale, an aspect ratio of about 3:3:400 was used (17 nm × 17 nm × 2.3 μm); the simulation domain at three different scales is shown in Fig. 6.Thus, the simulations are quasi-1D, while still allowing for fully 3D collisions and correlations at the microscale, where the transverse (to the interface) length scale is roughly 75 interparticle spacings.To mimic separate, cold materials before the heating pulse, the plastic (CHO) ablator had a mass density of 1.05 g=cm 3 and consisted of 42.3% C (15,764,054 particles) and 57.2% H (21,293,565 particles), with a small dopant of 0.5% O (187,981 particles).The fuel consisted of approximately 50% D (10,166,786 particles) and 50% T (10,187,614 particles), with a mass density of 0.25 g=cm 3 .Thus, a total of 57.6 million particles were used in the entire simulation.These conditions were chosen to be consistent with typical ICF capsules [64].Each region was initialized as separate body-centered cubic (BCC) lattices that were brought together into the main simulation domain.While real DT ice and CHO plastic are not truly distributed in BCC lattices, our current electronic-structure method is not robust enough to describe the molecular or polymer states of these materials; however, the matter in these simulations heats very quickly away from this initial state.

Input summary for the MOD-MD model
The nuclear equations of motion in Eq. ( 1) are evolved using MD, given the total force on each nucleus, which arises from the three contributions of other nuclei, the slow electronic structure, and the fast Langevin forces.This is the dominant contribution to the computational cost of simulating the MOD-MD model.Because of the typical screening lengths of the system of interest, a neighbor-table algorithm was used to truncate the (infinite but convergent) force contributions from the nuclei.The time discretization used a standard velocity-Verlet integrator modified for a Langevin process [65], with a time step of Δt ¼ 0.01 fs.All simulations were carried out for 10-15 ps (10 6 -1.5 × 10 6 time steps) using 2 16 MPI ranks (4 MPI ranks per core, 16 cores per node, 1024 nodes).
Some specific features of the MOD-MD model made the simulations far more efficient than previous models that directly solve for the total electronic density with a single 3D grid that must resolve both macroscopic gradients and atomic-scale gradients.Because interface problems (and others, such as planar shocks) are one dimensional at the macroscale, we solve for the macroscale electronic structure on a one-dimensional grid perpendicular to the initial interface.The time associated with the electronic structure solution with 2000 grid points was negligible compared with that required by the MD integrator, as was the on-thefly determination of the ionic charge states.For this reason, we did not explore nonuniform grids, which the interface geometry lends itself to.To explore different heating FIG. 6. Simulation setup.Several views of the simulation cell are shown at different scales.Note that the taper that appears in the images results from the viewing perspective; the main cell is cuboid.The periodic cell contains 57.6 million particles total in a cell with an aspect ratio of 3:3:400, as seen in the lower image; approximately half of the total cell is shown, and the cross denotes the center of the cell.Atomic-scale mixing at the plastic-fuel interface is readily seen in the upper-right image, which is shown at 300 fs for the slowheating case.A color code for each species is shown in the legend.profiles, we did not use the electron-temperature equations in the form (11); rather, we prescribed the electron temperature to emulate an external source of energy.Both charged particles and photons primarily deposit their energy into electrons, and we model that process through two prescribed T e ðtÞ profiles: instantaneous heating to 50 eV and a linear ramp-up to 50 eV over 10 ps.The purpose of having two heating rates is to understand how kinetic effects are impacted by the timescale of the heating process.Through the Langevin model, the ions heat on a delayed time scale corresponding to γ −1 ei , as described in Sec.II.The results of the macroscale electronic structure calculation are shown in Fig. 7, where the ionic, electronic, and total charge densities are shown for two different times and for both heating rates.The specific parameters of the simulation are listed in Table II.
Ablator-fuel mixing was simulated to 10 ps, at which time the periodic boundary conditions interfered with the reality of the model.The results for the two heating rates at early and late times are shown in Fig. 8.In this figure, and the ones that follow, a mild filter [66] was applied to the data to reduce large fluctuations in cells with very small particle numbers.The location of the initial interface at z ¼ 0 is shown as the gray, vertical line.We see that the pressure in the CHO plastic exceeds that of the DT fuel, and the interface moves toward the DT, creating a density pulse that propagates into the fuel.As expected, this phenomenon is more pronounced when the heating is faster, as seen in the lower row in the figure.Interestingly, the penetration of plastic species into the fuel region is sensitive to the heating rate.With the lower heating rate (top row), the C species precedes the H and the O species into the fuel; conversely, for instantaneous heating, the H greatly precedes the C and O species, which penetrate at roughly the same rate.However, at late times, the lighter H from the plastic region has penetrated deepest into the fuel.In the case of rapid heating, the H jets into the fuel to distances that greatly exceed those for the other species.This hydrogen jetting will be explored further below.Species separation of this kind must be modeled through a mixture-transport model [63] that is applicable to warm dense matter conditions.Note that the D and T densities remain mostly locked to each other.
The MOD-MD approach enables simulations at the hydrodynamic scale, allowing us to explore the validity of approximations made within specific hydrodynamic models.For example, the mixing profiles in Fig. 8 show the importance of using an extended hydrodynamics model for modeling interfaces in ICF targets.While multispecies transport can be handled in terms of multiple continuity FIG. 7. Charge densities.The ionic (upper, blue line), electronic (lower, green line), and total (middle, red line) charge densities are shown, denoted as q i ≡ Z i en i , q e ≡ −en e , and the sum q tot , respectively.The electronic charge density results from the macroscale solver and is used to construct the mesoscale interaction potentials.The full (unfolded) cell is shown with both interfaces; the locations of the original interfaces are denoted by vertical grey lines.The DT fuel is in the center region between the grey lines, whereas the plastic ablator material is in the left and right regions that adjoin the boundaries.The top and bottom rows correspond to slow and rapid heating, respectively, and the left and right columns show early and late (300 fs and 10 ps) times, respectively.Strong quasineutrality is seen in all cases.FIG. 8. Mixing.Interfacial mixing (as species density fields) is shown for the two different heating rates (slow, top row; instantaneous, bottom row) at early (first column) and late (second column) times.The gray vertical line denotes the original location of the interface (plastic on the left; fuel on the right), and the O density has been multiplied by 50 in each panel to make it visible.Note that the D and T densities remain locked together, whereas the C and H densities separate, with strong "jetting" of the H into the DT fuel region.

Ion heating
Ions "thermostated" by coupling with electrons using a Langevin algorithm.The ion-electron coupling γ ei was determined using a Skupsky model and varied spatially as a function of electron temperature, local electron density, ion charge, and mass.For these simulations, γ ei was in the range 0.3-1 ps −1 .
Force cutoff radius 5.5× local screening length, which varies spatially equations [63], the need for multiple momentum equations is less clear.To explore whether multiple species-separated momentum equations are needed, we show v z ðzÞ, the z component of the velocity, for all five species in Fig. 9.In this variable, the hydrogen jetting is very pronounced, especially in the case of instantaneous electron heating.As with the density profiles, the D and T velocity fields are locked together.However, all three of the plastic species have separate velocity profiles, suggesting the need for separate momentum equations for each species.In practice, this situation is extremely problematic because one does not typically have access to species-resolved mixture equations of state.The velocity profile of the H shows steady acceleration into the fuel, followed by a stopping region in which the fast H deposits its energy into the other ions (directly through the MD forces) and the electrons (through the Langevin model [67]).The accelerated ions, in turn, convert their kinetic energy into thermal energy, as shown in Fig. 10, generating local species temperatures that greatly exceed the prescribed electron temperature.The velocity fields in Fig. 9 show strong segregation by mass: The peaks of the velocity fields tend to appear in the order H, C, and then O.It should be noted that the "noise" in the curves is due to the sparsity of particles at the leading edge of the front.In other words, the fluctuations show the presence of only a few particles in a numerical cell.We explore this result further in Fig. 11; in this figure, we show the H temperatures, now separated into parallel and perpendicular (to the initial interface) components, and the total electric field (in red).An extremely small amount of anisotropy is observed in the H temperature, mostly in the leading edge of the temperature profile in the rapidly heated case, suggesting a very mild kinetic effect.Because the warm dense matter is very collisional, we do not expect strong anisotropies in such quantities; however, at much higher temperatures, temperature isotropy is not likely to occur.Strong electric fields are observed near the interface; these are responsible for the acceleration of the plastic species into the fuel region, and this acceleration leads to the heating seen in Fig. 10.The early-time peak electric field is twice as strong for the instantly heated case (note the change of scale).The electric field peak is located on the fuel side of the interface, consistent with the rapid movement of the interface because of the balance of pressure between the plastic and the fuel.Note that a small secondary peak appears to the right of the main peak at early times, and comparing with the profiles in Fig. 8, we see that this peak is associated with the density pulse driven into the fuel by pressure equilibration.
Finally, we can use MOD-MD to examine deviations from Maxwellian distributions, which precludes a key assumption of most hydrodynamic (i.e., Navier-Stokestype) models.Near equilibrium, the normalized distribution of each velocity component for a given species will be of the form 9. Velocity fields.The z component of the velocity field of each species is shown.Hydrogen jetting is reflected in the large values of v z ðzÞ for that species, relative to the others; the jetting is obviously more extreme with instantaneous heating.(In the lower-left panel, the H velocity field has been scaled by a factor of 0.5.)FIG. 10.Ion temperatures.Ion temperature profiles are shown for four cases, as in the previous figures.The hydrogen species attains temperatures that exceed the prescribed electron temperature by more than a factor of 2, which leads to slight heating of the DT fuel species away from the interface.This heating results from stopping the fast plastic species in the fuel region.FIG.11.Hydrogen temperatures and electric field.The electric field (red), which has been multiplied by a factor of 500 in each panel, and the parallel (blue) and perpendicular (green) temperatures are shown.A strong electric field is present near the location of the interface, which is slightly shifted from its original location (see Fig. 8).We see that the field begins to diminish by 10 ps.However, before that happens, the hydrogen is accelerated, as seen in Fig. (9), and begins to deposit its kinetic energy in the fuel, heating the fuel above the target temperature of 50 eV.Note that there are subtle, but insignificant, differences in the parallel and perpendicular temperatures, suggesting rapid thermal equilibration.
The associated central moments can then be readily calculated as This pattern reveals that the moments of a Gaussian are connected and that specific ratios of powers of the moments are constant; these findings yield a metric for quantifying how "kinetic" a system is [68].For example, one can define the quantity If the system is near equilibrium, then we have G k ¼ 1 for any integer k.Therefore, deviations from unity in G k reveal non-Maxwellian behavior.We have shown G 2 in Fig. 12 for H using the v z velocity component.Despite the high collisionality at these temperatures and densities, G 2 exhibits clear deviation from unity, with an overall trend below unity in the tail region in which the H is accelerated and eventually stopped.Recent results using a kinetic model [69] show stronger kinetic effects for interfaces at much higher temperatures.

V. CONCLUSIONS
We have developed a multiscale computational method for heterogeneous charged systems (MOD-MD).Our method works by separating the length scales associated with atomic core electrons, interparticle forces, and global density fluctuations; the method is based on the observation that the Poisson equation can be separated into smoothly and rapidly varying contributions that can be handled with separate physics modules.This method differs from other types of multiscale models because it employs a single computational method (MD) but addresses different scales within the electronic structure calculation.Ions are treated directly through MD, with forces computed from this rapid electronic structure calculation, as well as from a Langevin prescription used to treat fluctuations in the electron density.Through the use of finite-temperature OFDFT, we are able to include quantum mechanics at all temperatures and densities in our method.This multiscale decomposition offers several advantages, including that it allows different symmetries and gridding at various scales, provides a modular treatment of the physics at each scale, and permits a more dynamic interaction with the electronic bath.What is potentially interesting about this idea generically is that it allows people to think about bottlenecks caused by dominant symmetries at different scales.Here, we have shown that by making the electronic structure multiscale, we can use the appropriate coordinate systems at three different length scales.
We have used our model to explore interspecies mixing at a heated interface relevant to ICF experiments.In this demonstration, we employed several specific modeling choices, including a TF ionization model for mixtures (see Appendix A), a macroscale TF model, and a dipole expansion at the mesoscale, to examine atomic mixing near a rapidly heated interface.These choices have permitted simulations with 57.6 million particles, allowing for micron-scale (in one direction) and tens-of-picoseconds scale simulations.We have performed simulations using two choices of electronic heating rates for an interface between DT ice and CHO plastic.Because of the high collisionality of matter in this system (small mean-free path), the simulations take place near the hydrodynamic regime (small Knudsen number), allowing our model to examine deviations from various forms of hydrodynamics (e.g., Euler, Navier-Stokes, Burnett).
From our results, we are able to explore hydrodynamic and kinetic phenomena.In our hydrodynamic studies, we extract the lowest-order hydrodynamic moments (density, velocity, and temperature), as well as transient electric fields in the interface region.While we observe very strong electric fields, which can potentially greatly enhance ionic transport and therefore mixing, the fields decay rapidly (over tens of femtoseconds), typically faster than any reasonable heating rate (tens of picoseconds).However, steady electric fields are observed when the heating rate is large enough to produce a shock moving away from the interface.Kinetics were examined through the two metrics of temperature anisotropy and the moments of the ionic velocity distribution; we find mild kinetic effects with these metrics, which is expected for matter of this collisionality.At higher temperatures (e.g., postshock), conditions are likely to experience more significant kinetic effects.
Our implementation of the MOD-MD model can be improved and applied to other applications in many ways.Perhaps the most important improvement would be better internal consistency, without breaking the multiscale method.For example, improvements in the mesoscale OFDFT model could be used to improve the Langevin model, including an improved (e.g., nonuniform) stopping power model and electronic heat capacity C e .A natural extension of the application presented here is to the kinetics of shocked plasmas [68,70].Because our multiscale model is modular, parts of the model can be adapted to different applications.For example, a Saha ionization approach [71] could be preferable to the TF model used here for lower-density plasmas that involve molecular species.Alternatively, the TF-AA result could be replaced with a higher-fidelity AA model [45].Because our choice of TF allowed us to exploit an existing fit, more work is warranted on rapid, high-fidelity, atomic-physics solvers.The most obvious extension is the development of a rapid, higherfidelity average atom, on-the-fly extraction of an improved pseudopotential from that model, and the use of an improved mesoscale model, such as a Perrot-like kinetic energy functional and a finite-temperature exchangecorrelation functional, for the valence electrons that surround those pseudo-ions.Such improvements would provide a better electronic equation of state in the cool-to-warm dense matter regime; it would also obviate the Teller nonbinding theorem.Similarly, different geometries would suggest different symmetries at the macroscale, perhaps different from those observed with the slab geometry used for our interface application.For some applications, an important weakness of the current model is its lack of electronic heat conduction, which can occur in mesoscale heterogeneous systems.An extension to include the electronic energy equation of the form used in the "two-temperature model" [33,[41][42][43] would allow the electronic bath to conduct heat across the simulation domain without the need for the continuity and momentum equations.radii a α ¼ ð4πn α =3Þ −1=3 .The mean ionization in each sphere is given functionally by where T e is the local electron temperature.The function ZðxÞ is the single-species mean-ionization function for given plasma conditions x [63].Note that Z α and T e are given, and there is one equation of type (A1) for each species; however, we do not yet know the subvolumes v α (or, equivalently, the ion-sphere radii a α ).
Next, we note that each species must self-consistently produce the same free-electron number density.Put another way, each atomic subvolume places pressure on the other subvolumes until the pressures equilibrate among the species.In the TF-AA model, the electronic pressure at the cell boundary is obtained from the electronic density at that boundary, and two cells are considered to be in pressure equilibrium if they have the same electronic density at the cell boundary.This is shown schematically in Fig. 13 and can be expressed mathematically as the two conditions which can be combined to give the volumetric constraint It is important to note that the first condition (A2) corresponds to an ionization model in which the free electrons are assumed to be uniformly distributed throughout the atomic subvolume [45]; otherwise, the relation between Z Ã α and the subvolume v α would require a separate calculation.
These nonlinear, algebraic equations are readily solved using any robust root solver with a reasonable initial guess.One simple approach to solving Eqs.(A1) and (A4) is to use a multidimensional Newton solver.The Jacobian for this set of equations is diagonal with the exception of the last row; thus, it can be inverted analytically.Hence, the multidimensional Newton iteration can be written as We find that the initial guess (corresponding to equal ionsphere volumes for each species) is very reliable.It is worth noting that this atomic physics scheme can be built into other computational models, such as hydrodynamic and kinetic models [69].Figure 14 shows an example of this multispecies atomic physics module [72].For a CHODT mixture with a fixed density, the mean ionization state is shown versus temperature.As expected, all species are partially ionized at low temperature, as a result of pressure ionization (captured through the finite cell volumes), and they tend toward their respective nuclear charges at very high temperature.

APPENDIX B: ASYMPTOTIC ANALYSIS
To find asymptotic solutions to the one-body problem (41), we first introduce the following transformations: where κ 2 i ¼ K 2 ðr i Þ and uð0Þ ¼ 1 for charge conservation.These transformations result in the following system to be solved: where r is the radius in spherical coordinates.Given the slow variations in hðrÞ, we can find asymptotically exact solutions by taking into account only the direction in which hðrÞ varies the most, which will be along the vector ∇K 2 ðr ¼ r i Þ.Without loss of generality, the one-body problem (B3) can then be rotated such that the z axis is collinear with this gradient, which transforms the operator L expressed here in spherical coordinates as follows: where θ is the inclination angle from the z axis.Obviously, this rotation will have to be inverted when the full manybody solution is constructed.To approximate solutions of Eq. (B4), we expand uðr; θÞ and hðr; θÞ as where ε is assumed to be small in the limit of lengthscale disparities between interparticle spacings and variations in K 2 ðzÞ, and the coefficients h n are known from a Taylor expansion of K 2 ðzÞ about z i given by h n ¼ ðd n =dz n ÞK 2 ðz i Þ=n!.At the lowest order (ε → 0), we have the system The only bounded, isotropic solution that satisfies the appropriate boundary conditions is u 0 ¼ 1, which in turn yields the Oð1Þ solution around a given ion: Using the definition of L in Eq. (B4), we can gather higher powers of ε to obtain the hierarchy of equations At this point, we can no longer assume isotropy, and it is useful to express our solutions in terms of Legendre polynomials, which have the defined properties For example, the first three Legendre polynomials are P 0 ðxÞ¼1, P 1 ðxÞ¼x, and P 2 ðxÞ¼ð3x 2 −1Þ=2.Furthermore, a solution of the form uðr; θÞ ¼ RðrÞP n ( cosðθÞ) is an eigenfunction in θ under the operator L given by the relation where the primes denote differentiation with respect to r.Before proceeding, we must examine the homogeneous solutions R H ðrÞ to the above operator; these solutions are of the form As the modified Bessel function of the second kind has small-argument behavior ffiffi ffi r p e κ i r K nþ1=2 ðκ i rÞ ∼ r −n , we must then have c 2 ≡ 0 for all n ≥ 1.Similarly, the modified Bessel function of the first kind has large-argument behavior ffiffi ffi r p e κ i r I nþ1=2 ðκrÞ ∼ e 2κ i r , so we must also have c 1 ≡ 0 for all n ≥ 1.Therefore, only the trivial homogeneous solution R H ≡ 0 is allowed; hence, we need to calculate only particular solutions of Eq. (B9).At OðεÞ, we have Seeking a solution of the form u 1 ðr; θÞ ¼ R 1 ðrÞP 1 ( cosðθÞ) yields the differential equation We have set ε to unity, as it was only being used to track orders.Note that the correction induces a small well in the potential at large r for h 1 z > 0; this effect is spurious and results from the truncation of the asymptotic series.As this attraction will occur only for roughly r ≳ ð4κ i =h 1 Þ 1=2 , it will be negligible because the force decays exponentially and because h 1 is a small quantity.For example, if κ i ¼ 1 and h 1 ¼ 0.1, then e −κ i r =r 2 < 4.5 × 10 −5 in this range.Of course, this artificial attraction can potentially yield misleading results as the assumptions of the expansions break down.To recover the general result, we rotate the coordinates back to the original reference frame through the transformations z → ∇K 2 ðr i Þ j∇K 2 ðr i Þj • r; h 1 → j∇K 2 ðr i Þj ðB18Þ to obtain the solution A similar result for a polarized pair interaction can be found in Ref. [73].

APPENDIX C: HIGHER-ORDER CORRECTIONS
To calculate the next-order correction to Eq. ( 43), we examine the Oðε 2 Þ problem, We write u 2 ðr; θÞ ¼ R 0 ðrÞP 0 ( cosðθÞ) þ R 2 ðrÞP 2 ( cosðθÞ) and obtain the system of ordinary differential equations which have the particular solutions Combining the above results yields the second-order perturbation which, in turn, can be combined with Eq. (B17) to obtain the following next-order correction to the potential: Using Eq. (C9), the force on the ith ion that is correct to where k is the unit vector along the z axis, and h ðiÞ 1;2 denote the particular expansion coefficients about r ¼ r i .The many-body contributions to the force F i can be expressed in terms of the transformed variables in Eq. (B1) as where uðr ij Þ and its derivatives are given by with z ij ¼ z i − z j .The first term in Eq. (C12) can be somewhat simplified by writing e −κ j r ij ∂u ∂z ðr ij Þ k; Lastly, the screening cloud corrected to Oðε 2 Þ is given by n ðiÞ e ðrÞ ¼ Z Ã i 4πr e −κ i r ½κ 2 i þ g 1 ðr; zÞ þ g 2 ðr; zÞ; ðC17Þ where the functions g 1;2 ðr; zÞ are defined as As before, the general solution in the original reference frame can be recovered using the transformations (B18), as well as where fx i g are the various components of r.
FIG. 2. MOD-MD algorithm.Schematic showing the concurrent steps at each length scale during each time step within the MOD-MD method.

FIG. 3 .
FIG.3.Convergence of K. Convergence of the screening function KðzÞ with increasing number of bins.Profiles of KðzÞ are shown in the top row, while convergence of the solution is presented in the bottom row in terms of the relative error E c ¼ kK − K c k=kK c k, where K c is the converged function, and kfðzÞk ≡ R R dzjfðzÞj 2 .An example system of a planar interface is explored (see Sec. IV for more details), where the left column is the initial (t ¼ 0 ps) sharp interface and the right column is a late-time (t ¼ 15 ps) diffuse interface.

FIG. 13 .
FIG.13.Multispecies average-atom algorithm.Once an initial guess is chosen, Eqs.(A1) and (A4) are solved iteratively.Here, each atomic subvolume places pressure on the other subvolumes until the pressure equilibrates among the species.In the TF model, the pressure is obtained from the electronic density at the cell boundary, and pressure equilibrium is equivalent to exhibiting the same electronic density at the cell boundary.

R H ðrÞ ¼ ffiffi ffi r p e κr ½c 1 I
nþ1=2 ðκ i rÞ þ c 2 K nþ1=2 ðκ i rÞ; ðB13Þ where c 1;2 are arbitrary constants, and I ν and K ν are modified Bessel functions of the first and second kind.The perturbations u n ðr; θÞ, for n ≥ 1, must have the properties lim r→0 u n ðr; θÞ ¼ 0; lim r→∞ Z i r e −κ i r u n ðr; θÞ ¼ 0 ðB14Þ to satisfy the boundary conditions of the original problem.

TABLE II .
A summary of the computational choices made for simulating interfacial mixing.Ablator particle number N C ¼ 15 764 054, N H ¼ 21 293 565 and N O ¼ 187 981 initially distributed on a 60 × 60 × 10 346 lattice Fuel particle number N D ¼ 10 166 786 and N T ¼ 10 187 614 initially distributed on a 60 × 60 × 5654 lattice Step function at t ¼ 0 (0 → 50 eV) Case 2: Ramp from 0 eV at t ¼ 0 to 50 eV at t ¼ 10 ps Constant of 50 eV after t ¼ 10 ps has the particular solution R 1 ðrÞ ¼ −h 1 r 2 =4κ i .Combined with the Oð1Þ solution, we now have the corrected electrostatic potential