Ion Modes in Dense Ionized Plasmas through Non-Adiabatic Molecular Dynamics

We perform non-adiabatic simulations of warm dense aluminum based on the electron-force field (EFF) variant of wave-packet molecular dynamics. Comparison of the static ion-ion structure factor with density functional theory is used to validate the technique across a range of temperatures and densities spanning the warm dense matter regime. Differences in the dynamic structure factor and dispersion relation between adiabatic and non-adiabatic techniques suggest that the explicit inclusion of electrons is necessary to fully capture the low frequency dynamics of the response function.

The warm dense matter (WDM) regime defines a dense plasma state where strongly coupled classical ions coexist with partially or fully degenerate electrons [1]. It is a complex state of matter where multi-body particle correlations and quantum effects play an important role in determining the overall structure and equation-of-state [2]. A complete description of WDM is important for describing many physical phenomena, ranging from phase transitions within the interior of large astrophysical objects [3] to temperature relaxation rates during the internal processes of inertial confinement fusion [4].
Strong ion-ion coupling combined with the quantum behavior of the electron fluid make simulation and modeling challenging. One approach that is able to fully capture the strongly coupled behavior of the ions has been molecular dynamics (MD), i.e. a fully atomistic simulation of the ion dynamics that results from numerical integration of the equations of motion [5]. As MD simulations explicitly calculate the ion trajectories they are able to provide transport properties, such as viscosity and thermal diffusivity [6], acoustic properties, such as the sound speed [7,8], and thermodynamic variables, including the equation-of-state [5,7]. The accuracy and predictive power of such simulations is encompassed within the calculation of the interatomic potential, with increased accuracy coming at the expense of computational cost. Simulations using analytic model potentials are able to model the largest systems [9]; however, potentials calculated with orbital-free [7] (OF) or Kohn-Sham [10,11] (KS) density functional theory (DFT) have had the most success when compared to available experimental results [12,13].
MD simulations of dense plasmas typically employ the Born-Oppenheimer (BO) approximation, also known as the adiabatic approximation, usually justified through the disparate energy scales of the electron and ion motion [14]. This approximation is the cornerstone of almost all classical [9] and ab-initio [7,10,11] atomistic simulations; the electrons are assumed to instantaneously adjust to the ion fields, while the ions themselves are confined to a single adiabatic surface [15]. The BO approximation thus neglects important details of the interaction and, of particular relevance to high energy density matter, the dynamics of electron-ion collisions are ignored [16]. Although justified for many equilibrium properties, like the equation-of-state [8,17], the adiabatic treatment prohibits direct energy transfer between electrons and ions and is therefore problematic for the modelling of transport properties and non-equilibrium matter [18].
Recent work has investigated the applicability of a Langevin thermostat to capture the effect of electron-ion interactions on the ion dynamics [8,17,19]. In this phenomenological approach, an additional stochastic Gaussian force is added to the equations of motion with a single collision frequency, or friction parameter, that must be determined a-priori.
A number of classical (Rayleigh scattering [20]) and quantum (Born approximation [21]) models exist, but their applicability in the WDM state is unknown.
The inclusion of these collisions has been shown to have a profound effect on the ion dynamics [8], including a strongly decreased adiabatic sound-speed and increased diffusive mode around ω = 0.
In order to explicitly go beyond the BO approximation, the sytem must be treated on an electronic timescale while correctly modelling the quantum nature of the electrons. In quantum chemistry, time-dependent DFT is one approach that attempts this. However, for systems made up of the thousands of the ions necessary to correctly model the low-frequency ion dynamics, this method remains too computationally intensive [22]. More recently, a technique based on the Bohm theory of quantum mechanics has been implemented. This method propagates an array of thermally-coupled N-trajectories interacting with a thermally-averaged, linearized Bohm potential [15]. An adiabatic sound-speed in agreement with the Langevin approach was found, with a collision frequency between the quantum and classical limits. Interestingly, the dynamics were found to have a significantly reduced diffusive mode when compared to the Langevin technique.
Here, we utilize an alternative method, Wave Packet Molecular Dynamics (WPMD). It is a time-dependent quantum mechanical technique able to simultaneously simulate (1) the propagation of the ions as classical point particles, and (2) the electrons as quantum mechanical entities [23][24][25][26][27]. The direct inclusion of electrons, and thus the effects of electron-ion interactions, means that WPMD intrinsically goes beyond the adiabatic approximation and is capable of calculating observables in quantum many-body systems [24].
In WPMD, each electron is represented as a wavepacket, a spatially localized complex function. It is common to represent a single electron wave-packet using a basis set constructed from a sum of M Gaussians. i.e., x is a normalization factor. The remaining terms are dynamical parameters: complex amplitude (c α ), position (r α ), momentum (p α ), Gaussian width (s α ), and conjugate width momentum (p sα ) [23]. These wave-packets uniquely define the quantum state of a single electron, with the total many-body wavefunction constructed from either a Hartree product or Slater determinant [24].
Equations of motion for the dynamical parameters are easily derived from variation of the time-dependent Schrodinger equation [28,29]. One of the great accomplishments of WPMD as a method of modelling plasmas is maintaining the same computational efficiency as in many classical methods [23,25,30].
Here we utilize a simple variation of WPMD, where electrons are described by a single, floating Gaussian wave-packet combined into a many-body wavefunction using a Hartree product [31][32][33]. In this case, the equations of motions take on a simple Hamilton form. This method, known as the Electron Force Field (EFF), includes exchange and correlation through an additional pairwise energy term to account for the lack of antisymmetrization. These terms, introduced between same-spin and opposite-spin electrons, are constructed with solely the kinetic energy penalty upon pairwise antisymmetrization, and a number of scaling parameters, matched to a set of molecular test structures [31,34].
In addition, for high-Z systems, a fixed size wavepacket is attached to the otherwise classical ions to represent the strongly bound core electrons. These effective core potentials (ECPs) provide not only an increased computation speed, as the high-frequency oscillations associated with these electrons do not need to be captured, but also increased accuracy as they may exhibit non-Gaussian shapes [34,35]. Due to the strong binding energy and high-frequency oscillations of these core electrons we do not expect them to play a large role in the low-frequency ion dynamics. EFF is able to simulate systems consisting of thousands of particles for many picoseconds. It has been applied to the study of material FIG. 1. The static structure factor for warm dense aluminum at 3.5 eV and 5.2 g/cm 3 , calculated with the nonadiabatic electron force field (EFF) method. For comparison, results from orbital-free density functional theory (OF-DFT) and Kohn-Sham density functional theory (KS-DFT) are included [8].
properties in extreme environments [32,[35][36][37], temperature relaxation rates in warm dense hydrogen [38], and wave-packet spreading in electron-nuclear scattering [23]. However, the low frequency ion modes in dense plasmas have not yet been examined. The ion modes are typically investigated through the spatiotemporal Fourier transform of the density-density auto-correlation function, or the dynamic structure factor (DSF). This common metric, used to compare different theoretical techniques, is related to the dynamic and thermodynamic properties of the plasma, but also accessible experimentally, allowing for critical benchmarking of models [12,39]. In particular, collective ion modes are a prominent feature, whose spectra serve as an important tool to validate theoretical predictions for dense matter [40]. For a system in thermodynamic equilibrium, the DSF gives its response to fluctuations of frequency ω and wave-vector k, and is defined as where N is the total number of particles, and . . . refers to an ensemble average. Here, ρ is the Fourier transform of the real space time-dependent density distribution. For a system of particles with coordinates r j (t) this is given by, The static structure factor (SSF) is the frequency integrated version of this quantity that relates purely to the structure of the system.  Figure 1 shows the SSF calculated with EFF for warm dense aluminum at a temperature of 3.5 eV and a density of 5.2 g/cm 3 . This thermodynamic condition allows for direct comparison with a number of previously published results, both adiabatic and non-adiabatic. All EFF simulations were performed using the LAMMPS software package [41] using 216 ions and 648 electrons. Each simulation was run for 20 ps and utilized a 1 attosecond timestep. Additional details of the simulation, including the equilibration procedure, are provided in the supplemental material [42]. In figure 1 we have additionally plotted the SSFs calculated from the more computationally intensive OF-DFT and KS-DFT MD. Although adiabatic, such results can be used to validate EFF throughout the WDM regime, as they have been shown to accurately model the SSF [8]. Excellent agreement is found between the SSF calculated with EFF and with OF-DFT. Both results take similar computational times [42]; however, EFF explicitly includes the electron dynamics. The KS-DFT results are currently the most accurate available, but extremely computationally intensive. They exhibit an increased peak height around k = 2 a −1 B , suggesting that neither EFF or OF-DFT completely captures the strong ion-ion correlations in the system.
To validate the EFF technique across the entire WDM regime, we ran multiple simulations that spanned a range of thermodynamic conditions. The SSF for each condition was calculated using both EFF and OF-DFT, with the L1-difference between the two plotted in figure 2. The OF-DFT simulations used 108 ions and were performed in the software package ABINIT [43,44]. Each simulation was run for 5 ps and utilized a 0.5 fs timestep. Additional details can be found in the in the supplemen- tal material [42]. These results suggest an area of validity in the low-density/high-temperature region. The cross marks the thermodynamic condition plotted in figure 1 and represents the best agreement between the two techniques. The biggest disagreement between EFF and OF-DFT is found in the low-temperature/high-density regime; this difference is attributed to crystallization in the OF-DFT simulations that is not apparent with EFF. However, with neither KS-DFT simulations or experimental results available in this regime, the validity of both techniques may be questioned. Two modifications to the EFF technique were required to achieve the best agreement with the OF-DFT simulations for warm dense aluminum. The first, as also found in Refs. [25,38], is the inclusion of a harmonic potential to constrain the width of the Gaussian wave-packet. Such a constraint is needed to prevent the size of the wave-packet increasing without limit and leading to unphysical diminished electron-ion and electron-electron interactions. While such expansion is itself not unphysical, a wave-packet will naturally spread when not confined by a potential; the limited basis set prevents localization around ions. It is this lack of localization that leads to unphysically low interaction.
The second change made was to alter the width of the ECP Gaussian, the Gaussian attached to the ions representing the 10 inner shell electrons. For aluminum, EFF uses a radius of 1.66 a −1 B , determined from the ground state of face-centered cubic bulk aluminum [34]. We found that this number needed to be decreased to 1.285 a −1 B , much closer to the radius of an Al 3+ ion [45]. Colormaps showing the L1-difference without the harmonic constraint or the decreased ECP width are given in the supplemental material [46].
We now turn our attention to the dynamic plasma properties, which are expected to vary more with the inclusion of non-adiabatic effects. Figure 3 shows the DSF calculated with EFF for aluminum at 3.5 eV and 5.2 g/cm 3 . The EFF results are well matched by the much more computationally intensive KS-DFT results, agreeing well at the highest and lowest k values tested. However, differences are apparent in the shape of the spectra at the intermediate spatial scale (i.e., k = 0.76 a −1 B ). Interestingly, as with the Bohmian Mechanics technique, we see a strongly reduced diffusive mode in the DSF around ω = 0.
The peak position is excellently matched across all kvalues, as exemplified in the dispersion relation shown in Figure 4. We find excellent agreement between the EFF method and the KS-DFT results. Both techniques agree on the adiabatic sound speed (given by the low-k slope of the dispersion relation) as well as their prediction for the onset and strength of negative dispersion. For comparison, we have included in Figure 4 the results from OF-DFT, with a Langevin thermostat, and those calculated with from the Bohmian Mechanics technique. While the Bohmian mechanics results are closer to those with a high collision-frequency, highlighting the importance of non-adiabatic techniques, the EFF results suggest that the electrons may have less of an effect on the ion dynamics than previously thought.
Despite its success, the comparison with OF-DFT suggests that EFF breaks down at the extremum of temperature and density. Even at our trial condition of 3.5 eV and 5.2 g/cm 3 , there is a lack of correlation in the SSF when compared to that of KS-DFT. Furthermore, the discrepancy in shape of the DSF at k = 0.76 a −1 B between EFF and KS-DFT suggests there may be differences in transport properties between the two models. To address these discrepancies, we suggest improvements to the EFF formalism that tackle three fundamental approximations. These are (1) a limited basis set consisting of a single Gaussian wave-packet; (2) neglect of the electron-electron and electron-ion components of the exchange energy; and (3) applicability of the constants used in the correlation functional to the WDM regime. The first of these limits wave-packet breakup and has previously been identified as a flaw in WPMD [23]. Work to expand the region of applicability of EFF should focus on these three limitations of the current EFF model. Ultimately, the we have shown that the EFF technique is able to closely match the static properties of warm dense aluminum, obtained with the more robust density functional theory, across a large portion of the WDM phase space. This is achieved with less computational effort and full treatment of the electron dynamics. The effect of electron-ion correlations on the ion modes were investigated and found to have less effect than when utilizing Bohmian mechanics [15], with a dispersion relation in closer agreement to the adiabatic KS-DFT results. These results directly address concerns over WPMD, as an ab-initio many-body method for dense plasmas [24], by providing corroboration with other models.
Overall, there is still much we do not understand about the ion dynamics in warm dense matter. The results of the Langevin thermostat, Bohmian mechanics technique, and now WPMD techniques play an important role in our understanding of the ion modes, sound-speeds, and transport coefficients in dense plasmas. Of course, experimental verification along with further computational investigation is essential.

I. SIMULATION PARAMETERS AND EQUILIBRATION
All EFF simulations were performed using the LAMMPS open-source software. Standard periodic boundary conditions were used with N i = 216 aluminum ions and N e = 648 electrons. These simulations were run with a time-step of 1 attosecond for a total of 20 million time-steps, giving an overall simulation time of 20 ps. The inner core electrons were modeled using an effective core potential with a radius of 1.285 a −1 B . All other parameters utilized the default EFF values. Before the calculation of any structure factors, the simulation was equilibrated for 90 ps using a Nose-Hoover and Langevin thermostat until the electrons and ions were in equilibrium. After equilibration, each simulation took 9 hours on 32 CPUs.
Temperature equilibration in the simulation employed a Langevin thermostat coupled only to the ions, with the temperature of the electron sub-system allowed to equilibrate purely through electron-ion collisions. The built in EFF temperature control, which couples to both electrons and ions, was found to lead to non-equilibrium velocity distributions. Figure S1 shows the evolution of the energy in each sub-system throughout the T=0.5 eV simulation, with similar behaviour found at all temperatures and densities. At 90 ps, the thermostat is removed, and the static and dynamic structure factors calculated. It should be noted that the electron energy approaches 4 2 N e T = 23.8 Ha and the ion energy 3 2 N i T = 5.95 Ha. In the WPMD approximation, the Gaussian width is taken to be a classical fourth degree of freedom for each electron [24].
The OF-DFT simulations used to generate Figure 2 in the main document were performed in the open-source software ABINIT using the Thomas-Fermi module and the exchange-correlation terms in local density approximation of Perdew-Zunger-Ceperley-Alder [47]. For the interaction between ions and electrons, a bulk-local pseudopotential [48]  is used in which the ten inner electrons are frozen, while the three valence electrons are explicitly treated. Simulations employed a 108-ion cubic supercell with periodic boundary conditions. They ran for 5 ps and used a 0.5 fs timestep for 5 ps. A plane-wave energy cut-off of 30 eV was used with 880 bands. To control the temperature, a Nose-Hoover thermostat [49] with an inertia parameter corresponding to a temperature oscillation period of over 100 time-steps was used. This ensures weak coupling to the heat bath and a correctly sampled canonical distribution. After equilibration, each OF-DFT simulation took 5 hours on 16 CPUs.

II. MODIFICATIONS TO THE ELECTRON FORCE FIELD MODEL
Two modifications to the EFF technique, as described in the literature, were required to achieve the best agreement between EFF and OF-DFT for warm dense aluminum. The first is the inclusion of a harmonic potential to constrain the width of the Gaussian wave-packet. Secondly, the width of the effective core potential (ECP) was decreased from 1.66 a −1 B to 1.285 a −1 B , much closer to the radius of an Al 3+ ion. Figure S2 demonstrates the enhanced improvement between EFF and OF-DFT with these changes. Figure S2a uses the parameters implemented in the main manuscript, while Figure S2d contains the original EFF parameters. All of the SSFs that were used to produce these plots are shown in Figures S3-S6.