COMPUTATIONAL STUDY OF STRONGLY COUPLED CHARGED PARTICLE SYSTEMS By Gautham Dharuman A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Computational Mathematics, Science and Engineering – Doctor of Philosophy Electrical Engineering – Dual Major 2018 COMPUTATIONAL STUDY OF STRONGLY COUPLED CHARGED PARTICLE SYSTEMS ABSTRACT By Gautham Dharuman Exciting experiments in ultracold neutral plasmas, laser-matter interaction, charged particle stop- ping, mixing under extreme conditions etc., at academic facilities or at even larger facilities such as the National Ignition Facility, Z machine or the Linac Coherent Light Source, have necessitated the need for models that can simulate these systems at large length- and time-scales. This the- sis summarizes my research work, falling within the category of computational plasma physics, aimed at three aspects: effective quantum potentials based method for non-equilibrium quantum electron dynamics at scale, efficient force calculation method for molecular dynamics simulation with screened Coulomb interactions, and an avenue based on compressed gases for creation of laboratory-scale tunable strongly coupled plasmas as a platform for understanding large-scale ex- periments. Effective classical dynamics provide a potentially powerful avenue for modeling large-scale dy- namical quantum systems. We have examined the accuracy of a Hamiltonian-based approach that employs effective momentum-dependent potentials (MDPs) within a molecular-dynamics frame- work through studies of atomic ground states, excited states, ionization energies and scattering properties of continuum states. Working exclusively with the Kirschbaum-Wilets (KW) formula- tion with empirical MDPs [C. L. Kirschbaum and L. Wilets, PRA 21, 834 (1980)], leads to very accurate ground-state energies for several elements (e.g., N, F, Ne, Al, S, Ar and Ca) relative to Hartree-Fock values. The KW MDP parameters obtained are found to be correlated, thereby re- vealing some degree of transferability in the empirically determined parameters. We have studied excited-state orbits of electron-ion pair to analyze the consequences of the MDP on the classical Coulomb catastrophe. From the ground-state energies, we find that the experimental first- and second-ionization energies are fairly well predicted. Finally, electron-ion scattering was examined by comparing the predicted momentum transfer cross section to a semi-classical phase-shift cal- culation; optimizing the MDP parameters for the scattering process yielded rather poor results, suggesting a limitation of the use of the KW MDPs for plasmas. Efficient force calculation methods are needed for molecular dynamics simulation with medium- range interactions. Such interactions occur in a wide range of systems, including charged-particle systems with varying screening lengths. We generalize the Ewald method to charged systems de- scribed by interactions involving an arbitrary dielectric response function ε(kkk). We provide an error estimate and optimize the generalization to find the break-even parameters that separate a neighbor list-only algorithm from the particle-particle particle-mesh (PPPM) algorithm. We examine the implications of different choices of the screening length for the computational cost of computing the dynamic structure factor. We then use our new method in molecular dynamics simulations to compute the dynamic structure factor for a model plasma system and examine the wave-dispersion properties of this system. Laboratory-scale non-ideal plasmas with controllable properties over a wide range of densi- ties below solid density are needed for understanding large-scale plasma experiments. Based on a suite of molecular dynamics simulations, we propose a general paradigm for producing such con- trollable non-ideal plasmas. We simulated the formation of non-equilibrium plasmas from pho- toionized, cool gases that are spatially precorrelated through neutral-neutral interactions that are important at moderate to high pressures. We examined the plasma-formation process over orders- of-magnitude variations in the initial gas pressure to characterize variations in several physical properties, including Coulomb collisional rates, partial pressures, screening strengths, continuum lowering, interspecies Coulomb coupling, electron degeneracy and ionization states. We find that variations in the initial gas pressure lead to controllable variations in a wide range of plasma prop- erties, including the equation of state, collisional processes, atomic processes and basic plasma properties (coupling, screening and degeneracy). This paradigm has significant advantages over solid-density experiments because the collisional, collective and recombination timescales are re- duced by a factor of 3− 10, potentially broadening the efficacy of diagnostics. The paradigm also has advantages over ultracold plasma experiments because the trapping and cooling phases are avoided. Copyright by GAUTHAM DHARUMAN 2018 Dedicated to the eager minds. vi ACKNOWLEDGEMENTS I wish to express my deepest gratitude to my adviser Prof. Michael S. Murillo for your exceptional mentoring and for the years of patience and support you have given me. It is through your invalu- able guidance, constant encouragement and endless patience that I have been able to complete my journey through the challenging stages of graduate school. The information and skills you taught me, and motivated me to learn, have advanced me academically and professionally. I cannot thank you enough for your mentorship and friendship. I want to thank my co-adviser Prof. John Verboncoeur, for stimulating my interest in compu- tational plasma physics through his lectures that I attended during my Masters program at MSU and for your patience and the knowledge sharing and advice you shared with me in class or con- versation. I want to thank Prof. Andrew Christlieb, for his kind advice, knowledge sharing and encouragement. Together with my other committee members, Profs. Philip Duxbury, Vladimir Zelevinsky and Shanker Balasubramaniam, I want to thank you all for your time and patience with my dissertation and writing process. Finally, I want to thank my family for their unrelenting support that has pushed me, calmed me, and made all of this possible. Particularly, to my father, Dharuman, thank you for standing by me in hard times and encouraging me through the tired stages. To my mother and sister, thank you for your endless love, support, and motivation for me to try my hardest. Finally, to my aunt and uncle, thank you for your faith in me to grow as a researcher and for your endless love and support. vii TABLE OF CONTENTS LIST OF TABLES . LIST OF FIGURES . . . . . . . . . . . KEY TO ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x xi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 CHAPTER 2 ATOMIC BOUND STATE AND SCATTERING PROPERTIES OF EF- . . . . . . . 2.1 Formulation . 2.2 Ground state energies . FECTIVE MOMENTUM-DEPENDENT POTENTIALS . . . . . . . . . . 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.1 Constants of motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 . . . . . . . . . . . . . . . . . . . . . . . . 19 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3 Excited state orbits . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4 First and second ionization energies 2.5 Scattering properties of free electron states . . . . . . . . . . . . . . . . . . . . . . 24 2.6 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.2.1 Correlated core sizes (εH, εP) . . . . . INTERACTIONS . CHAPTER 3 A GENERALIZED EWALD DECOMPOSITION FOR SCREENED COULOMB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1 Linearly screened Coulomb systems . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2 Ewald decomposition for arbitrary dielectric response functions . . . . . . . . . . . 35 3.3 Error analysis and timing studies . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4 Application: Dynamic Structure Factor . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4.1 Numerical computation of the dynamic structure factor from MD . . . . . . 47 3.4.2 Kernel smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.4.3 Comparison of the theoretical and MD dispersion relations . . . . . . . . . 52 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.5 Summary and conclusions CHAPTER 4 CONTROLLED NON-IDEAL PLASMAS FROM PHOTOIONIZED COM- PRESSED GASES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.1 Mitigating DIH by photoionizing a precorrelated neutral state at room temperature . 58 4.2 Impact of precorrelation on the initial microfield distribution . . . . . . . . . . . . 60 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.3 Effective coupling parameter 4.3.1 Coulomb Coupling Estimation . . . . . . . . . . . . . . . . . . . . . . . . 62 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.4.1 Ionization potential depression . . . . . . . . . . . . . . . . . . . . . . . . 66 4.4.2 Transport phenomena through Coulomb logarithms . . . . . . . . . . . . . 68 4.4.3 Multi-temperature equation of state model . . . . . . . . . . . . . . . . . . 68 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.5 Summary and conclusions 4.4 Application of PIPPs . . viii CHAPTER 5 CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 APPENDICES . APPENDIX A APPENDIX B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 APPENDIX TO CHAPTER 3 . . . . . . . . . . . . . . . . . . . . . 76 APPENDIX TO CHAPTER 4 . . . . . . . . . . . . . . . . . . . . . 85 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 ix LIST OF TABLES Table 3.1 Comparison of times per time-step for force calculations using the PPPM and LCL algorithms for Yukawa interaction, for a range of κ and the number of unit charges (N). The entries correspond to an error of approximately 10−5 (e2/a2 i ). All times are in seconds. α is the Ewald splitting parameter, rc is the cutoff radius for the real-space calculation, M is the number of grid points along each direction in Fourier-space, and rD c is the cutoff radius for LCL. For all Fourier-space calculations, B-spline of order p = 6 was found to be optimal. TR and TF are the times per time-step of real-space and Fourier- space calculations, respectively. TPPPM refers to the total time per time-step for PPPM and is given by TPPPM = TR + TF . TLCL is the time per time-step for LCL. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 x LIST OF FIGURES Figure 1.1 Temperature-density plot showing different kinds of plasmas (For reference see: http://plasma.szfki.kfki.hu/ zoli/LFO_webpage/LFO_2011_EN.html). Con- tours of Γ = 0.1,1,10 indicate regions of strongly coupled and weakly cou- pled plasmas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 1.2 Timeline summarizing some of the essential developments in modeling non- equilibrium quantum electron dynamics at large length- and time-scales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (cid:34) Figure 1.3 Timeline summarizing some of the essential developments in force and po- tential energy calculation methods for MD simulations involving Coulomb and screened Coulomb interactions. One class of screened Coulomb interac- tions is the Yukawa interaction. 4αr2 eα 1−(rp/ε)4(cid:105)(cid:35) (cid:104) Figure 2.1 Plots of r2 (in atomic units) where α = 1 and ε = 1 (a) and α = 1 and ε = 2 (b) show that for rp (cid:28) ε, the potential becomes very repulsive. Since V H and V P have similar functional forms, common notations are used to denote the potentials and their corresponding variables: V (r, p) denotes V H or V P, α denotes αH or αP, ε denotes εH or εP, r denotes ri or ri j and p denotes pi or pi j. Therefore, the highly repulsive behavior of V (r, p) as rp (cid:28) ε enforces the Heisenberg and Pauli principles within the classical framework. V (r, p) = ε2 2 5 5 . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Figure 2.2 Projected surface of ∆HF (%) = 100|(E − EHF )/EHF| in (εH ,εP) plane for Al (a) and Ca (b) showing the region of minima. (εH,εP) points (white dots) corresponding to ∆HF (cid:46) 4% follow a curve indicating that there is a correla- tion between εH and εP values that result in ground state energies which are in very good agreement with HF values. The correlation between the param- eters reveal some degree of transferability to ground state energies of other elements that are not tested. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Figure 2.3 (a) Relative error ∆HF (%) of MDP prediction of ground state energies with respect to HF values (green circles connected by dashed green line) for N, Ne, Al, Ar and Ca show that the KW formulation is an excellent model for (b) Correlated (εH,εP) ground state energies (lines are to guide the eye). points (circles in shades of red from light to dark corresponding to increasing atomic number) extracted from the ∆HF surfaces. Also shown are curve fits (solid curves in shades of red increasing from light to dark corresponding to increasing atomic number) using εP = Aε2 H + BεH +C for N, Ne, Al, Ar and Ca. For all the elements considered, the curve fits match well with the points extracted from the corresponding ∆HF surfaces suggesting a possible transferability with respect to the atomic number. . . . . . . . . . . . . . . . . 16 xi Figure 2.4 Figure 2.5 (a) Coefficients A, B and C of the fit εP = Aε2 H + BεH +C vary as a function of atomic number (Z). The pattern in the data points of A (green circles), B (black star) and C (red squares) corresponding to N, Ne, Al, Ar and Ca are captured well by 4th degree polynomial fits for A (solid green curve), B (dot- ted black curve) and C (dashed red curve). (b) Correlated (εH,εP) curves for F (blue dashed) and S (blue solid) computed with A, B and C interpolated us- ing their corresponding 4th degree polynomial fits. Ground state energies for F and S from their correlated (εH,εP) are in good agreement with HF values (with a maximum ∆HF of ∼10% and ∼5% for F and S, respectively). This confirms the transferability of the trained (εH,εP) to ground state energies of elements that were not included in the training. . . . . . . . . . . . . . . . . . . 17 (a) Electron trajectory around ion (black dot) corresponding to different ini- tial conditions: Kepler-like motion from t = 0 to t = 50 a.u. (blue curve), zig-zag motion from t = 0 to t = 50 a.u. (dark green curve), zig-zag motion from t = 50 to t = 231 a.u. (light green curve) and reciprocating motion from t = 0 to t = 50 a.u. (red curve). Reciprocating motion corresponds to zero initial angular momentum which would result in an unstable trajectory in the absence of Heisenberg MDP. (b) Magnitude of position (r) and magnitude of momentum (p) corresponding to the trajectories in (a) from t = 0 to t = 50 a.u. are superimposed on the Hamiltonian surface. (c) Region 1 of (b) mag- nified to show the band-like structure in (r, p) dynamics. (d) Region 2 of (b) magnified to show a similar band-like structure in (r, p) dynamics. . . . . . . . 20 Figure 2.6 For N, F, Ne, Al, S, Ar and Ca, first and second ionization energies were com- puted using MDPs with the correlated (εH, εP) optimized to give accurate neutral ground-state energies. MDP prediction of first ionization energies (red curve) are in good agreement with NIST data (dashed blue curve); sec- ond ionization energies using MDPs (green curve) yield mixed results com- pared to NIST data (dashed cyan curve). Therefore, the parameters trained on neutral ground state energies transfer fairly well to the prediction of first and second ionization energies with some outliers. . . . . . . . . . . . . . . . . 21 Figure 2.7 Comparison of the classical and quantum MTCS for Z = 1, κB = 1, and E = 10 to 100 a.u. The classical MTCS using semi-analytic method (blue dashed) and the trajectory method (cyan dashed) are in very good agree- ment. They also match with the asymptotic limit of the classical MTCS (red dashed) given by Eq. 2.28. The quantum MTCS (black dashed) differ sig- nificantly from the classical result. Also shown is the asymptotic limit of the quantum MTCS (pink dashed) given by Eq. 2.29 with the difference from the numerical quantum MTCS decreasing as energy increases. . . . . . . . . . 24 xii Figure 2.8 Scattering angle θ vs. impact parameter b corresponding to Z = 1, κB = 1 and E = 10 a.u. for the purely classical case (black dashed curve) and with Heisenberg MDP (blue curve: αH = 1,εH = 0.125; green curve: αH = 1,εH = 0.25; red curve: αH = 1,εH = 0.5). Scattering at low impact pa- rameters is highly influenced by the MDP resulting in a structure in the scat- tering angles with the angle decreasing to zero and then increasing as impact parameter increases. Scattering angles of all the curves overlap with each other for large impact parameters. . . . . . . . . . . . . . . . . . . . . . . . . . 25 Figure 2.9 (a) Scattering angle vs. impact parameter obtained using MDP with αH = 1, εH = 0.125 for Z = 1, κB = 1 and E = 10 a.u. There is a structure in the scattering angles with the angles becoming nearly zero for an impact param- eter of b = 0.028 a.u. (b) Electron trajectory for b = 0.028 a.u. (red curve) reveals that though there is a strong interaction between the electron and the screened ion (black dot), the interaction is such that the scattering angle is nearly zero in the asymptotic limit of the trajectory. Electron trajectories for b = 0.026 a.u. (blue curve) and b = 0.03 a.u. (green curve) indicated by the vertical lines in (a) show that though the corresponding scattering angles are similar in magnitude, the nature of the trajectories are different with the electron scattered upward for b = 0.026 a.u. while the electron is scattered downward for b = 0.03 a.u. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Figure 2.10 (a) For Z = 1, κB = 1, and E = 10 to 100 a.u. the classical MTCS (blue curve) and the quantum MTCS (green curve) are compared with the MTCS using MDP (red dots) optimized with respect to its free parameters αH and εH to minimize the squared difference between the quantum MTCS and the MTCS using MDP. (b) Filled contour of MTCS using MDP for a range of αH and εH for Z = 1, κB = 1, and E = 25 a.u. The lowest value on the MTCS contour is 0.016 a.u. which is still large compared to the corresponding quantum MTCS value of 0.0122 a.u. Though the MDP significantly influences the electron- ion scattering for small impact parameters, its MTCS predictions are close to the classical MTCS suggesting the MDP’s limitation for modeling electron- ion scattering in dense plasmas. . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Figure 3.1 The logarithm (base 10) of the number of particles in a cutoff sphere with radius determined by rc = 5λs is shown. Because charged particle systems occur over very large ranges of temperature and density, the use of a simple neighbor-list algorithm may not be optimal. . . . . . . . . . . . . . . . . . . . 35 xiii Figure 3.2 The liquid portion of the Γ−κ phase diagram for a Yukawa system is shown. Three regions - caged, modulated decay and exponential decay - are shown according to the properties of the velocity autocorrelation function [200]. Five different types of charged particle systems are shown: liquid metals [147], ultracold neutral plasmas [113], dusty plasmas [145], warm dense matter [82], and conditions associated with heated capsules at the National Ignition Facility (NIF) [54]. This diagram shows the various types of behav- ior these disparate systems exhibit and their span of the Γ− κ plane. . . . . . . 36 Figure 3.3 The RMS force error is shown for the real-space (top row) and the Fourier- space (bottom row) for three values of κ = 0.1,1,2 (the columns). For dif- ferent values of α, the computed RMS errors (dots) agree well with the es- timates (curves) for the real-space part of the Ewald decomposition for (a) κ = 0.1, (b) κ = 1, and (c) κ = 2; and for the Fourier-space part of the Ewald decomposition with B-splines of order p = 2 to 7 for (d) κ = 0.1, (e) κ = 1, and (f) κ = 2. The computed RMS errors correspond to a system of 5× 103 unit charges placed at random positions in a cube. . . . . . . . . . . . . . . . . 39 Figure 3.4 Comparison of approximate RMS error estimates (curves) and full RMS er- ror estimates (dots) of the Fourier-space forces for (a) κ = 0.1 and (b) κ = 1.0 for a system of 103 unit charges. There is good agreement between the sim- plified (approximated) error estimates and the full error estimates for κ = 0.1. For κ = 1, the simplified estimates deviate from the full estimates as α de- creases below unity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Figure 3.5 RMS force errors versus the cutoff radius rc is shown. RMS force errors (dots) computed using 105 unit charges placed at random positions, with κ = 1, show very good agreement with the corresponding error estimate (dashed curve) given by (3.43). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Figure 3.6 Ratio of the cutoff radius to half the edge length of simulation cube (cid:19) (cid:18) rD c L/2 i ) and 10−6 (e2/a2 for a fixed RMS force error over a range of particle numbers. For κ = 0.1, the ratios for an error of 10−8 (e2/a2 i ) are given by the red line and the red dashed line, respectively. For κ = 0.5, the ratios for an error of 10−8 (e2/a2 i ) and 10−6 (e2/a2 i ) are given by the blue line and the blue dashed line, respectively. Also shown is the light gray region corresponding to beyond the minimum image convention (rD c > L/2) which implies that for κ = 0.1 and an error of 10−8 (e2/a2 i ), LCL would be less efficient if the system has less than about 107 unit charges, while for κ = 0.5, LCL would be less efficient for a system with less than about 105 unit charges. . . . . . . . 41 xiv Figure 3.7 Figure 3.8 S(q,Ω) obtained from PPPM-MD for a Yukawa system with Γ = 50 and κ = 1 is shown for a range of values of q. Each panel shows S(q,Ω) for one value of q, for q = kai =0.39, 0.78, 1.17, 1.56, 1.95, 2.31, which are the first six multiples of 2π/L, where L is the edge length of the simulation cube containing 103 unit charges. The grey curves are the noisy S(q,Ω) computed using the particle positions in the MD simulation. The red curves are S(q,Ω) obtained by kernel-smoothing the noisy S(q,Ω) with a Gaussian kernel of width σ = 0.01. The blue band surrounding the smooth S(q,Ω) curves shows the standard deviation in the smooth S(q,Ω) with respect to kernel smoothing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 S(q,Ω) obtained from PPPM-MD for a Yukawa system with Γ = 50 and κ = 1 is shown for a range of values of q. Each panel shows S(q,Ω) for one value of q, for q = kai =0.39, 0.78, 1.17, 1.56, 1.95, 2.31, which are the first six multiples of 2π/L, where L is the edge length of a simulation cube containing 103 unit charges. The grey curves are the noisy S(q,Ω) computed using the particle positions in the PPPM MD simulation. The red curves are S(q,Ω) obtained by kernel-smoothing the noisy S(q,Ω) with a Gaussian kernel of width σ = 0.005. The blue band surrounding the smooth S(q,Ω) curves shows the standard deviation in the smooth S(q,Ω) with respect to kernel smoothing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 Figure 3.9 Theoretical dispersion relations are compared with dispersion relation from dynamic structure factor of MD. In the left panel (a), dispersion relation from S(q,Ω) of MD smoothed with Gaussian kernels of widths σ = 0.01 (black dots) and σ = 0.005 (green dots). The green band corresponds to the full- width-half-maximum of smooth S(q,Ω) obtained using σ = 0.005 (the band for the other case with σ = 0.01 is similar to the band for σ = 0.005, hence, is not shown). Also shown are Ω(q) obtained using Eq. (3.60) with S(q) from RPA (solid red), with S(q) ≈ SHNC(0) (dashed red) and with full SHNC(q) (solid blue) where SHNC(q) refers to S(q) computed using the hypernetted chain (HNC) equations with a bridge function. In the right panel (b), full- width-half-maximum of smooth S(q,Ω) denoted by ∆/ωp with σ = 0.005 (green dots) and σ = 0.01 (blue dots). Both (a) and (b) correspond to the MD results described in Figs. 3.7 and 3.8. . . . . . . . . . . . . . . . . . . . . 50 xv Figure 3.10 Dynamic structure factor and dispersion relation for Γ = 2 and κ = 0.1. In the left panel (a), S(q,Ω) for the first six integral multiples of ∆q = 2π/L = 0.18 by averaging over data from 20 different MD simulations using 104 unit charges in a cube and force calculation using PPPM corresponding to an RMS force error of 10−6 (e2/a2 i ). For this error, LCL would be less efficient since the minimum image convention is not satisfied. The 20 different MD simulations differed in their initial conditions for particle positions and ve- locities. Each MD run was 800 ω−1 long with a time step ∆t = 0.05 ω−1 p p and frequency resolution of ∆ω (cid:39) 0.008 ωp. S(q,Ω) curves were obtained by smoothing S(q,Ω) from MD using Gaussian kernel smoothing with ker- nel width σ = 0.01. In the right upper panel (b), dispersion relation from peaks of kernel smoothed S(q,Ω) (red dots), QLCA (blue dots), Murillo’s MNSE formulation (green dots) and VE-DDFT formulation (orange dots). In the right lower panel (c), for a wider range of q, we can see significant deviation of the theoretical dispersion relations from the MD dispersion rela- tion for q greater than about 0.4. Lines connecting the dots in (b) and (c) are for guiding the eye. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 Figure 4.1 Impact of precorrelation. In the left four panels, we show the impact of precorrelating a neutral gas on the DIH process. The gray curves are the MD-predicted g(r) for a dense xenon gas (E6 potential) at four different pressures. The colored curves show g(r) for a plasma (Yukawa potential) following DIH; note that the g(r) curves at high pressures are not substan- tially different from those at low pressures prior to ionization. For the four different pressures, the right panel shows the impact of precorrelation on the coupling parameter, using several definitions of the parameter (dotted line, Γii; dashed line, Γs ii ); the colors of the curves in the right panel correspond to the colors of the g(r) curves in the left four panels. ii ; and solid line, Γeff ii; dashdot line, Γis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 (blue), 62 atm (green), and 143 atm (red). The corresponding(cid:112) Figure 4.2 Microfield distributions. P(F) for (cid:104)Z(cid:105) = 1 and initial pressures of 50 atm (cid:104)F2(cid:105) are 6.6 i , respectively. All the P(F) were converged i , and 0.1 e2/a2 i , 0.2 e2/a2 e2/a2 using 105 particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Figure 4.3 Plasma properties versus gas pressure. A wide range of physical prop- erties, including degeneracy (θ), screening (κ), species couplings (Γαγ), Coulomb logarithms (Λαγ) and equations of state from MD simulations (P) and using a mean-field approximation ((cid:101)P), are shown as a function of the ini- tial gas pressure for a xenon plasma. Note that by varying the gas pressure, many quantities of interest can be varied substantially. For all of the cases shown, the electron temperature Te ∼ 6 eV, corresponding to (cid:104)Z(cid:105) = 2; these physical properties can be controlled further by varying Te, as well. . . . . . . . 67 xvi Figure B.1 Xe ionization. The left panel shows the mean ionization (cid:104)Z(cid:105) as a function of electron temperature for different Xe densities. The right panel shows the SP prediction of the IPD energies for the first (dashed line) and the second (solid line) ionization levels of Xe as a function of electron temperature. The colors in the right panel correspond to the colors in the left panel that differentiate the Xe densities for different initial pressures. . . . . . . . . . . . . . . . . . . 87 Figure B.2 Partial radial distribution functions versus pressure. The functions gii(r) (blue), gei(r) (green), and gee(r) (red) were computed with MD using quantum- statistical potentials for 3×104 total particles (104 ions and (cid:104)Z(cid:105)=2). Note the extreme behavior of gee(r), which exceeds unity near r = 0, revealing substantial ion-mediated effective interactions as electrons localize near (rel- atively) isolated ions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 xvii KEY TO ABBREVIATIONS ICF OCP PIC VFP Inertial Confinement Fusion One-Component Plasma Particle-in-Cell Vlasov-Fokker-Planck BBGKY Bogoliubov-Born-Green-Kirkwood-Yvon MD QSP Molecular Dynamics Quantum Statistical Potential WPMD Wave Packet Molecular Dynamics MDP XRTS CP BOMD TDSE TDHF TDDFT PDE KW BFGS HF N Ne Al Ar Ca Xe Momentum-Dependent Potential X-ray Thomson Scattering Car-Parinello Born-Oppenheimer Molecular Dynamics Time-Dependent Schrödinger Equation Time-Dependent Hartree-Fock Time-Dependent Density Functional Theory Partial Differential Equation Kirschbaum-Wilets Broyden-Fletcher-Goldfarb-Shanno Hartree-Fock Nitrogen Neon Aluminum Argon Calcium Xenon xviii NIST MTCS FFT PME SPME PPPM DRF LCL TF EGS PP PM RMS MNSE HNC RPA QLCA National Institute of Standards and Technology Momentum Transfer Cross Section Fast Fourier Transform Particle-Mesh Ewald Smoothed Particle-Mesh Ewald Particle-Particle Particle-Mesh Dielectric Response Function Linked-Cell List Thomas-Fermi Exact Gradient-Corrected Screening Particle-Particle Particle-Mesh Root-Mean-Square Modified Navier-Stokes Equation Hypernetted Chain Random-Phase Approximation Quasilocalized Charge Approximation VE-DDFT Viscoelastic-Dynamic Density Functional Theory UCNP WDM DIH IPD EOS XFEL PIPP RDF E6 Ultracold Neutral Plasma Warm Dense Matter Disorder Induced Heating Ionization Potential Depression Equation of State X-ray Free-Electron Lasers Pressure-Induced Precorrelation Plasma Radial Distribution Function Exponential-6 xix AS DH SP EK CL Ashcroft-Stroud Debye-Hückel Stewart-Pyatt Ecker-Kröll Coulomb Logarithm xx CHAPTER 1 INTRODUCTION Strongly correlated charged particle systems are statistical systems of charged particles with high degree of correlation between the particles. The focus here is on a system of strongly correlated mobile charges that constitute a strongly coupled plasma [105]. The charges within the plasma interact though electromagnetic forces resulting in a number interesting phenomena. Strongly coupled plasmas occur in a wide range of scenarios such as the extremely cold and dilute ultracold plasmas [113], astrophysical plasmas like the interior of giant planets [82] and white dwarfs [23], various stages of inertial confinement fusion (ICF) [131] and dusty plasmas occurring in space and created in laboratory [207]. A simple parameter that quantifies the extent of coupling in a plasma is the Coulomb coupling parameter. The coupling parameter is defined as the ratio of average potential energy to the average kinetic energy. For simplicity, let’s consider a system of single species of charged particles embedded in a uniform background of neutralizing charges. This system is called a one-component plasma (OCP) [105] and serves as an idealized model system for real plasmas; although some plasmas in nature do satisfy the conditions for such idealization. For a system of charged particles obeying statistical physics, the kinetic energy per particle is approximately T , where T is temperature in energy units. The potential energy per particle of an OCP with number density n and electric charge Q is given by Q2/a where a is the characteristic inter-particle separation given by a = (3/(4πn))1/3; a is referred to as the Wigner-Seitz radius. The Coulomb coupling parameter is then given by Γ = Q2 aT . (1.1) When Γ exceeds unity, the plasma is characterized as a strongly coupled plasma. For values of Γ below unity, the plasma is considered to be weakly coupled. Figure 1.1 shows different kinds of plasmas in the temperate-density space with Γ contours marking regions of strongly coupled plasmas. 1 Figure 1.1 Temperature-density plot showing different kinds of plasmas (For reference see: http://plasma.szfki.kfki.hu/ zoli/LFO_webpage/LFO_2011_EN.html). Contours of Γ = 0.1,1,10 indicate regions of strongly coupled and weakly coupled plasmas. Modeling strongly coupled plasmas is significant from a fundamental science perspective of deepening our understanding of the statistical physics of correlated charged particle systems and from a technological perspective of better understanding and improving large-scale experiments like ICF experiments at the National Ignition Facility [144] and experiments on the Z-machine at Sandia National laboratory [45] to name a few. Some of the popular modeling tools for studying 3- dimensional plasmas in phase space include kinetic methods like particle-in-cell (PIC) and Vlasov- Fokker-Planck (VFP). However, these avenues have their limitations for modeling strongly coupled plasmas. A brief overview of the methods and their limitations are discussed below. PIC method allows the statistical representation of general distribution functions in phase space. PIC is used to model plasmas that have a sufficiently low collision frequency [213], that is, PIC evolves the Vlasov equation that captures the mean-field collective effects [24] in addition to 2 weak collisions modeled using Monte Carlo collision operator [213]. PIC employs fundamental equations that have the fully nonlinear effects. Space charge and other collective effects can be included self-consistently by coupling the equation evolving the phase space distribution function and the electromagnetic field equations via the source terms. Further, PIC is also a computational efficient method with implementations scaling well on massively parallel computers. PIC is widely employed for a variety of plasmas that fall within its domain including wave launchers to labora- tory plasmas with external circuits providing complex boundary conditions. For a comprehensive review of PIC see Ref. [213]. Despite the advantages offered by PIC, its inherent limitation of capturing only weak collisions renders it as an inadequate method for modeling strongly coupled plasmas. VFP system of equations have contributed to the understanding of hot plasmas generated by interaction of intense lasers with solid matter [206] such as occurring in ICF experiments. VFP provides a fully kinetic description that is particularly important when nonthermal populations play a dominant role in the system dynamics. Scenarios that are modeled well by VFP include nonlo- cal thermal electron transport, magnetic field generation and dynamics and fast electron transport for fast ignition. VFP is nothing but a linear combination of two equations: Vlasov equation that captures the mean-field collective effects characteristic of a plasma and the Fokker-Planck equa- tion that describes a combination of advection and diffusion in velocity space and is nothing but a truncated Taylor series expansion of the statistically averaged small-angle collisions. For a com- prehensive review of VFP see Ref. [206]. Despite the advantages offered by VFP for modeling ICF and similar experiments, it still suffers from the limitation of only capturing small-angle collisions. Hence, VFP is not an adequate method for modeling strongly coupled plasmas. Because the average potential energy exceeds the average kinetic energy in a strongly coupled plasma, the interaction between the charged particles plays a dominant role in determining the spatial and temporal properties of the plasma. The particle interactions, also referred to as colli- sions, are intrinsically complex due to the many-body nature and medium- to long-range nature of the interactions. Long- and medium-range interactions are a result of the Coulomb and screened- 3 Coulomb interactions, respectively, where screening is a consequence of the spatial structure of the charged particles [200]. Many-body nature of the interactions is a result of the strong corre- lation between the particles in the system. The complex collisions resulting from the many-body interactions cannot be adequately captured in the modeling methods of PIC and VFP mentioned above. When the dynamics resulting from the many-body interactions is described in phase space, it leads to a hierarchy of equations called the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [130]. Evolving the BBGKY hierarchy of equations involves solving for a 6N dimen- sional distribution function which is impossible to solve even on the fastest supercomputers of this time. This leads to the use of molecular dynamics (MD) tools for modeling strongly coupled plasmas owing to MD’s capability to solve the BBGKY hierarchy of equations with a natural dis- cretization of phase space using real particles. MD is a very popular method used in many fields beyond strongly coupled plasmas. There are a number of articles and books providing an excellent review of MD; for details of the different components of MD and its implementation see Refs. [84, 175, 212]. Although MD provides significant advantages over kinetic methods for modeling strongly correlated systems, it is limited by the length- and time-scales that are currently achiev- able [192, 71]. These scales fall well short of the experimental scales of interest, hence, there is a need for developing computational methods that can enhance the length- and time-scales of MD by orders of magnitude. Exciting experiments such as ultracold neutral plasmas, warm dense matter, x-ray free elec- tron laser experiments, nonequilibrium x-ray Thomson scattering, experiments on charged particle stopping in dense plasmas, and multispecies mixing under extreme conditions as in ICF experi- ments have further motivated the need for high fidelity MD simulations that can reach large length- and time-scales of experimental relevance. Such a capability is crucial for understanding the ex- periments to enable predictive capabilities. My thesis focuses on three areas of my research work aimed at developing such computational capabilities and to computationally identify an alternative paradigm that enables creation of small- scale controllable strongly coupled plasmas in the laboratory. This pursuit is summarized below 4 Figure 1.2 Timeline summarizing some of the essential developments in modeling non-equilibrium quantum electron dynamics at large length- and time-scales. Figure 1.3 Timeline summarizing some of the essential developments in force and potential energy calculation methods for MD simulations involving Coulomb and screened Coulomb interactions. One class of screened Coulomb interactions is the Yukawa interaction. 5 Methods for non-equilibrium quantum electrons at scale:Quantum statistical potentialsWave-packet molecular dynamicsMomentum-dependent potentialsUhlenbeck & Gropper (1932): Bose and Fermi gas Kelbg (1963): Coulomb interaction Hansen & McDonald (1981): Molecular dynamics simulation of dense hydrogen plasma 1930196020002020Dunn & Broyles (1967): Electron gas Lado (1967): density effects in Fermi and Bose gas 1980Glosli et al. (2008): Molecular dynamics simulation of temperature equilibration in dense hydrogen Time-dependent density functional theoryRunge - Gross theorem (1984)Kerman & Koonin (1976): Hamiltonian formulation of time-dependent variational principle Klakow et al. (1994): Coulomb interaction (strongly coupled systems) Grabowski et al. (2013): Wave packet spreading and localization Baczewski et al. (2015): Electronic dynamic structure factor of warm dense matter Wilets et al. (1977): Heavy ions collisions with Pauli principle; Kirshbaum & Wilets (1980): Atomic ground statesBeck et al. (2006): Protonium formationZhou et al. (2012): Sequential double ionizationEwald decomposition based methods and molecular dynamics of charged particle systems:Ewald decomposition for periodic boundary conditions (1921) (Fermi-Pasta-Ulam-Tsingou (1955): First form of molecular dynamics with 1D chain of particles with nonlinear interactions) Perram et al.(1988): O(N3/2) Ewald1920195020201980Hansen et al. (1974): Molecular dynamics for classical one-component plasmaHockney & Eastwood’s PPPM (1981): O(N log N) due to use of FFTRosenfeld (1996): Ewald for Yukawa interaction20001970Salin & Caillol (2000): Ewald for Yukawa interaction with further details of the work provided in the subsequent chapters. Some of the above listed experiments require models to resolve nonequilibrium and correlated quantum electron dynamics over large length- and time-scales involving ∼106-109 particle simu- lations over ∼1µs. An existing approach that can, in principle, be employed is the time-dependent density functional theory. But, its O(N3) scaling is prohibitive for the large length- and time-scales of interest even with a many-core multi-CPU implementation. Methods, that tend to scale well, include quantum statistical potentials (QSPs) and wave packet molecular dynamics (WPMD). But, each have their own limitations. QSPs are temperature-dependent restricting them to equilibrium dynamics and WPMD has the undesirable property of wave packet spreading and splitting. An alternative approach is a classical Hamiltonian treatment using effective quantum potentials that are momentum-dependent. Fig. 1.2 shows a timeline to summarize some of the essential devel- opments in modeling non-equilibrium quantum electron dynamics to date. Momentum-dependent potentials (MDPs) were the focus of one aspect of my research work. This work led to a publica- tion: G. Dharuman et al., Phys. Rev. E 94, 043205, (2016); the results and discussions in this publication constitute Chapter 2. We now move to Chapter 3. For some processes that occur over time-scales where electrons have equilibrated, it is more efficient (due to the large ion-electron mass ratio) to explicitly model only the correlated ions with effective screened interactions that implicitly incorporate the effects of electrons. MD is a highly effective tool in studying dynamical (equilibrium or non-equilibrium) correlated systems. The critical step in MD is the many-body force and potential energy calcula- tion. If the particles interact in a pairwise manner, the force or potential energy calculation has a naive O(N2) scaling. Fig. 1.3 provides a timeline that summarizes some of the essential develop- ments in methods for force and potential energy calculation in MD simulations involving Coulomb and screened Coulomb interactions. For many of the experimental scenarios listed above, a sys- tem of ions with screened Coulomb interactions serve as an adequate model. Since the need is for efficient simulations at scale, another aspect of my research work focused on identifying an efficient force calculation method for MD simulations with screened Coulomb interactions that 6 are medium-range. This work led to a publication: G. Dharuman et al., J. Chem. Phys. 146, 024112, (2017); the results and discussions in this publication constitute Chapter 3. We now move to Chapter 4. Having developed efficient computational methods for MD simu- lations with screened Coulomb interactions, the third aspect of my research work involved identify- ing a new paradigm (by computational means) for creation of laboratory plasmas with controllable coupling strength. With strongly coupled plasmas created in laboratory, it is not possible to control their physical state and simultaneously diagnose the same plasma using any existing experimental platform. For ex. warm dense matter is generated at very high pressures and yield transient states that are difficult to probe without large-scale facilities. Dusty plasmas tend to form monolayers and are difficult to create in moderately coupled regime. Ultracold neutral plasmas require extremely low temperatures and are difficult to produce in the strongly coupled regime. Among the exisit- ing avenues, ultracold neutral plasmas seemed to be promising, but are plagued by the process of disorder-induced heating which prevents strong coupling despite plasma creation from ultracold gases (∼ 10−6K). Since the need is for a novel avenue for creation of plasmas with controllable coupling strength in a small-scale experimental facililty with easy diagnosis, another aspect of my research focussed on identifying this new paradigm through computational means. The paradigm involves plasma creation from compressed gases at room temperature. The results and discussions of this work constitute Chapter 4. 7 CHAPTER 2 ATOMIC BOUND STATE AND SCATTERING PROPERTIES OF EFFECTIVE MOMENTUM-DEPENDENT POTENTIALS Large-scale simulations are needed to model non-equilibrium electronic dynamics in a wide vari- ety of scenarios, including stopping power experiments in dense plasmas [83, 228], multi-species mixing under extreme conditions [11, 12], non-equilibrium x-ray Thomson scattering (XRTS) [38, 120], laser-matter experiments (e.g., core ionization in x-ray free electron laser experiments [173, 215]), and ultracold neutral plasmas [113]. The need for such modeling stems from the emergence of recent large-scale experimental facilities, such as the Z machine [117, 45], Na- tional Ignition Facility [144, 122], Linac Coherent Light Source [152, 28], Deutsches Elektronen- Synchrotron [208, 35], to name a few. Further, diagnostic capabilities such as imaging XRTS [89, 88] will provide unprecedented information about the dynamical evolution of electronic states in these experiments. Coupled with recent advances in computational power that allows molec- ular dynamics simulations to span unprecedented length (multi-trillion particles) and time scales (pico to micro-seconds)[134, 116, 192, 71, 133, 163], a detailed knowledge of the non-equilibrium dynamics of such systems is, in principle, obtainable; however, it is currently not possible to per- form such large scale simulations for electronic dynamics because of the computational overhead in modeling quantum systems. Most computational approaches to electronic structure fall into three broad categories. Histor- ically, the Car-Parinello (CP) [36] method provided an avenue for coupling an electronic struc- ture calculation to ion dynamics, albeit with a fictitious electron dynamics. Similarly, Born- Oppenheimer Molecular Dynamics (BOMD) [123, 219], a limiting case of the CP method for massless electrons, forces the electronic evolution to track the (potentially non-equilibrium) ion dy- namical scales. BOMD can also be approximately extended to some electronic dynamical quanti- ties, such as the AC electrical conductivity in the Kubo-Greenwood formulation [224, 40], through the use of the Kohn-Sham orbitals and energy eigenvalues. In all three cases, the true electronic 8 dynamics is not modeled. Conversely, more direct approaches to dynamical evolution employ Time-Dependent Hartree- Fock (TDHF) [67, 136] or Time-Dependent Density Functional Theory (TDDFT) [181, 29]. TDDFT has been employed for calculating excitation spectra of atoms and molecules [135, 76] and the dy- namic structure factor of warm dense matter [8]. Conventionally, for a N particle system, TDDFT has an unfavorable O(N3) scaling that results in a few seconds per propagation step on a multi-core implementation (∼ 8000 cores) for a system of a few thousand atoms [3]; thus, TDDFT is currently quite limited to small-scale systems over short times. Moreover, incorporating finite temperature states in TDDFT remains a challenge despite recent progress in this area [142, 169]. The complete dynamics of the non-equilibrium electrons is described by a 6N-dimensional partial differential equation (PDE) (the complex, time-dependent Schrödinger equation (TDSE) in 3 spatial dimensions). Simpler alternative approaches that balance physics fidelity with lower computational cost have also been proposed. The general idea is based on mapping the quantum problem onto a framework computable in terms of a classical approach. Mapping to classical- like dynamics offers the advantage of using classical MD techniques with O(N) or O(NlogN) scaling [101] that enable large-scale simulations of interest [110, 94]. There are several avenues for constructing a classical framework for solving the time-dependent quantum problem. For ex- ample, Remacle and Levine [176] construct a classical-like framework based on ordinary differ- ential equations for the occupancies and phases. Similarly, the Gaussian-based time-dependent variational principle [191, 93] yields classical-like equations of motion. Alternatively, Schiff and Poirier [188] build an effective Lagrangian method that contains higher-order derivatives, which in turn yields classical-looking equations with extra degrees of freedom [166]. Quantum Statis- tical Potentials (QSPs) [126, 92, 109] and empirical potentials for molecular systems [189] are purely classical in their form, with effective potentials; many of these methods have been reviewed elsewhere [75]. However, the WPMD method has several undesirable properties [93], whereas the QSP method suffers from a reliance on statistical properties (e.g., temperature [109]) not well suited for describing non-equilibrium phenomena. 9 Here, we wish to replace the original TDSE with a smaller computational problem using 6N ordinary differential equations. We will employ a Hamiltonian formulation that retains the classi- cal phase space variables, but introduces a momentum-dependent potential (MDP) that contains a non-separable term to account for quantum commutator and Pauli properties. The MDP method represents the full problem in terms of a well chosen model and empirical parameters. A “well chosen” model is one that satisfies as many constraints as possible; here, the Hamiltonian formula- tion was specifically chosen because of its natural classical limit and its conservation properties (as discussed in Chapter 2). The empirical parameters to this model must be used to train the model to match a finite (and usually small) set of known properties, preferably from accurate experimental data. For these properties, the MDP model can be considered to be exact. Unfortunately, very few exact dynamical properties for quantum systems are known, and limited training of the MDP pa- rameters is possible. Given that set of parameters, the most important issue is then transferability: do the parameters chosen to match some “exactly known” property also describe those properties for which we have no prior knowledge? In fact, this strategy is similar to other approaches, such as the wavepacket approach of [204] and the machine learning approach of [197]. MDPs have been quite successful in atomic [227, 115, 48], molecular [47] and nuclear physics [220, 174, 25, 26]; however, little work has been done for bulk (plasma-like) systems [127]. In such finite-temperature electronic systems electrons undergo excitation, deexcitation, ionization and recombination. Therefore, the MDP approach for plasma-like systems needs to be established because the identity of an electron in a non-equilibrium system varies from (1) being in the ground state, (2) being in an excited state, (3) ionizing into the continuum and (4) performing free-free scattering important to transport processes. Our goal here is to establish the efficacy of the MDP approach for modeling large-scale plasma-like system through a careful examination of these four properties. For simplicity, we utilize the best known MDP, the Kirschbaum-Wilets (KW) MDP [217, 17, 52, 53, 227], which has proven very successful for bound states. A general formulation of the Hamiltonian approach for arbitrary pair MDPs is presented in Section 2.1 along with the examination of four basic MDP properties, from ground state energies to excited state properties 10 to ionization energies and, finally, free electron scattering properties. Optimization of ground state energies is discussed in Section 2.2, and transferrability of the parameters to atomic systems not in the training set is tested. Next, we turn to the examination of excited state properties in Section 2.3. The transferability of parameters among the properties is examined by using optimized ground state properties to predict first and second ionization energies; this is discussed in Section 2.4. Free-electron properties are then examined in terms of electron-ion scattering in Section 2.5. 2.1 Formulation In this section we present the Hamiltonian formulation for a non-separable MDP of the form V (rrr, ppp) that is otherwise arbitrary, including a discussion of the implied constants of motion that serve as constraints. We assume that the equations of motion derived from this Hamiltonian retain their familiar classical form. Consider an effective Hamiltonian for a system of Ne electrons and Ni ions of the form where HC is the purely classical contribution, given by N ∑ i< j ppp2 i 2mi N ∑ i=1 HC = + and HQ incorporates quantum corrections through interactions of the form (cid:104) (cid:16) rrri −rrr j, pppi − ppp j (cid:17) HQ = N ∑ i< j V H i j H = HC + HQ , (2.1) (2.2) (2.3) (cid:17)(cid:105) , , ZiZ je2 |rrri −rrr j| (cid:16) rrri −rrr j, pppi − ppp j + δsis j V P i j where N = Ne + Ni is the total number of particles, i and j are particle indices referring to particles of electron (e) or ion (I) subsystems. Ze = −1 for electron and ZI is the nuclear charge. V H is a general Heisenberg MDP between all particles and V P is a general Pauli MDP between identical particles (selected by δsis j in Eq. 2.3 where si and s j are spins of particles i and j, respectively). V P prevents two identical particles from occupying the same regions of phase space. It is important to note that the V H and V P terms do not correspond to purely kinetic or potential energies, consistent 11 with the usual commutator properties of quantum operators. The form for HQ is not arbitrary, but should be chosen for stability reasons (mitigating the Coulomb catastrophe) and to have empirical flexibility that allows its parameters to be tuned to experimental values. All forms in use have these two properties. The Hamilton equations for particle i are given by ∂ V H i j (cid:104) (cid:35) (cid:17) (cid:16) rrri −rrr j, pppi − ppp j (cid:104) ∂ V H i j + δsis jV P i j ∂ pppi (cid:16) rrri −rrr j, pppi − ppp j (cid:17) N ∑ j(cid:54)=i + (cid:34) ZiZ je2 |rrri −rrr j| − N ∑ j(cid:54)=i (cid:16) rrri −rrr j, pppi − ppp j (cid:17)(cid:105) , (cid:16) rrri −rrr j, pppi − ppp j (2.4) (cid:17)(cid:105) . + δsis jV P i j ∂rrri drrri dt = ∂ H ∂ pppi = pppi mi dpppi dt =− ∂ H ∂rrri =− ∂ ∂rrri N ∑ j(cid:54)=i (2.5) We would like to point out that ˙rrri (cid:54)= pppi; that is, the velocity is not proportional to the canonical momentum because of the addition of the non-separable MDP. An important note is that this framework naturally captures finite temperature aspects through the initial conditions for the phase space coordinates; therefore, this formulation does not suffer from the same issues as TDDFT [142], which neglects natural thermal fluctuations [8]. 2.1.1 Constants of motion Due to the non-separable terms in the Hamilton equations (Eqs. 2.4 and 2.5) as a result of the MDPs, it becomes necessary to check if the fundamental constants of motion like total energy and total angular momentum are conserved. The equation of motion of any function of phase space coordinates A(rrr1, ppp1,rrr2, ppp2, ...,rrrN, pppN) for a N particle system is given by dA dt = {A,H} + ∂ A ∂t , (2.6) where H is the Hamiltonian of the system as given by Eq. 2.3 and {A,B} denotes the Poisson bracket defined as (cid:19) {A,B} = ∂ B ∂ pppi − ∂ A ∂ pppi ∂ B ∂rrri . (2.7) (cid:18) ∂ A ∂rrri N ∑ i=1 12 dt = ∂ H When the function A is the Hamiltonian itself, then the equation of motion is given by dH ∂t , since {H,H} = 0. If the Hamiltonian is time-independent, ∂ H dt = 0, that is, a time-independent Hamiltonian even with MDPs is a constant of motion. Therefore, the Hamil- ∂t = 0 resulting in dH tonian with MDPs acts as a conserved energy, making it a useful theoretical concept and also an important tool in numerical implementations. Now, let’s consider function A to be the total angular momentum of the system given by LLLT = N ∑ i=1 LLLi = N ∑ i=1 ∑ ν Liν ˆννν , (2.8) where Liν is the νth component of particle i’s angular momentum LLLi = rrri × pppi with ν = (x,y,z). Since the total angular momentum is a time-independent quantity, its equation of motion is given by dLLLT dt = {LLLT ,H} = ∑ ν { N ∑ i=1 Liν ,H} ˆννν . (2.9) If the potential is spherically symmetric in position and momentum space, expanding the Poisson bracket of the x-component of the total angular momentum and Hamiltonian gives N ∑ { i=1 Lix,H} = = + − + pi ziyi ri N ∑ (cid:18) pyipzi i=1{Lix,H} , N ∂ H ∑ (cid:20) pzi(pyi − py j) ∂ pi i=1 N ∑ (cid:20) pyi(pzi − pz j) i=1 N ∑ i=1 ∑ j(cid:54)=i ∑ j(cid:54)=i pi j pi j (cid:19) ∂ H ∂ ri − ∂ H ∂ pi j ∂ H ∂ pi j − + ∂ H ziyi ∂ ri − ri zi(yi − y j) ri j yi(zi − z j) ri j ∂ H ∂ pi pyipzi (cid:21) (cid:21) pi ∂ H ∂ ri j ∂ H ∂ ri j . There are some obvious cancellations due to terms of equal and opposite sign resulting in N ∑ i=1 { Lix,H} = + N ∑ i=1 N ∑ i=1 ∑ j(cid:54)=i ∑ j(cid:54)=i 1 pi j 1 ri j ∂ H ∂ pi j ∂ H ∂ ri j (pz j pyi) + ∑ j(cid:54)=i 1 ri j 1 pi j ∂ H ∂ ri j N ∑ i=1 ∑ j(cid:54)=i ∂ H ∂ pi j (−pzipy j) (−ziy j) . (z jyi) + N ∑ i=1 13 (2.10) (2.11) (2.12) (cid:34) 1−(rp/ε)4(cid:105)(cid:35) (cid:104) V (r, p) = ε2 4αr2 eα Figure 2.1 Plots of r2 (in atomic units) where α = 1 and ε = 1 (a) and α = 1 and ε = 2 (b) show that for rp (cid:28) ε, the potential becomes very repulsive. Since V H and V P have similar functional forms, common notations are used to denote the potentials and their corresponding variables: V (r, p) denotes V H or V P, α denotes αH or αP, ε denotes εH or εP, r denotes ri or ri j and p denotes pi or pi j. Therefore, the highly repulsive behavior of V (r, p) as rp (cid:28) ε enforces the Heisenberg and Pauli principles within the classical framework. Since pi j = p ji, 1st and 2nd terms cancel each other. Similarly, ri j = r ji results in 3rd and 4th terms canceling each other. Therefore, Similar steps lead to {∑N i=1 Liy,H} = 0 and {∑N N ∑ i=1 { Lix,H} = 0 . i=1 Liz,H} = 0, resulting in dLLLT dt = {LLLT ,H} = 0 . (2.13) (2.14) Therefore, the total angular momentum is also a conserved quantity for a Hamiltonian with MDPs that are spherically symmetric in position and momentum space. 14 0.00.51.01.5r(a.u.)0.00.51.01.5p(a.u.)(a)r2V(r,p)0.000.070.140.210.280.350.420.490.560.630.00.51.01.5r(a.u.)0.00.51.01.5p(a.u.)(b)r2V(r,p)0.000.280.560.841.121.401.681.962.242.52 Figure 2.2 Projected surface of ∆HF (%) = 100|(E − EHF )/EHF| in (εH ,εP) plane for Al (a) and Ca (b) showing the region of minima. (εH,εP) points (white dots) corresponding to ∆HF (cid:46) 4% follow a curve indicating that there is a correlation between εH and εP values that result in ground state energies which are in very good agreement with HF values. The correlation between the parameters reveal some degree of transferability to ground state energies of other elements that are not tested. 2.2 Ground state energies In the next four sections we will examine ground state, excited state, ionization and scattering properties of MDPs, beginning in this section with ground state energies. This requires the choice of a specific MDP and we have chosen the KW MDP because ground states of many-electron atoms have been modeled with KW MDP [46] (therefore, we drop the KW designation in what follows). Before we continue it is important to specify more precisely the strategy. The empirical parameters in this model are used to train the model to match Hartree-Fock (HF) ground state energies of a subset of the elements 1. The trained model is then tested for its transferability to ionization energies and ground state energies of other elements that were not used in the training. 1This step of the training could have used experimental values; because HF is very accurate for this property, no substantive difference would have been observed. 15 1.01.11.21.31.41.51.6εH0.81.01.21.41.61.82.02.2εP(a)∆HF(%)0.06.513.019.526.032.539.045.552.058.51.01.11.21.31.41.51.6εH0.81.01.21.41.61.82.02.2εP(b)∆HF(%)0.09.619.228.838.448.057.667.276.886.4 Figure 2.3 (a) Relative error ∆HF (%) of MDP prediction of ground state energies with respect to HF values (green circles connected by dashed green line) for N, Ne, Al, Ar and Ca show that the KW formulation is an excellent model for ground state energies (lines are to guide the eye). (b) Correlated (εH,εP) points (circles in shades of red from light to dark corresponding to increasing atomic number) extracted from the ∆HF surfaces. Also shown are curve fits (solid curves in shades of red increasing from light to dark corresponding to increasing atomic number) using εP = Aε2 H + BεH +C for N, Ne, Al, Ar and Ca. For all the elements considered, the curve fits match well with the points extracted from the corresponding ∆HF surfaces suggesting a possible transferability with respect to the atomic number. The MDP interaction between an electron and the nucleus is stabilized by the Heisenberg MDP V H(ri, pi) = H ¯h2 ε2 4αHmer2 i eαH (cid:104) 1−(ri pi/εH ¯h)4(cid:105) , (2.15) where me is mass of electron, ri and pi are the magnitudes of position and momentum of ith electron relative to the nucleus and εH and αH are parameters. Because V H is more repulsive, scaling as r−2, than the attractive electron-ion attraction, scaling as −r−1, the electron resides in a potential at a finite distance from the nucleus. The ion is taken to be at rest relative to the electron due to the large ion-electron mass ratio. However, the KW prescription only includes a Heisenberg interaction between species and no statistical interaction between ions by its definition 2. For a many-electron atom, the Pauli exclusion principle is incorporated through the Pauli MDP 2In principle, Heisenberg and statistics arguments apply between all species, including the ions. 16 1.051.151.251.351.45εH0.751.001.251.501.752.002.25εP(b)NNeAlArCa6810121416182022Z0.0050.0100.0150.0200.025∆HF(%)(a) Figure 2.4 (a) Coefficients A, B and C of the fit εP = Aε2 H + BεH + C vary as a function of atomic number (Z). The pattern in the data points of A (green circles), B (black star) and C (red squares) corresponding to N, Ne, Al, Ar and Ca are captured well by 4th degree polynomial fits for A (solid green curve), B (dotted black curve) and C (dashed red curve). (b) Correlated (εH,εP) curves for F (blue dashed) and S (blue solid) computed with A, B and C interpolated using their corresponding 4th degree polynomial fits. Ground state energies for F and S from their correlated (εH,εP) are in good agreement with HF values (with a maximum ∆HF of ∼10% and ∼5% for F and S, respectively). This confirms the transferability of the trained (εH,εP) to ground state energies of elements that were not included in the training. expressed as V P(ri j, pi j) = P¯h2 ε2 4αPmer2 i j eαP 1−(ri j pi j/εP¯h)4(cid:105) (cid:104) , (2.16) where ri j and pi j are the magnitudes of relative position and momentum of ith and jth electron respectively and εP and αP are parameters. V H tries to impose the condition ripi ≥ εH ¯h which is analogous to the Heisenberg uncertainty principle and V P tries to impose the condition ri j pi j ≥ εP¯h which is analogous to the Pauli exclusion principle. The region excluded from phase space by Pauli and Heisenberg MDPs is referred to as the ’core’. εH and εP are parameters that decide the size of the core while αH and αP are parameters that decide the hardness of exclusion from the core. Note that V H and V P have similar functional forms. Therefore, they are illustrated in a combined Here, we deliberately use the KW MDP in its original form to retain a strong connection with the previous literature on MDPs. The development of new MDPs would naturally include a more general form. 17 1.051.151.251.351.45εH0.751.001.251.501.752.002.25εP(b)FS68101214161820Z−200−150−100−500A,-B,C(a)AA:polynomialfit-B-B:polynomialfitCC:polynomialfit manner as r2V (r, p) in Fig. 2.1 where V denotes V H or V P, r denotes ri or ri j and p denotes pi or pi j. From Fig. 2.1 we can infer that for ripi (cid:28) εH ¯h and ri j pi j (cid:28) εP¯h the potentials become very repulsive, thereby enforcing the Heisenberg and Pauli principles within the classical framework. The properties of the ground state depend on the values of the core sizes and the hardness parameters. For simplicity 3, we try to quantify the influence of core sizes on the ground state energies; therefore, we fix the values of the hardness parameters αH and αP to be 2 and 1, respec- tively, as suggested by Beck et al. [18] based on their stopping power studies. The ground state of a many-electron atom is obtained by minimization of the Hamiltonian with respect to positions and momenta of the electrons keeping their spins fixed which requires a si- multaneous solution for the set of equations ∂ H ∂rrri = 0 for i = 1 to Ne. The minimized = 0 and ∂ H ∂ pppi Hamiltonian would result in a frozen configuration for the ground state with zero electron veloc- ities but non-zero momenta, as mentioned below Eq. 2.5. Thus, the MDP model has the desired non-classical behavior. The energy of the minimized Hamiltonian gives the ground state energy E. The electrons are assigned successive spin values of 1/2 and −1/2 and are initialized with certain position and momentum values prior to minimization. Atomic units were used for the calcula- tion. Following [46], minimization was performed using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [81] as implemented in the MATLAB solver ‘fminunc’. BFGS is an uncon- strained optimization method belonging to the class of quasi-Newton methods. Despite the use of BFGS we cannot be sure that the optimization results in a global minimum. To identify the core sizes that resulted in ground state energies in close agreement with HF (Hartree-Fock) values we performed the minimization for a range of core sizes for some elements. For an atom with M elec- trons, the Hamiltonian is a function of 6M variables, therefore, the cost of minimization increases as M becomes larger. On a laptop with Intel Core i3 processor (2.53 GHz and 4 GB RAM) it took about 4 seconds and 5 minutes to obtain the ground states of nitrogen and calcium respectively. 3Although all the four parameters could have been varied for the training, we restricted our training to only the core sizes for simplicity and to identify any underlying pattern in the core sizes that could help towards transferability. 18 A search in the parameter space requires a number of minimizations, therefore, we performed the search only for some elements, namely, nitrogen (N), neon (Ne), aluminum (Al), argon (Ar) and calcium (Ca). 2.2.1 Correlated core sizes (εH, εP) The percentage deviation from HF defined as ∆HF (αH ,αP) = 100|(E − EHF )/EHF| was traced as a surface for N, Ne, Al, Ar and Ca for εH spanning from about 1 to 2 and εP spanning from about 1 to 2.5. The space was discretized with a grid spacing of 0.01, therefore about 15000 minimizations were performed for each element considered. After parallelizing, minimizations to trace the surface for Ca on a node with about 10 cores took about 2 days. The search was crucial because the search led us to identify correlated εH and εP values that resulted in ground state energies in very good agreement with HF values for all the elements considered. Fig. 2.2 shows the projected surface of ∆HF for Al and Ca and the correlated (εH, εP) points that correspond to ∆HF (cid:46) 4%. Using the same condition, the correlated set of (εH, εP) were extracted from ∆HF surfaces of other elements considered (N, Ne and Ar) as shown in Fig. 2.3(b) (points). Relative error with respect to HF values ∆HF (%), are shown in Fig. 2.3(a). For each of these elements the ground state energy with minimum ∆HF are comparable with ground state energy experimentally measured [121]. We observed that the correlated (εH, εP) follow a pattern with respect to atomic number (Z), that is, the correlation extends to higher values of εH with increasing atomic number, as seen in Fig. 2.3(b). The reason for this pattern is not understood. The correlation between εH and εP is captured well by a parabolic relation expressed as εP = A(Z)ε2 H + B(Z)εH +C(Z). (2.17) There is no particular reason for expressing εP as a function of εH; a similar fit would yield εH as a function of εP. We quantified the pattern with respect to atomic number by fitting the coefficients A, B and C to fourth degree polynomials as shown in Fig. 2.4(a). We then interpolated 19 Figure 2.5 (a) Electron trajectory around ion (black dot) corresponding to different initial condi- tions: Kepler-like motion from t = 0 to t = 50 a.u. (blue curve), zig-zag motion from t = 0 to t = 50 a.u. (dark green curve), zig-zag motion from t = 50 to t = 231 a.u. (light green curve) and reciprocating motion from t = 0 to t = 50 a.u. (red curve). Reciprocating motion corresponds to zero initial angular momentum which would result in an unstable trajectory in the absence of Heisenberg MDP. (b) Magnitude of position (r) and magnitude of momentum (p) corresponding to the trajectories in (a) from t = 0 to t = 50 a.u. are superimposed on the Hamiltonian surface. (c) Region 1 of (b) magnified to show the band-like structure in (r, p) dynamics. (d) Region 2 of (b) magnified to show a similar band-like structure in (r, p) dynamics. 20 1.51.00.50.00.51.01.5x(a.u.)1.51.00.50.00.51.01.52.0y(a.u.)(a)Coulomb-liketrajectory:t=0to50(a.u.)Zig-zagtrajectory:t=0to50(a.u.)Zig-zagtrajectory:t=50to231(a.u.)Reciprocatingtrajectory:t=0to50(a.u.)Ion0.80.91.01.11.21.3r(a.u.)0.80.91.01.11.21.3p(a.u.)(b)H(a.u.)(r,p):Zig-zagtrajectory(r,p):Reciprocatingtrajectory(r,p):Coulomb-liketrajectory0.60.50.40.30.20.10.00.10.20.3210.00010.00020.00030.00040.00050.0006r(a.u.)+9.998⇥1010.000100.000150.000200.00025p(a.u.)+9.998⇥101(c)H(a.u.)(r,p):Zig-zagtrajectory(r,p):Reciprocatingtrajectory(r,p):Coulomb-liketrajectory0.60.50.40.30.20.10.00.10.20.30.8450.8500.8550.860r(a.u.)1.1241.1261.1281.1301.1321.134p(a.u.)(d)H(a.u.)(r,p):Zig-zagtrajectory(r,p):Reciprocatingtrajectory0.60.50.40.30.20.10.00.10.20.3(c)(d)H (a.u.)H (a.u.)H (a.u.) Figure 2.6 For N, F, Ne, Al, S, Ar and Ca, first and second ionization energies were computed using MDPs with the correlated (εH, εP) optimized to give accurate neutral ground-state energies. MDP prediction of first ionization energies (red curve) are in good agreement with NIST data (dashed blue curve); second ionization energies using MDPs (green curve) yield mixed results compared to NIST data (dashed cyan curve). Therefore, the parameters trained on neutral ground state energies transfer fairly well to the prediction of first and second ionization energies with some outliers. the coefficients for fluorine (F) and sulphur (S) to obtain their correlated (εH, εP) that resulted in ground state energies in good agreement with HF. The maximum ∆HF was ∼10% and ∼5% for F and S respectively. This is remarkable because it confirms the transferability of the trained (εH,εP) to ground state energies of F and S that were not included in the training. Fig. 2.4(b) shows the interpolated (εH, εP) curves of F and S along with the parabolic fits of (εH, εP) for N, Ne, Al, Ar and Ca. 21 0510152025Z0.00.51.01.52.02.53.0IonizationEnergy(a.u.)NISTdata:FirstionizationenergyKWMDP:FirstionizationenergyNISTdata:SecondionizationenergyKWMDP:SecondionizationenergyNFNeAlSArCaNFNeAlSArCa 2.3 Excited state orbits In this section we turn to excited state properties, defined as ˙rrr (cid:54)= 0, ˙ppp (cid:54)= 0, and H < 0. Again, we make use of the MDP, but for an electron-ion pair with the simpler Hamilton equations given by ∂Veff(r, p) , ppp me drrr dt = dppp dt = − + ∂ ppp ∂Veff(r, p) ∂rrr , (2.18) where Veff(r, p) contains both the attractive Coulomb potential and the (repulsive) Heisenberg MDP. In contrast with the minimization procedure used in the previous section, we now examine several initial value problems for (2.18). The Hamilton equations (2.18) with ionic charge Z = 1, αH = 5, and εH = 0.9535 were numerically integrated (after conversion to atomic units) using MATLAB’s RK45 integrator. We considered the following initial conditions: (i) rrr = 1ˆxxx, ppp = 1ˆyyy, (ii) rrr = 1ˆxxx, ppp = 1√2 As shown in Fig. 2.5(a), the trajectory is Kepler-like (blue curve) for initial condition (i), is of ˆyyy, (iii) rrr = 1ˆxxx, ppp = 1ˆxxx; these choices yield quite different behaviors. ˆxxx + 1√2 zig-zag nature (green curve) for (ii) and is a reciprocating pattern (red curve) for (iii). Note that condition (iii) corresponds to zero initial angular momentum, which is an important test of any MDP; zero initial angular momentum would result in a Coulomb catastrophe in the absence of the proper MDP. The stabilizing nature is due to the dominant 1/r2 contribution in the MDP in the limit of p = 0 or rp = εH, compared to the (infinitely deep) attractive Coulomb potential −1/r. We found the energy and angular momentum to be conserved for all three trajectories as expected from Section 2.1. One way of interpreting the trajectories’ nature is by superimposing their (r, p) dynamics on the total energy contour as shown in Fig. 2.5(b). Depending on the initial condition, the electron executes a trajectory that confines its energy to a contour marked by its initial energy. This implies that the minimum of the confining potential Veff(r, p) changes with the electron’s momentum in accordance with energy conservation. Further, though Fig. 2.5(b) gives the impression that (r, p) dynamics for the three different trajectories overlap with each other, 22 Figs. 2.5(c) and 2.5(d) showing the magnified regions 1 and 2 of Fig. 2.5(b) reveal a kind of band structure in their (r, p) dynamics that do not exactly overlap with each other. 2.4 First and second ionization energies The identity (being bound or free) of an electron in a finite-temperature plasma is continuously changing. Moreover, in processes such as charged-particle stopping, much of the energy loss can be due to ionization [138]. Therefore, in addition to atomic properties, MDPs must accurately capture transitions between bound and free states. In this subsection we examine these properties through comparisons of predicted first and second ionization energies with experimental values. We do this in a manner that allows us to assess the transferable properties of the MDPs by using parameters previously trained on ground-state properties. The nth ionization energy is given by In = En − E0 , (2.19) where En is the energy of an ion with n electrons removed and E0 is the (ground state) energy of the neutral atom. En for n =1, 2 were obtained by minimizing the corresponding Hamiltonian using the same minimization algorithm employed above for every correlated εH and εP pair trained on the ground state energy of the corresponding neutral atom. The corresponding set of first and second ionization energies were then computed. From this set, we chose those that minimized the combined error defined as ∆EGS,1,2 = (EGS − EHF)2 + (I1 − I1,expt.)2 + (I2 − I2,expt.)2 , (2.20) where I1,expt. and I2,expt. are experimentally measured first and second ionization energies given by National Institute of Standards and Technology (NIST) [121]. As shown in Fig. 2.6, first ionization energies are in good agreement with the NIST data, while the second ionization energies yield mixed results. Therefore, the parameters εH and εP trained to give very accurate neutral ground 23 Figure 2.7 Comparison of the classical and quantum MTCS for Z = 1, κB = 1, and E = 10 to 100 a.u. The classical MTCS using semi-analytic method (blue dashed) and the trajectory method (cyan dashed) are in very good agreement. They also match with the asymptotic limit of the classical MTCS (red dashed) given by Eq. 2.28. The quantum MTCS (black dashed) differ significantly from the classical result. Also shown is the asymptotic limit of the quantum MTCS (pink dashed) given by Eq. 2.29 with the difference from the numerical quantum MTCS decreasing as energy increases. state energies transfer well to the prediction of first and second ionization energies, with some outliers. 2.5 Scattering properties of free electron states In this section we examine our fourth criteria for plasma behavior: scattering properties of contin- uum states. We quantify the MDP’s ability to accurately describe continuum properties through tests based on the momentum transfer cross section (MTCS), an important quantity related to stop- ping power [6, 160, 172] and other transport properties [180]. As in previous sections, we reduce the many-body Hamiltonian to a simpler system of an electron scattered by a screened ion. The 24 101102E(a.u.)10−310−210−1σtr(a.u.)ClassicalσtrusingtrajectorymethodClassicalσtrusingsemi-analyticmethodClassicalσtr:asymptoticlimitQuantumσtrusingWKBphaseshiftsQuantumσtr:asymptoticlimit Figure 2.8 Scattering angle θ vs. impact parameter b corresponding to Z = 1, κB = 1 and E = 10 a.u. for the purely classical case (black dashed curve) and with Heisenberg MDP (blue curve: αH = 1,εH = 0.125; green curve: αH = 1,εH = 0.25; red curve: αH = 1,εH = 0.5). Scattering at low impact parameters is highly influenced by the MDP resulting in a structure in the scattering angles with the angle decreasing to zero and then increasing as impact parameter increases. Scattering angles of all the curves overlap with each other for large impact parameters. Hamiltonian in the reference frame of the ion is given by HC = p2 2me Ze2 VY (r) = − r +VY (r) , e−r/λ (2.21) (2.22) where Z is the ionic charge and λ is the screening length, which is chosen appropriate to a plasma; for example, for dense plasmas, λ is typically a finite temperature Thomas-Fermi screening length [199]. For the rest of this section atomic units have been used and λ is expressed through an inverse screening parameter κB = aB/λ where aB is the Bohr radius. We begin by comparing the quantum and classical MTCS to identify the conditions where they differ appreciably. The quantum MTCS is given by σ QM tr = 4π k2 ∞ ∑ l=0 (l + 1) sin2(δl − δl+1) , (2.23) where k is the magnitude of the wavevector of a free electron and δl is the phase shift for angular momentum quantum number l. For simplicity, we chose the WKB approximation for δl [6, 180] 25 0.00.10.20.30.40.5b(a.u.)0.00.51.01.52.02.53.03.5θ(rad.)αH=1,H=0.125αH=1,H=0.25αH=1,H=0.5PurelyClassical Figure 2.9 (a) Scattering angle vs. impact parameter obtained using MDP with αH = 1, εH = 0.125 for Z = 1, κB = 1 and E = 10 a.u. There is a structure in the scattering angles with the angles becoming nearly zero for an impact parameter of b = 0.028 a.u. (b) Electron trajectory for b = 0.028 a.u. (red curve) reveals that though there is a strong interaction between the electron and the screened ion (black dot), the interaction is such that the scattering angle is nearly zero in the asymptotic limit of the trajectory. Electron trajectories for b = 0.026 a.u. (blue curve) and b = 0.03 a.u. (green curve) indicated by the vertical lines in (a) show that though the corresponding scattering angles are similar in magnitude, the nature of the trajectories are different with the electron scattered upward for b = 0.026 a.u. while the electron is scattered downward for b = 0.03 a.u. which is given by (cid:115) k2 − The classical MTCS [180] is given by δl = (cid:90) dr (cid:90) (cid:115) k2 − dr (l + 1/2)2 r2 . (l + 1/2)2 r2 − 2VY (r)− (cid:90) ∞ 0 db (1− cosθ (b))b , σCL tr = (2.24) (2.25) where θ (b) is the scattering angle for an impact parameter b. For a Hamiltonian without an MDP the scattering angle is given by the semi-analytic formula (cid:113) 1− b2/r2 −VY (r)/E where rm is the classical turning point given by the largest root of the equation θ (b) = π − 2b dr r2 rm (cid:90) ∞ 1 , b2 r2m − 1− VY (rm) E = 0 . 26 (2.26) (2.27) 0.000.010.020.030.040.050.06b(a.u.)0.00.51.01.52.02.53.03.5✓(rad.)(a)↵H=1,✏H=0.1250.030.020.010.000.010.020.03z(a.u.)0.010.000.010.020.030.040.05x(a.u.)(b)b=0.026b=0.028b=0.03ScreenedIon Figure 2.10 (a) For Z = 1, κB = 1, and E = 10 to 100 a.u. the classical MTCS (blue curve) and the quantum MTCS (green curve) are compared with the MTCS using MDP (red dots) optimized with respect to its free parameters αH and εH to minimize the squared difference between the quantum MTCS and the MTCS using MDP. (b) Filled contour of MTCS using MDP for a range of αH and εH for Z = 1, κB = 1, and E = 25 a.u. The lowest value on the MTCS contour is 0.016 a.u. which is still large compared to the corresponding quantum MTCS value of 0.0122 a.u. Though the MDP significantly influences the electron-ion scattering for small impact parameters, its MTCS predictions are close to the classical MTCS suggesting the MDP’s limitation for modeling electron- ion scattering in dense plasmas. Now, for a Hamiltonian with an MDP, Eq. 2.26 cannot be applied; therefore, we computed the scattering angle from the electron’s trajectory which was obtained by numerically integrating its Hamilton equations until the electron-ion interaction became negligible. We refer to this as the trajectory method for MTCS. For the purely classical case with Z = 1, κB = 1 and energy range E = 10 to 100 a.u., we compared the MTCS from trajectories with the semi-analytic MTCS (Eq. 2.26). They are in good agreement as shown in Fig. 2.7, thereby validating our implementation. For the same conditions, the quantum MTCS differ appreciably from the classical MTCS as shown in Fig. 2.7. Also shown in Fig. 2.7 are asymptotic limits of the classical and quantum MTCS values (denoted by σCM tr,asym. respectively) derived in [59] tr,asym. and σQM that have analytic expressions given by σCM tr,asym. = 4π (cid:18) Z 2E (cid:19)2 (cid:20) ln (cid:19) (cid:18)4E Z λ − γ − 1 2 (cid:21) , (2.28) 27 101102E(a.u.)103102101tr(a.u.)(a)ClassicaltrusingtrajectorymethodQuantumtrusingWKBphaseshiftstrusingKWMDPoptimizedw.r.t.(✏H,↵H)0.050.100.150.200.25✏H0.10.20.30.40.50.60.70.80.91.0↵H(b)tr(a.u.)0.0160.0180.0200.0220.0240.0260.0280.030 (cid:18) Z 2E (cid:19)2 (cid:20) ln (cid:16) 2√2E λ (cid:17) (cid:21) , 1 2 − σQM tr,asym. = 4π (2.29) where γ = 0.577 [59]. These expressions reveal that in the asymptotic limit, the quantum MTCS qualitatively differs from the classical MTCS due to a lower limit set by the deBroglie wavelength on the distance of closest approach for quantum scattering. Therefore, the test for MDP is if it can include the necessary quantum effects to bridge the gap between classical and quantum MTCS values. Since the interaction is between an electron and a screened ion, Heisenberg MDP is added to HC (Eq. 2.21). Using the trajectory method, we computed scattering angles for a range of impact parameters for αH = 1 and εH = 0.125, 0.25, and 0.5, as shown in Fig. 2.8. We found a structure in the scattering angles for low impact parameters which is due to the MDP’s influence on the electron-ion interaction as evident in the electron trajectories for b = 0.026, 0.028, and 0.03, as shown in Fig. 2.9(b). Those trajectories correspond to the impact parameter range where the scattering angle decreases to zero and then increases (as marked by the vertical dotted lines in Fig. 2.9(a)). We performed an optimization with respect to the free parameters, αH and εH, to minimize the squared difference between the quantum MTCS and the MTCS using MDP. The optimized MTCS using MDP are close to the purely classical values, as shown in Fig. 2.10(a), implying that the MDP doesn’t incorporate the required quantum effects despite having a strong effect on the electron-ion scattering for small impact parameters. Fig. 2.10(b) shows the filled contour of MTCS for a range of αH and εH for Z = 1, κB = 1 and E = 25 a.u. The lowest MTCS value on the MTCS contour is about 0.016 a.u. which is larger than the corresponding quantum MTCS value of 0.0122 a.u. Similar observations were made from MTCS contours for energies in the range of 10 to 100 a.u. for Z = 1 and κB = 1. Therefore, though MDP serves as a good model for ground state energies and first ionization energies of many-electron atoms, it is unable to incorporate the quantum effects in scattering of a free electron by a stationary screened ion, suggesting its limitation for plasmas. 2.6 Summary and conclusions In summary, we have examined a time-dependent, quantum-mechanical method for large-scale simulations of non-equilibrium systems as an alternative to more expensive methods such as TD-DFT, TD-HF etc. and which relaxes most limitations associated with the relatively fast WPMD [93] and QSP [92] methods. In 28 particular, our focus has been on the use of a classical-like, Hamiltonian-based framework based on effective momentum-dependent interactions [220] that has desirable conservation properties and the computational scaling of standard classical molecular dynamics. For simplicity, we employed the KW MDP form [220] since they have been quite successful for many atomic properties [17, 47, 49, 50, 53, 217, 129, 127, 143, 51, 227]. We examined their strengths and weak- nesses for use in non-equilibrium dense plasma simulations using four criteria. We first trained the parameters based on HF calculations to give accurate ground state energies for neutral atoms with atomic number less than 20; the excellent agreement with HF calculations is shown in Fig. 2.3(a). We found a correlation between the parameters and fitted the correlation to a parabolic relation as shown in Fig. 2.3(b), revealing a level of transferability of these parameters to previously untrained systems. Next, we computed the properties of excited state orbits of electron-ion pair and found disparate proper- ties of the trajectories depending on the initial angular momentum, including unusual reciprocating patterns. An important feature of the KW MDP is the stabilization of the Coulomb catastrophe for the special case of zero initial angular momentum (see reciprocating motion case in Fig. 2.5(a)). Because plasma electrons persistently undergo ionization and recobination events, we then turned to the ionization process itelf, with the ionization energy as our quality metric for the MDP. From the fixed, previously-determined ground state parameters we predicted the first- and second-ionization energies and found first ionization ionization energies to be in good agreement with experimental values, again suggesting good transferability. However, the second ionization energies were mixed, with half accurate and half with an error as large as a factor of two, as shown in Fig. 2.6. The reason for this behavior in the second ionization energy is currently unknown. Finally, we examined continuum states responsible for electronic transport processes, such as stopping power [6] and electrical and thermal conductivity [180]. We chose to examine the ability of the KW MDP to reproduce the MTCS versus energy for electron-ion collisions. A screened interaction was chosen because it more realistically represents the dense plasma environment and because the Rutherford MTCS has patho- logical properties (i.e., a large impact parameter divergence) in this context. Using a WKB approach [180], we computed the classical and quantum MTCS and compared with predictions from KW MDP trajectories. As we showed in Fig. 2.10, KW MDP cannot yield the correct MTCS despite training the parameters, sug- 29 gesting that the functional form itself is responsible. Thus, although KW MDPs were successful for ground state properties, they cannot capture scattering properties important for plasma simulations. The results presented here were for isolated atoms and ions. To quantify the implications of our findings on many-body properties of a plasma, large-scale molecular dynamics simulation are needed; this is beyond the scope of the present work and will be explored in a future work. However, our results have suggested several areas of improvement. First, the KW forms did not include a Heisenberg interaction between elec- trons, which is unphysical. Second, the functional form of the KW appears to have emerged to recover a Bohr picture of the ground state, which yields poor properties for scattering states. Third, alternate training methods should be explored, including optimizing on several properties (e.g., ground state energy and cross sections) simultaneously. Such explorations are underway and will be the subject of a future work. 30 A GENERALIZED EWALD DECOMPOSITION FOR SCREENED COULOMB INTERACTIONS CHAPTER 3 The interactions employed in MD simulations are typically considered to be either short- or long-range. In d dimensions, interactions of the form 1/rn are considered short range if n > d. Short-range interactions are typically treated with a fast neighbor search within a cutoff radius rc; the interaction is truncated thereafter, and a tail correction can be added [193, 140]. Conversely, if n ≤ d, then the interactions are considered long range; such interactions are considerably more difficult to treat. In periodic systems, most methods for handling long-range interactions are based on the Ewald sum [74, 210], which decomposes the Coulomb potential into two sums that are rapidly convergent. This decomposition can be implemented through a wide variety of algorithms, of which a naive implementation yields O(N2) scaling; a simple optimization reduces this scaling to O(N3/2) [164]. Faster, grid-based methods have been developed that allow the fast Fourier transform (FFT) to be used[101], resulting in O(N logN) scaling; among these methods are the particle-mesh Ewald (PME) [57], smoothed particle-mesh Ewald (SPME) [73] and particle-particle particle-mesh (PPPM) methods [101, 61, 201, 61, 167, 9]. Alternatively, non-Ewald methods [86] have also been developed, such as the Wolf summation method[222] and fast multipole methods [124], and hybrid methods [194, 66] have been developed, as well. However, this categorization of interactions as either short range or long range reflects the mathematical issue of force-sum convergence. In fact, many interactions can be considered to be medium range when their convergence is guaranteed yet either there are a prohibitive number of particles in the neighbor radius or there are unacceptable truncation errors. For example, power-law potentials of the form ∼ 1/r2 (charge- dipole), ∼ 1/r3 (dipole-dipole) and 1/r6 (dispersion) fall into this category [158, 60], as well as Sutton- Chen [205] potentials relevant for metallic glasses [171]. Isele-Holder et al. [106] have extended the PPPM algorithm to interactions other than 1/r in their work on 1/r6 dispersion interactions. Similarly, prior work on screened Coulomb interactions includes extending the Ewald summation to Yukawa (exponentially screened) interactions [177, 185]. Here, we examine a particular class of medium-range interactions that form with linearly screened Coulomb potentials, where the screening is specified through an arbitrary dielectric response function (DRF) ε(kkk). The Ewald sum is formulated for this screened interaction and implemented with both a direct force 31 approach (linked-cell list (LCL) [218] alone) and the PPPM method. An error analysis is carried out to examine the crossover between the direct force calculation and the PPPM implementation as a function of the range of the interaction and the particle number. In Section 3.1, we describe a class of interactions for screened charged systems obtained through the Fourier representation of the dielectric properties of background species. Several choices for the screening function are given in terms of the dielectric response function ε(kkk), and their ranges are discussed. In Section 3.2, the Ewald decomposition is developed for charged particle systems interacting through an effective interaction involving ε(kkk). Section 3.3 includes error estimates and timing studies for the Yukawa dielectric response function. Finally, we apply our results to compute the dynamic structure factor for a screened Coulomb system and discuss the implications of the choice of ε(kkk) and the computational costs that result from the temperature and density dependence of the effective interaction. 3.1 Linearly screened Coulomb systems Consider N point particles with charges {Qi} (e.g., ions), which can be partial, in the presence of a back- ground charge density −enb (e.g., electrons). The electrostatic forces in such a system can be calculated from the Poisson equation 1 4π − ∇2Φ(rrr) = N ∑ i=1 Qiδ (rrr−rrri(t))− enb(rrr), (3.1) where Φ(rrr) is the total electrostatic potential, δ (rrr) is the delta function, and {rrri(t)} are the particle co- ordinates at time t. For simplicity, in the subsequent steps, the dependence of {rrri} on time is not explic- itly shown. Once this equation is solved, the corresponding force on the ith particle is calculated to be FFFi = −Qi∇Φ(rrr = rrri). Expression (3.1) can be equivalently represented in Fourier space as k2 4π Φ(kkk) = N ∑ i=1 Qie−ikkk·rrri − enb(kkk). (3.2) By assuming that the background density evolves instantaneously with the point particles, which is a valid assumption for a very light species, such as electrons, and approximating the structure of the background density as a linear response to the potential fluctuations induced by the presence of the point particles, we 32 can express the background density in terms of a response function χ(kkk) as nb(kkk) = −v(k)χ(kkk) N ∑ i=1 Qie−ikkk·rrri, v(k) = 4π k2 . We can write the electrostatic potential as Φ(kkk) = v(k) ε(kkk) N ∑ i=1 Qie−ikkk·rrri, where the reciprocal of the DRF can be expressed in terms of the response function as Finally, the force on the ith particle can be represented in terms of the pair interactions as ε−1(kkk) = 1 + v(k)χ(kkk). where the effective interaction ui j(rrr) is given in Fourier space as FFFi = − ∑ j(cid:54)=i ∇rrriui j(rrri −rrr j), ui j(kkk) = QiQ j v(k) ε(kkk) . (3.3) (3.4) (3.5) (3.6) (3.7) When the background does not respond, χ(kkk) = 0 in (3.3). Such a completely uniform background cor- responds to a dielectric response of unity in (3.5), and inverting (3.7) recovers the usual Coulomb interaction ui j(r) = QiQ j/r. Clearly, in the limit ε → 1, Ewald-class methods are needed. To continue further, we must make a specific choice of the DRF; we begin with the most widely used form ε(k) = 1 + (λsk)−2, (3.8) where λs is the screening length associated with the background fluctuations. The associated pair interaction is given in real space by ui j(r) = QiQ j r e−r/λs. (3.9) This interaction decays exponentially and is therefore “short" range. However, as λs → 0, the number of particles in a neighbor-sphere diverges. It is therefore convenient to characterize the strength of the screening in terms of the dimensionless quantity κ = ai λs , 33 (3.10) where ai = (3/4πni)1/3 is the ion-sphere radius, which is a measure of the interparticle distance. The κ = 0 case is the pure Coulomb (one-component plasma) limit; thus, despite the exponential decay for finite κ, (3.9) has a long-range interaction as one of its limits. This suggests that there may be a finite range of κ for which the interaction (3.9) is medium range in the sense that there is no clear, optimal choice of algorithm. To provide a concrete example, we can examine this choice in the context of the well-known Thomas- Fermi (TF) model [200]1, for which λT F = (cid:118)(cid:117)(cid:117)(cid:117)(cid:116)(cid:16) (cid:17)1/2 , (3.11) T 2 + 4E2 F /9 4πnbe2 where T is the background temperature in energy units, EF is the Fermi energy of the electron background in energy units, e is the elementary charge, and nb is the background number density. Let us suppose that, for example, an acceptable error is given by a choice of cutoff radius rc = 5λs, which yields the factor e−5 ∼ 0.674· 10−2 in (3.9). Allowing physical conditions (T and nb) to vary, we can examine the number of particles in the neighbor cell. This is shown in Fig. 3.1, where contours show the logarithm of the number of particles in such a sphere as a function of temperature and density. Condensed systems (e.g., ionic liquids, liquid metals, electrolytic solutions, etc.) tend to lie near the lower right region of the graph. Laser-produced plasmas originate as cold solids, in that same region, but can be rapidly heated to the very high temperatures of the upper right region. Roughly speaking, the lower right half of the diagram corresponds to ∼ 103 particles in the neighbor cell for this choice of error. In Fig. 3.2, we show the Γ− κ phase diagram. The coupling parameter Γ is given by Γ = Z2e2 aiTi , (3.12) where Z is the degree of ionization and Ti is the ionic temperature. This diagram reflects the regimes of different charged particle systems of interest such as liquid metals [147], ultracold neutral plasmas [113], dusty plasmas [145], warm dense matter [82] and conditions at the National Ignition Facility [54]. We also show the regions corresponding to different liquid-like properties based on the properties of their velocity autocorrelation function [200]. Note that these disparate systems cover a wide range of κ values. Of course, the long-wavelength form of (3.8) is not unique, and interactions other than (3.9) are possible [199]. Moreover, the method is applicable to a system of partial charges as well. For this reason, we will 1A slightly more accurate fit for λT F is discussed in Ref. [200]. 34 Figure 3.1 The logarithm (base 10) of the number of particles in a cutoff sphere with radius de- termined by rc = 5λs is shown. Because charged particle systems occur over very large ranges of temperature and density, the use of a simple neighbor-list algorithm may not be optimal. retain the generality of an arbitrary ε(kkk) throughout the Ewald development in the next section. In fact, medium range interactions of different nature occur in biological and chemical systems where the linear screening approach is not applicable [87, 68, 170]. 3.2 Ewald decomposition for arbitrary dielectric response functions To proceed with the development of the Ewald composition of (3.4), we note that the central quantity is the potential-like quantity φ (kkk) = v(k) ε(kkk) . (3.13) As is common, we implement the Ewald decomposition by adding and subtracting a screening-charge dis- tribution of the form ζρ(rrr), where ζ is a constant to be determined, and ρ(rrr) is chosen to have the Gaussian form ρ(rrr) = α3 π3/2 e−α2r2 , (3.14) where α is the Ewald splitting parameter. Fourier-space representation of the decomposition takes the form + v(k) ε(kkk) ζ e−k2/4α2 (3.15) (3.16) (cid:20) 1− ζ e−k2/4α2(cid:21) φ (kkk) = v(k) ε(kkk) = φR (kkk) + φF (kkk), 35 1012141618202224log(n∗cc)−3−2−10123log(T/eV)-1.01.03.05.07.09.0NumberofparticlesinfiveThomas-Fermilengths(log)−1.50.01.53.04.56.07.59.010.5 Figure 3.2 The liquid portion of the Γ− κ phase diagram for a Yukawa system is shown. Three regions - caged, modulated decay and exponential decay - are shown according to the properties of the velocity autocorrelation function [200]. Five different types of charged particle systems are shown: liquid metals [147], ultracold neutral plasmas [113], dusty plasmas [145], warm dense matter [82], and conditions associated with heated capsules at the National Ignition Facility (NIF) [54]. This diagram shows the various types of behavior these disparate systems exhibit and their span of the Γ− κ plane. where the real-space contribution (cid:20) 1− ζ e−k2/4α2(cid:21) φR (kkk) = v(k) ε(kkk) (3.17) (3.18) can be highly short range (because there are effectively two sources of cutoff). Analogously, φF (kkk) = v(k) ε(kkk) ζ e−k2/4α2 is short range in Fourier-space and hence can also be computed efficiently; this is the subject of the next section. In the absence of screening, ε(k) → 1, and we recover the Ewald decomposition for the standard Coulomb interaction [210]. The real-space representation of φR (rrr) is given by the inverse Fourier transform 1 (cid:90) (cid:90) ∞ (2π)3 2 πr 0 φR (rrr) = = (cid:20) 1− ζ e−k2/4α2(cid:21) 1− ζ e−k2/4α2(cid:21) (cid:20) d3k v(k) ε(kkk) sin(kr) ε(k) dk k eikkk·rrr , (3.19) (3.20) 36 ultracold neutral plasmasNIF capsule NIF burnNIFdusty plasmaswarm dense matterexponential decaycagedliquid metalsmodulated decay where an isotropic dielectric function ε(k) is assumed in the second line. We evaluate this first with the long-wavelength (k → 0) form of ε(k) given by ε−1 Y (k) = (3.21) k2 k2s + k2 . Here, ks is the magnitude of the inverse screening vector given by ks = 1/λs, where λs is the long- wavelength screening length associated with properties of the background (i.e., classical, quantum, mul- tispecies, etc.). This form yields what is generically referred to as the Yukawa interaction. In what follows, we will be using k−1 s in place of λs. Substituting ε−1 Y (k) in (3.20) yields (cid:20) (cid:19) (cid:18) αr− ks 2α φY R (r) = 1 2r e−ksrerfc (cid:18) (cid:19)(cid:21) ks 2α + eksrerfc αr + , (3.22) with ζ ek2s /4α2 [177, 185, 137]. = 1 (see Appendix A.1). Alternative derivations of φY R (r) have been given previously Gradient corrections to the Yukawa result can be captured by the form ε−1 EGS(k) = k2(k2 s + 1 4νk2) k4s + k2(k2s + 1 4νk2) , (3.23) where leading-order corrections to the Yukawa form appear as k4 terms. An example of this type of DRF is the exact gradient-corrected screening (EGS) potential [199], in which k4 terms arise from either electronic correlations or Heisenberg uncertainty (see Appendix A.1). Interestingly, the form of the potential is qual- itatively different for ν < 1 and ν > 1, with the latter corresponding to an oscillatory decay [199]. Hence, the Ewald decomposition for the two cases will be different; in this paper, we will restrict our discussion to the ν < 1 case. Substituting ε−1 EGS(k) in (3.20) yields φ EGS R (r) = A−φ−(r)− A+φ+(r) φ±(r) = ek2 ζ 1 e−k±r + r 2r − e−k±rerfc (cid:18)k±2α + αr ±/4α2(cid:20) (cid:18)k±2α − αr (cid:19)(cid:21) ek±rerfc (cid:16) 1±√1− ν k2 ± = 2 ν , (cid:17) k2 s . (cid:19) (3.24) (3.25) (3.26) where A± = s − 1 4νk2 k2 ± √1− νk2s , 37 Because k− < k+ (for ν < 1), the optimal choice is ζ = e−k2 −/4α2 in Appendix A.1. . We have given the detailed derivation The real-space part of the Ewald decomposition of the screened Coulomb interaction between two point charges for an arbitrary DRF is, therefore, given by ui j(R)(rrri −rrr j) = QiQ j φR (rrri −rrr j). (3.27) This completes the formulation of the real-space contribution, which will be common to all Ewald-based methods. In the next section, we first turn to the Fourier-space contribution, which we will treat with the PPPM algorithm for the particle-mesh contribution. Next, we provide an error analysis and a time-per-time- step comparison of the force calculations using PPPM and LCL algorithms for the Yukawa interaction. For the remainder of this paper, when we discuss LCL for a Yukawa interaction, this should be understood to refer to a force calculation using the LCL algorithm for a plain Yukawa interaction that is not subject to Ewald decomposition. 3.3 Error analysis and timing studies In this section, we begin with a brief explanation of the PPPM algorithm, followed by an error analysis for the computed force for the Yukawa DRF, with errors due to the truncation and discretization that are inherent to the algorithm. Then, we compare the execution times of the PPPM and LCL algorithms for the Yukawa potential to distinguish between the two algorithms in terms of computational efficiency for a certain level of error in the force calculation. The results of this section and the next section were obtained using the codes we implemented 2. The particle-particle (PP) part of the PPPM algorithm refers to the real-space part of the Ewald de- composition, which can be highly short range; therefore, interactions beyond a certain cutoff radius can be neglected. If the number of particles within a cutoff sphere Nc is small compared to the total number of particles N in a simulation box, then the total number of interactions O(NNc) tends to O(N) as Nc becomes 2The codes are written in high performant Python using the optimizing compiler Numba which compiles Python code to machine code giving similar performance to C and Fortran. More details about Numba can be found here: http://numba.pydata.org. 38 Figure 3.3 The RMS force error is shown for the real-space (top row) and the Fourier-space (bottom row) for three values of κ = 0.1,1,2 (the columns). For different values of α, the computed RMS errors (dots) agree well with the estimates (curves) for the real-space part of the Ewald decomposition for (a) κ = 0.1, (b) κ = 1, and (c) κ = 2; and for the Fourier-space part of the Ewald decomposition with B-splines of order p = 2 to 7 for (d) κ = 0.1, (e) κ = 1, and (f) κ = 2. The computed RMS errors correspond to a system of 5× 103 unit charges placed at random positions in a cube. much smaller than N (Nc (cid:28) N). Such an O(N) scaling is achieved by the LCL algorithm [101], in which the simulation box is divided into a number of cells with dimensions proportional to the cutoff radius. The particle-mesh (PM) portion of the PPPM algorithm refers to the Fourier-space part of the interaction that is computed on a grid using one of the two algorithms in [[201]]; these algorithms differ in the manner in which forces are computed on a grid. Both versions of the PPPM algorithm employ the highly efficient FFT, which scales as O(M logM), where M is the number of mesh points in the Fourier-space grid. Errors are associated with computing the energy and forces in real-space and Fourier-space. Because we are interested in MD simulations in which the trajectories of a system of particles evolve primarily because of the forces acting on the particles, we restrict our error analysis to the errors in the forces. We chose to quantify the error in a force using the root-mean-square (RMS) error [201, 61]. The RMS errors in the real-space and Fourier-space components of the forces, which are acting on particles in a simulation box, 39 24681012rc(ai)10−1310−1210−1110−1010−910−810−710−610−510−4∆FYR(e2/a2i)κ=0.1(a)α=0.5α=1α=1.5α=20.51.01.52.02.5α(a−1i)10−910−810−710−610−510−410−310−210−1∆FYF(e2/a2i)κ=0.1(d)p=2p=3p=4p=5p=6p=724681012rc(ai)10−1310−1210−1110−1010−910−810−710−610−510−4κ=1(b)α=0.5α=1α=1.5α=20.51.01.52.02.5α(a−1i)10−1010−910−810−710−610−510−410−310−210−1κ=1(e)p=2p=3p=4p=5p=6p=724681012rc(ai)10−1310−1210−1110−1010−910−810−710−610−510−4κ=2(c)α=0.5α=1α=1.5α=20.51.01.52.02.5α(a−1i)10−1210−1110−1010−910−810−710−610−510−410−310−210−1κ=2(f)p=2p=3p=4p=5p=6p=7 Figure 3.4 Comparison of approximate RMS error estimates (curves) and full RMS error estimates (dots) of the Fourier-space forces for (a) κ = 0.1 and (b) κ = 1.0 for a system of 103 unit charges. There is good agreement between the simplified (approximated) error estimates and the full error estimates for κ = 0.1. For κ = 1, the simplified estimates deviate from the full estimates as α decreases below unity. Figure 3.5 RMS force errors versus the cutoff radius rc is shown. RMS force errors (dots) computed using 105 unit charges placed at random positions, with κ = 1, show very good agreement with the corresponding error estimate (dashed curve) given by (3.43). 40 0.51.01.52.0α(a−1i)10−1010−910−810−710−610−510−410−310−2∆FF(e2/a2i)(a)κ=0.1approx.,p=3full,p=3approx.,p=4full,p=4approx.,p=5full,p=5approx.,p=6full,p=6approx.,p=7full,p=70.51.01.52.0α(a−1i)10−1010−910−810−710−610−510−410−310−2(b)κ=1510152025303540rc(ai)10−1710−1510−1310−1110−910−710−510−3∆FY(e2/a2i)computedwithparticlesestimate (cid:19) (cid:18) rD c L/2 i ) and 10−6 (e2/a2 for a fixed Figure 3.6 Ratio of the cutoff radius to half the edge length of simulation cube RMS force error over a range of particle numbers. For κ = 0.1, the ratios for an error of 10−8 (e2/a2 i ) are given by the red line and the red dashed line, respectively. For κ = 0.5, the ratios for an error of 10−8 (e2/a2 i ) are given by the blue line and the blue dashed line, respectively. Also shown is the light gray region corresponding to beyond the minimum image convention (rD c > L/2) which implies that for κ = 0.1 and an error of 10−8 (e2/a2 i ), LCL would be less efficient if the system has less than about 107 unit charges, while for κ = 0.5, LCL would be less efficient for a system with less than about 105 unit charges. i ) and 10−6 (e2/a2 have a common form given by (cid:118)(cid:117)(cid:117)(cid:116)∑N i=1 ∆Fξ P = (cid:104) FFFapprox. i,P −FFFexact i,P N (cid:105)2 , (3.28) where P = R, F refers to real-space or Fourier-space, FFFapprox. i, FFFexact i,P is the approximate force acting on particle i,P is the actual force (or numerically exact force) acting on particle i, and ξ is used to distinguish be- tween computed RMS errors and analytical estimates of RMS errors. Moreover, because the computational costs associated with real-space and Fourier-space force calculations are determined by the levels of error that are allowed, it is essential to understand how the errors vary with respect to the parameters. For this purpose, we have derived analytical estimates for the errors in the real-space and Fourier-space components of the forces for a system of randomly distributed particles. The error estimates were verified by comparing with the RMS errors computed for a system of particles placed at random positions within a simulation box. For the remainder of this paper, we use ai as the unit of distance and the elementary charge e as the unit of charge; thus, the RMS force error is in units of e2/a2 i . 41 103104105106107108109N102101100101102rc/(L/2)FY=106,=0.1FY=108,=0.1FY=106,=0.5FY=108,=0.5beyond minimum imageminimum image The error in the real-space component of the force is due to the truncation of particle interactions beyond a certain cutoff radius. For a system of N charged point particles interacting pair-wise with an arbitrary short- range interaction expressed as us(rrri −rrr j) = QiQ j φs(rrri −rrr j), the estimate of the RMS force error due to a cutoff radius rc can be derived by adapting the error analysis in [118], which considers the real-space portion of the Ewald decomposition for the Coulomb potential. With the assumption that beyond rc, particles are randomly distributed, the estimate of the RMS force error due to us(rrr) is given by ∆Fs = Vc(r≤rc<∞) Q √NV i=1 Q2 i , and V is the volume enclosing the particles. R (rrr) in (3.29) and making several approximations, we get the following estimate of the RMS force error for the real-space portion of the Ewald decomposition for the Yukawa where fff s(rrr) = −∇rrrφs(rrr), Q = ∑N R (rrr) = −∇rrrφY d3r |fff s(rrr)|2 Substituting fffY (3.29) (cid:20)(cid:90) (cid:21)1/2 , potential: ∆FY R (cid:39) 2Q √NV e−κ2/4α2 e−α2r2c √rc (3.30) (see Appendix A.2 for details). Note that as κ → 0 (Coulomb limit), we recover the corresponding error estimate for the Coulomb case [118]. Figs. 3.3a-c show very good agreement between the error estimates and the computed RMS error in the real-space force component for a system of 5× 103 unit charges with κ = 0.1, 1, and 2 and varying values of α. The choice of these κ values is consistent with the κ range of interest as shown in Fig. 3.2. The particles were randomly distributed in a cube of edge length L given by L = (4πN/3)1/3 in units of ai. In the Fourier-space part of the PPPM algorithm, the force calculation for a system of charged point particles with arbitrary spatial distribution is mapped to a real-space grid by assigning particle positions to the real-space grid points using B-splines of some order p. Calculations are then performed in Fourier-space starting with the FFT of the charge density. This is followed by computing the potential energy and forces on the Fourier-space grid following the steps described in [[201]]. Forces on the Fourier-space grid are then transformed to real-space by an inverse FFT. Next, the forces are interpolated from the real-space grid to the actual particle positions using B-splines of the same order, in a manner similar to the earlier step of assigning particle positions to the real-space grid. FFT involves an inherent truncation in the maximum allowable length of Fourier-space vectors, which is one cause for error in the Fourier-space component of the forces. In addition, errors occur because of aliasing, which arises because calculations are being made 42 on a grid, and because of the interpolation of particle positions to the real-space grid and of forces computed on the real-space grid to actual particle positions. Critically, the PPPM algorithm employs an optimal Green function [101] that is derived by minimizing the error in the forces for a random distribution of particles. The optimal Green function can be interpreted as a weighted average of Fourier modes of the actual Green function. For the version of the PPPM algorithm in which forces are computed on a grid [201], the optimal Green function corresponding to an arbitrary Green function is given by [201] ˆGkn+mkn+mkn+mknknkn ·kn+mkn+mkn+m ∑mmm ˆU2 ˆG(cid:48)knknkn = (cid:20) kn+mkn+mkn+m ∑mmm ˆU2 (cid:21)2 kn+mkn+mkn+m |knknkn|2 , (3.31) where nnn refers to the grid vector denoted by the triplet of grid indices (nx,ny,nz), ˆGknknkn is the actual Green function evaluated at the Fourier-space vector knknkn corresponding to grid vector nnn, mmm refers to the triplet of grid indices (mx,my,mz) that contribute to aliasing, and ˆUknknkn is the Fourier transform of the B-spline of order p. ˆUknknkn is given by (cid:20)sin(πnx/Mx) (cid:21)p (cid:34)sin(cid:0)πny/My(cid:1) (cid:35)p (cid:20)sin(πnz/Mz) (cid:21)p , (3.32) πnx/Mx πny/My πnz/Mz ˆUknknkn = where Mx, My and Mz refer to the number of grid points along the respective directions in Fourier-space. The triplet of grid sizes (hx,hy,hz) are then given by (Lx/Mx, Ly/My, Lz/Mz). The estimate of the optimal RMS error in the Fourier-space force is given by [201] where F V 2/3 = ∑ χ2 kkk(cid:54)=0 ˆGkkk|kkk|2 −∑ nnn ∆FF (cid:39) Q √NV 2/3 χF , (cid:104) ∑mmm ˆU2 (cid:20) kn+mkn+mkn+m ∑mmm ˆU2  (cid:21)2 ˆGkn+mkn+mkn+mknknkn ·kn+mkn+mkn+m kn+mkn+mkn+m |knknkn|2 (cid:105)2  . (3.33) (3.34) Figs. 3.3d-f show very good agreement between the RMS error estimates for the Fourier-space force with order of B-splines ranging from 2 to 7 and the corresponding RMS error computed using the same system of 5× 103 unit charges that was employed for the real-space error analysis. Because the simulation volume is a cube, the numbers of grid points M along each direction in Fourier-space are equal. Eq. (3.33) is not of 43 κ 0.3 0.5 1 1.5 2 N 1×105 1×106 1×105 1×106 1×105 1×106 1×105 1×106 1×105 1×106 ) i α (a−1 0.625 0.579 0.614 0.569 0.61 0.57 0.684 0.635 0.684 0.683 rc (ai) M rD 5.2 5.6 128 256 c (ai) TR (s) TF (s) TLCL (s) TLCL/TPPPM 35.48 35.62 324.2 7138.24 71.56 136 2.61 25.65 1.92 26.81 5.3 5.7 5.33 5.8 4.6 5.2 4.5 4.7 128 256 128 256 128 256 128 256 22.24 22.33 12.08 12.34 7.99 8.38 6.4 6.47 1.97 24.79 1.36 23.16 133.6 1173.05 1.89 26.1 1.35 21.2 1.28 18.66 1.22 26.34 1.2 13.2 2.17 22.8 11.7 119.9 3.07 34.76 1.83 16.3 40.1 24.5 3.6 2.53 1.23 0.77 0.54 0.45 Table 3.1 Comparison of times per time-step for force calculations using the PPPM and LCL algorithms for Yukawa interaction, for a range of κ and the number of unit charges (N). The entries correspond to an error of approximately 10−5 (e2/a2 i ). All times are in seconds. α is the Ewald splitting parameter, rc is the cutoff radius for the real-space calculation, M is the number of grid points along each direction in Fourier-space, and rD c is the cutoff radius for LCL. For all Fourier-space calculations, B-spline of order p = 6 was found to be optimal. TR and TF are the times per time-step of real-space and Fourier-space calculations, respectively. TPPPM refers to the total time per time-step for PPPM and is given by TPPPM = TR + TF . TLCL is the time per time-step for LCL. a form that explicitly reveals the influence of the parameters h, α and p on the error incurred in the Fourier- space force calculation. Therefore, (3.33) requires further simplification, which is possible if the arbitrary Green function decays sufficiently rapidly in Fourier-space that the alias contributions can be neglected. Following the approach in [[61]], (3.33) is approximated to ∆F(approx) F Q √NV (cid:39) A1/2 F , where AF is given by (2π)2 β (p,m) denotes an integral given by AF (cid:39) 3 p−1 ∑ m=0 C(p) m β (p,m) = (cid:19)2(p+m) (cid:18)h (cid:90) ∞ 2 0 2 [1 + 2(p + m)] β (p,m) ; dk ˆG2 k k2(p+m+2) , 44 (3.35) (3.36) (3.37) m are coefficients listed in Table I of [[61]] (see Appendix A.2 for further detail about the simpli- results in a and C(p) fication of (3.33)). For the Yukawa interaction, substituting ˆGk = modified form for β (p,m) given by e−(k2+κ2)/4α2 (k2+κ2) 4π β (p,m,κ/α) = (4π)2α2(p+m)+1 ψ(p,m,κ/α) , where ψ(p,m,κ/α) denotes an integral of the form (cid:90) ∞ κ/α e−u2/2 u3 du u2 − κ2/α2(cid:17)p+m+3/2 (cid:16) . ψ(p,m,κ/α) = (3.38) (3.39) Finally, a simplified form of the RMS error estimate for the Fourier-space component of the force for the Yukawa interaction is given by ∆FY (approx) F = where (cid:18)hα 2 BY F = p−1 ∑ m=0 C(p) m 2Q √NV (cid:19)2(p+m) (3α)1/2(cid:16) BY F (cid:17)1/2 , 2 [1 + 2(p + m)] ψ(p,m,κ/α) . (3.40) (3.41) F Fig. 3.4 shows good agreement between ∆FY (approx.) and the full error estimate ∆FY 1 for a system of 103 unit charges. For κ = 1 as α decreases below unity, ∆FY (approx.) ∆FY of interest. Nevertheless, the simplified form of ∆FY (approx.) approximately varies as the product of the grid size and the Ewald splitting parameter (hα) raised to the F for the α values is useful because it reveals that the error F . Further, we found that for κ > 1, ∆FY (approx.) F for κ = 0.1 and deviates from is not a good approximation of ∆FY F F F power given by the order of B-spline (p). Further, finding an optimal set of parameters requires considerable effort using (3.33), whereas using a simplified estimator can expedite the process of finding the optimal parameters. An alternative approach for estimating the errors is described in Ref. [154]. The total RMS error in the force computed using the PPPM algorithm is defined as (cid:113) ∆Ftot = ∆F2 R + ∆F2 F . (3.42) To compare the computational costs of the PPPM and LCL algorithms, the execution times for the two algorithms must be compared for the same RMS force errors. To this end, an estimate of the RMS force 45 error for Yukawa interaction is needed. Substituting fffY (rrr) = −∇rrrφY (rrr) in (3.29), where φY (rrr) = e−κr/r, results in the following estimate of the RMS force error for Yukawa interaction: ∆FY (cid:39) This equation, in turn, implies that rD c (cid:39) − 1 κ log 2πκ Q √NV √2πκ e−κrc . (cid:33) (cid:32) (cid:114) NV ∆FY Q (3.43) (3.44) (see Appendix A.2 for further detail). As shown in Fig. 3.5, the RMS error estimate in (3.43) has been verified by comparison with the RMS error computed using 105 unit charges placed at random positions in a box. Based on the error analysis for the Yukawa interaction, Fig. 3.6 shows that for κ = 0.1 and an RMS force error of 10−8 (e2/a2 i ), the minimum image convention rc ≤ L/2 is not satisfied below a system size of 107 unit charges. Therefore, for a κ value and particle numbers for which LCL would be less efficient, PPPM serves as an excellent alternative. As κ increases for approximately the same error of 10−8 (e2/a2 i ), the threshold above which the particle number satisfies the minimum image convention decreases; for example, for κ = 0.5, the minimum image convention is satisfied beyond 105 unit charges, as shown in Fig. 3.6. Importantly, for some conditions where the LCL algorithm is applicable, we find the PPPM algorithm to be computationally advantageous. This observation is based on the comparison of the time per time step of PPPM and LCL algorithms for a fixed RMS force error. The timing results for a range of κ and N (with unit charges) are provided in Table I. For example, with 106 unit charges and an RMS force error of approximately 10−5 (e2/a2 i ), we find that the PPPM algorithm was approximately 25 and 130 times faster than the LCL algorithm for κ = 0.5 and 0.3, respectively (as shown in Table I). As κ approaches unity, we find the LCL and PPPM algorithms to be comparable in performance, as summarized by the entries in Table I for κ = 1. As κ increases beyond unity, we find that the PPPM algorithm loses its computational advantage over the LCL algorithm, as can be seen in the entries in Table I for κ = 2. We should point out, however, that execution times are architecture dependent and our conclusions are based on that fact. All calculations used a 3.5 GHz Intel Xeon(R) processor and 32 GB RAM. 46 3.4 Application: Dynamic Structure Factor In this section, we apply our results to the computation of the dynamic structure factor S(k,ω). Our interest in S(k,ω) falls into two categories. Theoretically, S(k,ω) contains all of the many-body processes associated with thermal density fluctuations, including collective modes and transport [119, 108, 187, 147]. All of these processes are obviously impacted by the choice of ε(k), and we wish to use MD for model validation. The ionic S(k,ω) is also connected directly with neutron scattering [30] and, indirectly, to X-ray Thomson scattering (XRTS) [148]. In the context of warm dense matter physics, for example, recent XRTS studies have developed a modified Navier-Stokes equation (MNSE) to model the ionic S(k,ω) [147]. Within the framework of the Yukawa model, we are interested in density and temperature conditions for which MD is required to validate approximate theoretical models and where LCL would be computationally less efficient than PPPM for force calculations. Consider the case Γ = 2 and κ = 0.1, for 104 unit charges the PPPM algorithm can have a force error of approximately 10−6 (e2/a2 i ) for an optimal choice of α, rc, h and p. For the same error and particle number, LCL would be less efficient because the minimum-image convention is not satisfied. This case corresponds to an effective coupling of Γeff = Γe−κ (cid:39) 0.74, which is of an order that is relevant to systems of recent interest, such as warm dense matter [148]. 3.4.1 Numerical computation of the dynamic structure factor from MD The dynamic structure factor for an N-particle system is defined as S(kkk,ω) ≡ 1 N dt eiωt (cid:104)n(kkk,t) n(−k−k−k,0)(cid:105). (cid:90) ∞ −∞ (cid:90) T Here, the brackets (cid:104)···(cid:105) refer to the stationary time average (cid:104)n(kkk,t)n(−k−k−k,0)(cid:105) = 1 T dτ n(kkk,t + τ)n(−k−k−k,τ), 0 (3.45) (3.46) where T is the length of the run. We note that for a system in equilibrium, the stationary time average is equivalent to the ensemble average as T → ∞. The Fourier-transformed densities n(kkk,t) are obtained from the definition of the density, n(rrr,t) = δ (rrr−rrr j(t)) , (3.47) N ∑ j=1 47 Figure 3.7 S(q,Ω) obtained from PPPM-MD for a Yukawa system with Γ = 50 and κ = 1 is shown for a range of values of q. Each panel shows S(q,Ω) for one value of q, for q = kai =0.39, 0.78, 1.17, 1.56, 1.95, 2.31, which are the first six multiples of 2π/L, where L is the edge length of the simulation cube containing 103 unit charges. The grey curves are the noisy S(q,Ω) computed using the particle positions in the MD simulation. The red curves are S(q,Ω) obtained by kernel- smoothing the noisy S(q,Ω) with a Gaussian kernel of width σ = 0.01. The blue band surrounding the smooth S(q,Ω) curves shows the standard deviation in the smooth S(q,Ω) with respect to kernel smoothing. to yield (cid:90) n(kkk,t) = Substituting (3.46) into (3.45), we obtain S(kkk,ω) = 1 NT (cid:20)(cid:90) T d3r n(rrr,t) e−ikkk·rrr = (cid:90) ∞ (cid:90) T −∞ dτ n(−k−k−k,τ) 0 e−ikkk·rrr j(t) . N ∑ j=1 dt eiωt n(kkk,t + τ) , (3.48) (3.49) which, through the change of variables y = t + τ, can be written as dτ n(−k−k−k,τ) e−iωτ We use the usual definition of the temporal Fourier transform S(kkk,ω) = 0 1 NT −∞ (cid:21)(cid:20)(cid:90) ∞ (cid:21) dy eiωy n(kkk,y) . (3.50) n(kkk,ω) = dy eiωy n(kkk,y) . (3.51) (cid:90) ∞ −∞ For numerical purposes, the infinite domain must be approximated as finite, and as we are calculating sta- tionary processes, we can take this domain to be [0,T ] without loss of generality for sufficiently large T . 48 0.00.20.40.60.81.01.21.40.00.51.01.52.02.5S(q,Ω)q=kai=0.39σ=0.010.00.20.40.60.81.01.21.40.00.51.01.52.02.5q=kai=0.78σ=0.010.00.20.40.60.81.01.21.40.00.51.01.52.02.5q=kai=1.17σ=0.010.00.20.40.60.81.01.21.4Ω=ω/ωp0.00.51.01.52.02.5S(q,Ω)q=kai=1.56σ=0.010.00.20.40.60.81.01.21.4Ω=ω/ωp0.00.51.01.52.02.5q=kai=1.95σ=0.01rawdatasmoothS(q,Ω)witherrorband0.00.20.40.60.81.01.21.4Ω=ω/ωp0.00.51.01.52.02.5q=kai=2.34σ=0.01 Figure 3.8 S(q,Ω) obtained from PPPM-MD for a Yukawa system with Γ = 50 and κ = 1 is shown for a range of values of q. Each panel shows S(q,Ω) for one value of q, for q = kai =0.39, 0.78, 1.17, 1.56, 1.95, 2.31, which are the first six multiples of 2π/L, where L is the edge length of a simulation cube containing 103 unit charges. The grey curves are the noisy S(q,Ω) computed using the particle positions in the PPPM MD simulation. The red curves are S(q,Ω) obtained by kernel-smoothing the noisy S(q,Ω) with a Gaussian kernel of width σ = 0.005. The blue band surrounding the smooth S(q,Ω) curves shows the standard deviation in the smooth S(q,Ω) with respect to kernel smoothing. This yields the approximate transform n(−kkk,−ω) = (cid:90) T 0 dτ n(−k−k−k,τ) e−iωτ . Therefore, the final form for S(kkk,ω) is given by (3.52) (3.53) (3.54) as n(−kkk,−ω) = n∗(kkk,ω). The Fourier-transformed density n(kkk,ω) is given by S(kkk,ω) = 1 NT |n(kkk,ω)|2 , (cid:90) T (cid:105) (cid:104) kkk·rrr j(t)−ωt dt e−i (cid:104) kkk·rrr j(nt∆t)−ω∆t e−i . (cid:105) n(kkk,ω) = N ∑ j=1 0 n(kkk,ω) = N ∑ j=1 Nt−1 ∑ nt =0 Particle positions from MD are available only at certain times in [0,T ]; hence, n(kkk,ω) is computed using , (3.55) where Nt is the number of steps after the equilibration phase of MD, and ∆t is the step size for time inte- gration. Particle positions are stored at intervals corresponding to this step size. Long runs are preferred 49 0.00.20.40.60.81.01.21.40.00.51.01.52.02.5S(q,Ω)q=kai=0.39σ=0.0050.00.20.40.60.81.01.21.40.00.51.01.52.02.5q=kai=0.78σ=0.0050.00.20.40.60.81.01.21.40.00.51.01.52.02.5q=kai=1.17σ=0.0050.00.20.40.60.81.01.21.4Ω=ω/ωp0.00.51.01.52.02.5S(q,Ω)q=kai=1.56σ=0.0050.00.20.40.60.81.01.21.4Ω=ω/ωp0.00.51.01.52.02.5q=kai=1.95σ=0.005rawdatasmoothS(q,Ω)witherrorband0.00.20.40.60.81.01.21.4Ω=ω/ωp0.00.51.01.52.02.5q=kai=2.34σ=0.005 Figure 3.9 Theoretical dispersion relations are compared with dispersion relation from dynamic structure factor of MD. In the left panel (a), dispersion relation from S(q,Ω) of MD smoothed with Gaussian kernels of widths σ = 0.01 (black dots) and σ = 0.005 (green dots). The green band corresponds to the full-width-half-maximum of smooth S(q,Ω) obtained using σ = 0.005 (the band for the other case with σ = 0.01 is similar to the band for σ = 0.005, hence, is not shown). Also shown are Ω(q) obtained using Eq. (3.60) with S(q) from RPA (solid red), with S(q) ≈ SHNC(0) (dashed red) and with full SHNC(q) (solid blue) where SHNC(q) refers to S(q) computed using the hypernetted chain (HNC) equations with a bridge function. In the right panel (b), full-width-half-maximum of smooth S(q,Ω) denoted by ∆/ωp with σ = 0.005 (green dots) and σ = 0.01 (blue dots). Both (a) and (b) correspond to the MD results described in Figs. 3.7 and 3.8. because the noise in S(kkk,ω) decreases as T increases. The temporal Fourier transform is computed by FFT. S(kkk,ω) computed from MD data tends to be noisy, and the extent of the noise increases as k increases. To determine the location of the peak of S(kkk,ω) and its width, the noisy curves are usually smoothed by aver- aging over data from different MD runs that differ in their initial conditions. In the following section, we present a potentially more efficient smoothing approach that employs kernel smoothing [196, 85]. 3.4.2 Kernel smoothing Kernel smoothing refers to weighted averaging of noisy observations to obtain a smooth function [196, 85]. To smooth our data using kernel smoothing, we use a kernel regression technique [85] defined by ˜y( ˜x) = ∑B i=1 Kσ ( ˜x,xi) y(xi) ∑B i=1 Kσ ( ˜x,xi) , (3.56) 50 0.00.51.01.52.02.5q=kai0.00.20.40.60.81.01.2Ω=ω/ωp(a)Γ=50,κ=1ΩwithSRPA(q)ΩwithS(q)≈SHNC(0)ΩwithSHNC(q)PPPM-MD(σ=0.01)PPPM-MD(σ=0.005)0.00.51.01.52.02.5q=kai0.000.050.100.150.200.250.300.350.400.45∆/ωp(b)σ=0.01σ=0.005 Figure 3.10 Dynamic structure factor and dispersion relation for Γ = 2 and κ = 0.1. In the left panel (a), S(q,Ω) for the first six integral multiples of ∆q = 2π/L = 0.18 by averaging over data from 20 different MD simulations using 104 unit charges in a cube and force calculation using PPPM corresponding to an RMS force error of 10−6 (e2/a2 i ). For this error, LCL would be less efficient since the minimum image convention is not satisfied. The 20 different MD simulations differed in their initial conditions for particle positions and velocities. Each MD run was 800 ω−1 and frequency resolution of ∆ω (cid:39) 0.008 ωp. S(q,Ω) p curves were obtained by smoothing S(q,Ω) from MD using Gaussian kernel smoothing with kernel width σ = 0.01. In the right upper panel (b), dispersion relation from peaks of kernel smoothed S(q,Ω) (red dots), QLCA (blue dots), Murillo’s MNSE formulation (green dots) and VE-DDFT formulation (orange dots). In the right lower panel (c), for a wider range of q, we can see significant deviation of the theoretical dispersion relations from the MD dispersion relation for q greater than about 0.4. Lines connecting the dots in (b) and (c) are for guiding the eye. long with a time step ∆t = 0.05 ω−1 p 51 0.00.51.01.52.0⌦=!/!p0.00.51.01.52.02.5S(q,⌦)q=2⇡/L=0.18(a)q=qq=2qq=3qq=4qq=5qq=6q0.10.20.30.40.50.60.800.850.900.951.001.05⌦=2,=0.1(b)QLCAMurilloMNSEPPPM-MD(N=104)VE-DDFT0.00.20.40.60.81.01.2q=kai0.00.20.40.60.81.01.2⌦(c) where ˜y2( ˜x) is defined as ˜y2( ˜x) = (cid:113) ˜y2( ˜x)− [ ˜y( ˜x)]2 , ∆ ˜y( ˜x) = ∑B i=1 Kσ ( ˜x,xi) y2(xi) ∑B i=1 Kσ ( ˜x,xi) (3.58) (3.59) . where y(xi) is the observation of y at xi, B is the number of observations, and Kσ ( ˜x,xi) is the smoothing kernel. This yields a weighted average to give a smooth function ˜y( ˜x) evaluated at ˜x. We use the Gaussian kernel where σ is the kernel width. For Γ = 50 and κ = 1, Figs. 3.7 and 3.8 show the smooth dynamic structure fac- Kσ ( ˜x,xi) = e( ˜x−xi)2/2σ2 , (3.57) tor curves (red curves) obtained by kernel smoothing of noisy dynamic structure factor curves obtained from MD simulations (grey curves) with σ=0.01 and σ=0.005, respectively. MD simulations of 103 unit charges were performed using the PPPM algorithm, with an RMS force error of approximately 10−4 (e2/a2 parameters used were α = 1, rc = 3, M = 64, and p = 2. i ). The We have defined a simple way of quantifying the smoothing error as follows: More rigorous ways of quantifying the error in smoothing can be found in [[196]] and [[85]]. Figs. 3.7 and 3.8 show the error in smoothing (blue band) about the smooth S(q,Ω) curves (red curves). Because the kernel width influences the extent of smoothing and the corresponding error in the smoothing, the width must be chosen carefully such that the smoothed data are consistent with the physically expected result. 3.4.3 Comparison of the theoretical and MD dispersion relations In this section, we compare dispersion relations obtained using PPPM-MD with those obtained using dif- ferent theoretical models. For convenience, k and ω are expressed as dimensionless quantities q = kai and Ω = ω/ωp, where ωp = (4πnie2/M)1/2 is the ion plasma frequency for ionic mass M and ionic number density ni. One simple dispersion relation obtained from a static structure factor S(q) [147] is given by In the random-phase approximation (RPA) for a Yukawa system, S(q) takes the approximate form Ω2 = q2 3ΓS(q) . SRPA(q) = q2 + κ2 q2 + κ2 + 3Γ . 52 (3.60) (3.61) More accurate theoretical forms of S(q) are discussed in [[147]] and [[65]]. Fig. 3.9a shows the com- parison of a dispersion relation obtained using PPPM-MD and dispersion relations obtained using different theoretical approximations. The purpose of this comparison is to show that the MD result agrees well with the theoretical dispersion relations in the low-q limit, where the theory is considered to work well [147]. As q increases, the MD curve begins to deviate from the theoretical curves; this occurs because of the lack of required many-body effects in those theoretical models. A very accurate S(q) can be obtained from hypernetted chain (HNC) equations with a bridge function [58] (denoted by SHNC(q)). Although particle correlations are captured in Ω(q) by SHNC(q) (solid blue curve), this model does not incorporate viscoelastic effects [65], and hence, it underestimates Ω(q) compared to MD. Fig. 3.9b shows the full- width-half-maximum (FWHM) of a smooth S(q,Ω) constructed by kernel smoothing of S(q,Ω) obtained from MD data using a Gaussian kernel of widths σ =0.01 (blue dots) and σ =0.005 (green dots). A smaller width of σ =0.005 results in a curve that approaches zero (as expected from theory) more accurately than does a curve constructed using a width of σ =0.01. More rigorous theoretical models for computing the dispersion relation are the quasilocalized charge ap- proximation (QLCA) [111], modified Navier-Stokes (MNSE) [147] and viscoelastic-dynamic density func- tional theory (VE-DDFT) [65]. Using QLCA, the dispersion relation takes the form Ω2 = q2 q2 + κ2 + dr r [g(r)− 1] K (qr,κr) , (cid:90) ∞ 0 where K (x,y) is given by K (x,y) = −e−y − e−y y2 3 (cid:20)(cid:18) (cid:18)sinx 2 3 2 + 2y + (cid:19) x − 1 (cid:19)(cid:21) 3cosx x2 − 3sinx x3 (cid:19)(cid:18)sinx x + y2 . The dispersion relation from the Navier-Stokes equation modified to incorporate low-frequency corrections (3.62) (3.63) (3.64) results in the dispersion relation for MNSE given by (cid:33)2 , q2(cid:101)η 4√3ωE 3ωp Ω2 = q2 3ΓS(q) − where ωE is the Einstein frequency, and (cid:101)η is the dimensionless viscosity given by (cid:101)η = η/√3ωEMnia2 where η is the shear viscosity. Values of (cid:101)η can be found in [[184]]. VE-DDFT incorporates viscoelastic i , (cid:32) 53 effects that are not incorporated in QLCA and MNSE formulations. The dispersion relation calculated using VE-DDFT is given by where τ is given by Ω2 = q2 3ΓS(q) − i ¯ηq2Ω 1− iΩτωp , (3.65) (3.66) ωpτ = ; 15Ec 3 ¯ηΓ 1− γiµ + 4 (cid:16)4 3η + ε (cid:17) here, γi is the adiabatic index, µ is the compressibility, Ec(Γ,κ) is the correlation energy in units of temper- ature, and ¯η is the normalized viscosity given by ¯η = i ), where ε is the bulk viscosity. S(q,Ω) curves from MD were obtained by averaging over data from 20 different PPPM-MD simulations /(Mniωpa2 for Γ = 2 and κ = 0.1 using 104 unit charges and an RMS force error of 10−6 (e2/a2 i ). The parameters used long with a time step ∆t = 0.05 ω−1 were α = 0.46, rc = 7.8, M = 64, and p = 6. Each MD run was 800 ω−1 p p and frequency resolution ∆ω (cid:39) 0.008 ωp. S(q,Ω) curves shown in Fig. 3.10a were obtained by smoothing the noisy S(q,Ω) curves from MD simulations using Gaussian kernels of width σ = 0.01. Figs. 3.10b-c show the corresponding dispersion relations obtained using the smoothed S(q,Ω) curves and comparisons with the theoretical dispersion relations obtained for QLCA, MNSE and VE-DDFT, with S(q) given by SHNC(q). We can see that for values of q greater than approximately 0.3, the MD result increasingly deviates from the theoretical curves, while the VE-DDFT curve deviates from the MNSE and QLCA curves for q greater than approximately 0.4. Further, the QLCA and MNSE curves also deviate from each other for q greater than approximately 0.7. Hence, the differences between the dispersion relations for the theoretical models and their differences from MD results clearly demonstrate the need for model validation with accurate and computationally efficient MD simulations, which can be performed using the PPPM algorithm as discussed in this paper. 3.5 Summary and conclusions We have generalized the Ewald decomposition to screened Coulomb interactions with arbitrary DRFs. The resulting effective interactions span from short to long range, depending on the physics contained in the DRF. This extension allows us to identify parameters for which the decomposition offers a computational 54 advantage over efficient neighbor-list-only algorithms (e.g., LCL-only), thereby providing an operational definition of “medium range." We have demonstrated an elegant way of performing the Ewald decomposition in Fourier-space and applied the decomposition to two choices of the DRF, Yukawa and EGS interactions, that span a wide range of temperatures and densities of interest. However, this Ewald decomposition is quite general and can be applied to non-Coulomb interactions, such as those that employ electron-ion pseudopotentials [63]. For the Yukawa interaction, we implemented a standard PPPM [101, 201] algorithm. We verified our implementation with an extensive analysis of force errors that are inherent to the PPPM algorithm. The analysis involved deriving, for a random distribution of particles, analytical error estimates that agreed very well with the computed error. The error estimates provide insight into the effects of truncation and discretization on the errors in the computed forces, and this information is needed when choosing optimal parameters for calculating forces efficiently. The error analysis was also extended to the plain Yukawa interaction, which allows for a cutoff radius because of its exponentially decaying form and, hence, can be implemented using the LCL algorithm. For the same force errors computed using the PPPM and LCL algorithms, we compared the time required for one step of a force calculation for a range of κ values and particle numbers. We found that the PPPM algorithm offers a significant computational advantage over the LCL algorithm for a medium-range Yukawa potential, which corresponds to screening lengths exceeding the ion-sphere radius (κ < 1). For screening lengths less than the ion-sphere radius (κ > 1), the computational advantage of PPPM over LCL was found to diminish (see Table 3.1). As screening lengths decreased further to values less than half of the ion-sphere radius (κ ≥ 2), the LCL algorithm became more advantageous than the PPPM algorithm. We then applied our method to the computation of the dynamic structure factor for a Yukawa system. By varying the screening length and the system size, we found cases for which the PPPM algorithm is computationally advantageous over the LCL algorithm. One such case is κ = 0.1 and 104 unit charges, corresponding to an RMS force error of approximately 10−6 (e2/a2 i ), which precludes the use of the LCL algorithm because the minimum image convention is not satisfied. The computed dynamic structure fac- tor tends to be noisy and, therefore, must be smoothed, which is usually accomplished by averaging over the results of MD runs that vary in their initial conditions. We have presented an alternative and potentially 55 more efficient way of smoothing S(k,ω) data using kernel smoothing. With Gaussian smoothing kernels, the dispersion relation obtained from the smoothed dynamic structure factor curves were compared with theo- retical dispersion relations obtained using the QLCA, MNSE and VE-DDFT models, as shown in Fig. 3.10. The differences between the dispersion relations obtained using MD and those obtained using theoretical models, as shown in Fig. 3.10, highlight the need for model validation with accurate and computationally efficient MD simulations, which are possible using the PPPM algorithm for a wide range of temperatures and densities of interest. We note that a number of techniques, such as analytical differentiation [10], interlacing [153], oversam- pling [154], and others could be used to design PPPM algorithms. The particular choice of the technique has consequences for the form of the optimal Green function, the errors and the computational costs. Any of these techniques could be applied in a straightforward way to arbitrary dielectric response functions. 56 CONTROLLED NON-IDEAL PLASMAS FROM PHOTOIONIZED COMPRESSED GASES CHAPTER 4 The physics of non-ideal plasmas is important in a wide range of applications, including exploding wires [20, 62, 96], late-stage stellar evolution [31, 23, 32, 162, 104], interiors of giant planets [82], warm dense matter (WDM) experiments [214, 198, 80, 91], inertial-confinement fusion (ICF) experiments [54], dusty plasmas [207, 165, 155, 156, 211, 99], condensed matter [95, 77, 209], non-neutral plasmas [141, 107, 69, 2], and ultracold neutral plasmas (UCNPs) [112, 125, 195, 39, 97, 113, 37, 100, 13, 203]. Unfortunately, it is not currently possible to control the physical state of a plasma and simultaneously diagnose the same plasma using any existing experimental platform. For example, dusty plasmas tend to form monolayers and are difficult to create in the moderately coupled regime [207]; similarly, WDM is generated at very high pressures, which yield highly transient states that are difficult to probe without large-scale facilities [91]. UCNPs require extremely low temperatures and are difficult to produce in the strongly coupled regime [150, 146]. Based on a large suite of molecular dynamics (MD) simulations, we propose an experimental procedure for creating three-dimensional non-ideal plasmas with controllable properties at energy densities between those of UCNPs and those of WDM. By operating at a wide range of densities lower than solid, physical scales of interest (e.g., collisional and collective) are effectively stretched (relative to WDM experiments [214, 198, 80]), and a wider range of diagnostics becomes available. The paradigm builds upon previously successful examples of similar methods applied to UCNPs that employ a cold, photoionized gas at very low density [113]. UCNP experiments, however, are plagued by the process of disorder-induced heating (DIH) [150], which prevents strong coupling despite plasma creation from ultracold gases (Tion ∼ 10−6K). DIH occurs as a consequence of creating a plasma in a non-equilibrium state; the excess potential energy in the system is rapidly converted into kinetic energy. Mitigating DIH requires a precorrelated state (prior to ionization) [150], and strong coupling between ions has been shown to be theoretically possible in degenerate Fermi gases [150], Rydberg-blockaded gases [14] and optical lattices [151], all of which require very low temperatures. Here, we propose a novel paradigm for precorrelation: compressing gases to high pressures at room temperature, thereby reducing the burden of trapping and cooling a gas. As with other plasmas produced through photoionization, this procedure allows the electron temperature to be tuned 57 independently of the ion temperature, thereby permitting two-temperature studies in the non-ideal regime. Because pressure can be continuously varied, this approach allows the systematic variation of plasma conditions in order to study a wide range of phenomena, including transport and relaxation, ionization- potential depression (IPD), the non-ideal multi-temperature equation of state (EOS), and screening involv- ing strongly coupled electrons. Transport properties [221, 139] and relaxation processes [90, 149, 92, 41] are important to quantify for modeling ICF and WDM experiments. IPD is the lowering of energy levels (“continuum lowering") of ions embedded in dense plasma environments [70, 202, 55]; the advent of x-ray free-electron lasers (XFEL) has renewed interest in IPD, and WDM experiments have been used to validate different IPD models [43, 216, 42, 79]. An accurate EOS is essential for examining plasma evolution on the hydrodynamic scale, such as the evolution of fluid instabilities in ICF experiments [131, 22] or UCNP expansion into a vacuum [125]. Strongly coupled electrons, which are difficult to create at high densi- ties because of degeneracy, have been of theoretical and experimental interest for several decades [95, 77], including, more recently, for their role in the conductance of liquid-helium surfaces [209]. Finally, our ap- proach could also serve as a platform for measuring coupled electron-ion modes, which have been of recent interest in dense two-temperature plasmas [21]. A few of these plasma properties have been examined in isolation [128, 15, 16]; here, we present a unified approach for studying all of these properties. 4.1 Mitigating DIH by photoionizing a precorrelated neutral state at room temperature Our procedure is as follows. Because of DIH, cooling a gas before ionization has only modest advantages, so we consider plasmas formed at room-temperature. To mitigate DIH from that state, we propose pre- correlating the neutral gas through strong neutral-neutral interactions. This is most easily accomplished by increasing the gas pressure to achieve various levels of precorrelation before ionization; photoionization is then used to form a two-temperature plasma. In addition, the plasma formation and photoionization must be faster than the nucleation time for liquid-droplet formation. For conciseness, we will refer to a plasma formed through this procedure as a “pressure-induced precorrelation plasma" (PIPP). We have chosen to study neutral xenon (Xe) gas in this work, owing to the large size of the Xe atom, which optimally facilitates precorrelation at increasing pressures. We consider Xe at a temperature of 300K, which is above its critical 58 λ 2 e ≈ T 2e + ( 2 3 EF )2/(4πe2ne), (4.2) where Te is the electron temperature (in energy units), and EF = ¯h2(3π2ne)2/3/2me is the electron Fermi energy for electron density ne. We choose Zi = (cid:104)Z(cid:105) for all ions, where (cid:104)Z(cid:105) is the mean ionization calculated using a non-ideal Saha equation [182, 183, 223]; see the appendix for details. The temperatures before and after DIH are related to the correlation in the initial and final states by the following relation [150]: (cid:114) (cid:90) point [226], thereby avoiding liquid-droplet formation. DIH is a non-equilibrium phenomenon involving very strongly correlated ions; thus, MD simulations are well suited to explore this process. A sufficient MD model for cool, dilute plasmas [150] utilizes the Yukawa pair potential, given in Gaussian-cgs units by uY αγ (r) = Zα Zγ e2 r e−r/λe, (4.1) where r is the distance between two ions, Zα is the charge number of the αth ion, e is the fundamental charge, λe is the degeneracy-corrected electron-screening length [200] approximated by Tf = T0 + 1 3 dkkk (2π)3 v(k) ε(k) [S0(kkk,T0)− Sf(kkk,Tf)] , (4.3) where T0 and Tf are the initial and final temperatures, S0(kkk,T0) and Sf(kkk,Tf) are the initial and final structure factors, respectively, v(k) = 4πe2/k2 is the Coulomb potential in Fourier space, and ε(k) is the dielectric response function. Thus, to achieve large ion-ion coupling after DIH, an initial state with correlations similar to those of the desired plasma is needed. In a PIPP, precorrelations in the gas are controlled by varying the pressure. MD simulations using 104 Xe atoms interacting through the exponential-6 (E6) potential [178, 19, 159] (cid:20) eη(1−r/σ ) − η 6 (cid:17)6(cid:21) (cid:16)σ r uE6 αγ (r) = 6ε η − 6 (4.4) were performed over several orders of magnitude in pressure, where ε/kB = 243.1K, η = 13, and σ = 4.37 rrrA [159]. After equilibrating the gas, the E6 interactions were replaced by Yukawa potentials to model photoionization; the plasma was then evolved to equilibrium. Correlations in the initial state were quantified through the atom-atom radial distribution function (RDF) g(r). Results for four different initial gas pressures P = {50, 62, 143, 1393} atm are shown in Fig. 4.1. In the left four panels, we show the initial (gray) and final (colored) RDFs. It is evident that higher pressures 59 Figure 4.1 Impact of precorrelation. In the left four panels, we show the impact of precorrelating a neutral gas on the DIH process. The gray curves are the MD-predicted g(r) for a dense xenon gas (E6 potential) at four different pressures. The colored curves show g(r) for a plasma (Yukawa potential) following DIH; note that the g(r) curves at high pressures are not substantially different from those at low pressures prior to ionization. For the four different pressures, the right panel shows the impact of precorrelation on the coupling parameter, using several definitions of the parameter (dotted line, Γii; dashed line, Γs ii ); the colors of the curves in the right panel correspond to the colors of the g(r) curves in the left four panels. ii ; and solid line, Γeff ii; dashdot line, Γis result in increased correlation in the initial gas, which in turn leads to smaller differences in g(r) between the gas and plasma states. Thus, DIH can be controlled by varying the initial gas pressure, thereby allowing the Coulomb coupling of the resulting PIPP to be controlled. 4.2 Impact of precorrelation on the initial microfield distribution As discussed in the previous section, initial precorrelation in gii(r) of the dense gas leads to stronger corre- lations in the resulting plasma. Following the analysis of Ref. [146], we introduce the microfield distribution (cid:42) N δ ∑ j=1 (cid:17)(cid:43) (cid:16) FFF −FFF j(0) P(F) = 4πF2 , (4.5) which represents the distribution of the forces FFF (in units of e2/a2 i ) on the initial energy landscape of the system immediately after ionization. We can thus describe the early-time temperature evolution as T(cid:48)(t) ≈ T(cid:48)(0) + (t/τ2)2, where T(cid:48) is the temperature in units of e2/ai, t is time in units of inverse plasma frequency ω−1 (cid:113) (cid:104)F2(cid:105), and (cid:104)F2(cid:105) =(cid:82) dFF2P(F). p , τ2 = 3/ 60 0123456780.00.51.01.52.02.5g(r)P(t = 0) = 50 atmbefore DIHafter DIH0123456780.00.51.01.52.02.5P(t = 0) = 62 atmbefore DIHafter DIH012345678r (ai)0.00.51.01.52.02.5g(r)P(t = 0) = 143 atmbefore DIHafter DIH012345678r (ai)0.00.51.01.52.02.5P(t = 0) = 1393 atmbefore DIHafter DIH01020304050t (ω−1pi)100101102103104Γii , Γsii , Γeffii , ΓisiiΓiiΓsiiΓisiiΓeffii 62 atm (green), and 143 atm (red). The corresponding(cid:112) Figure 4.2 Microfield distributions. P(F) for (cid:104)Z(cid:105) = 1 and initial pressures of 50 atm (blue), i , and 0.1 e2/a2 i , respectively. All the P(F) were converged using 105 particles. (cid:104)F2(cid:105) are 6.6 e2/a2 i , 0.2 e2/a2 In the case of PIPPs, we can process our MD data to obtain microfields immediately after ionization of the correlated gas. Fig. 4.2 shows P(F) for (cid:104)Z(cid:105) = 1 and initial pressures of 50 atm, 62 atm and 143 atm; all distributions P(F) were converged using 105 particles. As the initial pressure is increased from 50 to 143 atm, the peak of P(F) shifts towards smaller F, and the width of P(F) decreases, thereby decreasing (cid:104)F2(cid:105) from ∼ 6.6 e2/a2 i . As with other quantities, note that the change in P(F) from 50 to 62 atm is more pronounced compared with the change from 62 to 143 atm. This pronounced change in P(F) for i to ∼ 0.1 e2/a2 the 50-to-62 atm transition could be attributed to a Fisher-Widom or Kirkwood line crossing, as suggested by the change in the neutral gii(r) from monotonic to oscillatory (see Fig.4.1) [4]. 4.3 Effective coupling parameter While it is common to characterize non-ideal plasmas through the bare ionic Coulomb coupling parameter Γii = (cid:104)Z(cid:105)2e2/aiTi, PIPPs are two-temperature electron-ion mixtures that are not well characterized by Γii. Here, ai = (3/4πni)1/3 is the ion-sphere radius for ion density ni, and Ti is the ion temperature. The modern interpretation of the ion-ion coupling is the ratio of the mean potential energy to the mean kinetic energy; however, this parameter omits the effects of screening by the electrons. For this reason, an effective coupling parameter can be used to capture the true coupling of the system more accurately. An effective coupling 61 0.00.51.01.52.02.53.0F (e2/a2i)012345678P(F)P(t = 0) = 50 atmP(t = 0) = 62 atmP(t = 0) = 143 atm parameter that incorporates screening [161] is given by where κ = ai/λe(Te) is the screening parameter. A more rigorous definition is given by (cid:12)(cid:12)(cid:12)(cid:12) (cid:104)Vi(cid:105) (cid:104)Ki(cid:105) Γeff ii ≡ ii = Γiie−κ , Γs (cid:12)(cid:12)(cid:12)(cid:12) = (cid:104)Z(cid:105)2 (cid:90) 1− (1 + κ)e−κ(cid:17) (cid:16) dkkk (2π)3 3Ti v(k) ε(k,Te) [1− Sii(kkk)] , (4.6) (4.7) (4.8) which can be approximated in the strongly coupled regime by the analytic formula Γis ii = Γii κ−2, see the following subsection for more details. The time evolution of these coupling parameters is shown in the right panel of Fig. 4.1, where (cid:104)Z(cid:105) = 2 in each case. Over a range of initial gas pressures, the post-DIH bare Γii ranges from ∼ 10 to 364, which covers most of the strong-coupling regime. When screening is included via Γeff , effective couplings from ∼ 1 to 55 are obtained, with κ varying from 1.7 to 2.2 (see ii Fig. 4.3). The simpler versions vary as Γs ii ∼ 1− 47. Thus, PIPPs can be used to create ii ∼ 1− 40 and Γis much larger variations in the effective coupling strength than can their UCNP counterparts [132, 14]. Formula (4.7) gives the effective ionic coupling; however, it can be useful to examine the explicit cou- pling between the ions and electrons through an interspecies coupling parameter Γαγ = |Zα Zγ e2|/aαγ Tαγ , where Tαγ = (mα Tγ + mγ Tα )/(mα + mγ ) [64], species masses are mα and mγ , the ion-sphere radius is aii = ai, and the electron-sphere radius is aee = aei = ae = (3/4πne)1/3. Species couplings Γαγ are shown in the upper-right panel of Fig. 4.3, where we see that the largest variations in the species couplings occurs in the 50− 100 atm range. We also examined two other basic plasma properties. In the lower-left panel of Fig. 4.3, we show the variations in the screening strength κ and degeneracy parameter θ = Te/EF for the PIPPs we examined; while the screening parameter remained approximately constant, the degeneracy varied by a factor of more than three. 4.3.1 Coulomb Coupling Estimation Because PIPPs allow one to control plasma properties by controlling disorder-induced heating, various levels of non-ideality are potentially achievable. Non-ideality is typically measured in terms of a coupling 62 parameter whose definition is not unique for mixtures. In this section, we describe some of the choices and their differences, as we use several definitions to characterize PIPPs. The traditional definition of the ion-ion coupling parameter is given by Γii. While Γii was historically derived as a non-dimensional parameter that characterized the celebrated, yet simplistic, one-component plasma model [33], its modern interpretation is as the ratio of the mean potential energy (∼ (cid:104)Z(cid:105)2e2/ai) to the mean kinetic energy (∼ T ). This interpretation plays a central role in characterizing different types of plasmas, and these characterizations can be used to impact experimental design. Recently, there has been disagreement [132, 44] about what form the coupling parameter should take in non-ideal plasmas. As the ionic interactions of a real plasma will be effectively modified by both bound and free electrons, we would like to incorporate these effects into our description of the coupling parameter. For example, a common ap- ii = (cid:104)Z(cid:105)2e2 proach [132, 161] to include screening by free electrons is to use the screened Coulomb (Yukawa) potential aiT e−κ , however, it can be shown that this approximation is phe- and simply write by analogy Γs nomenological at best. Nevertheless, because of its wide use, we use this definition in the main text. Rather than augmenting Γii to capture the desired physics of our problem, as in Γs ii, our starting point will instead be to calculate the ratios of the relevant mean energies. To approximate the energetic contributions of a real plasma, we follow a derivation similar to that given by Ashcroft-Stroud [7] (AS). The general Hamiltonian of a Coulomb system can be written as H = ∑ α + e2 |rrrα −rrrα(cid:48)| − E0 + E1 + E0 p2α 2m + 1 2 ∑ α(cid:54)=α(cid:48) (cid:104)Z(cid:105)2e2 1 2 ∑ |RRRγ −RRRγ(cid:48)| γ(cid:54)=γ(cid:48) (cid:104)Z(cid:105)e2 |rrrα −RRRγ| − E1 P2γ 2M , − ∑ α,γ +∑ γ (4.9) (4.10) (4.11) (4.12) where lines (4.9-4.12) have been arranged to yield components that are individually finite in the thermody- namic limit. For example, line (4.9) is simply the Hamiltonian of a free-electron gas, Hfeg, embedded in a neutralizing uniform background. Here, ne is the mean electron number density, E0 is the self-energy of a 63 uniform background of charge density ene given by (cid:90) E0 = 1 2 and E1 is the interaction between the ions and this background given by (cid:90) drrr(cid:48) (ene)2 |rrr−rrr(cid:48)| , drrr (cid:90) E1 = ∑ γ drrr(cid:104)Z(cid:105)e2ne |rrr−RRRγ| . (4.13) (4.14) The electronic and ionic coordinates are denoted by {rrrα} and {RRRγ}, with {pppα} and {PPPγ} being their momenta and m and M being their masses, respectively. Lines (4.10-4.11) can be expressed in Fourier space as (cid:90) (cid:90) Vii = (cid:104)Z(cid:105)2 2 Vei = −(cid:104)Z(cid:105) dkkk (2π)3 v(k) [ ˆni(kkk) ˆni(−kkk)− Ni] , dkkk (2π)3 v(k) ˆni(kkk) ˆne(−kkk), (4.15) (4.16) where ˆni(kkk) and ˆne(kkk) are the ion and electron density operators, respectively, Ni is the total number of ions within the domain, and v(k) = 4πe2/k2. Following AS, we next approximate the ensemble average electron density in terms of the ion density operator through a response function as (cid:104) ˆne(kkk)(cid:105) ≈ −λ(cid:104)Z(cid:105)v(k)χ(k)ni(kkk), (4.17) where χ(k) is the response function, λ is a thermodynamic integration parameter [114], and the angle brackets on the density denote the ensemble average in this context. By then taking the thermodynamic integration over the ensemble average of the Hamiltonian, the total internal energy can be expressed as (cid:104)H(cid:105) = Ufeg +(cid:104)Vii(cid:105) +(cid:104)Vei(cid:105) +(cid:104)Ki(cid:105). (4.18) The first term is the total internal energy of the free-electron gas interacting with a neutralizing uniform background. Using the definition of the ion-ion static structure factor Sαγ (kkk) ≡ (cid:104) ˆnα (kkk) ˆnγ (−kkk)(cid:105)/(cid:112)Nα Nγ , the second term, which represents the Coulomb interactions between the ions and the other compensating uniform background, can be expressed as (cid:104)Vii(cid:105) = Ni(cid:104)Z(cid:105)2 2 (cid:90) dkkk (2π)3 vee(k) [Sii(kkk)− 1] . 64 (4.19) The third term is the potential energy of the electron-ion interactions given by (cid:104)Vei(cid:105) = (cid:104)V(cid:48)ei(cid:105) +(cid:104)V(cid:48)(cid:48)ei(cid:105), (cid:104)V(cid:48)ei(cid:105) = Ni(cid:104)Z(cid:105)2 dkkk (2π)3 v2 (cid:104)V(cid:48)(cid:48)ei(cid:105) = Ni(cid:104)Z(cid:105)2 dkkk (2π)3 v2 2 2 (cid:90) (cid:90) ee(k)χ(k), ee(k)χ(k) [Sii(kkk)− 1] , (4.20) (4.21) (4.22) where we have separated out the interaction between the electron clouds and the uniform background. Fi- nally, the last term represents the kinetic energy of the ions (cid:104)Ki(cid:105) = 3 function ε−1(k) = 1 + v(k)χ(k), we can express the total potential interactions between the ions as 2NiT . Next, by introducing the dielectric (cid:90) (cid:104)Vi(cid:105) = (cid:104)Vii + V(cid:48)(cid:48)ei(cid:105) = (cid:104)Z(cid:105)2 2 Ni dkkk (2π)3 v(k) ε(k) [Sii(kkk)− 1] . (4.23) (4.24) We can thus use the above expression and the ionic kinetic energy to compute the effective ion-ion coupling of the system: (cid:12)(cid:12)(cid:12)(cid:12) (cid:104)Vi(cid:105) (cid:104)Ki(cid:105) Γeff ii ≡ (cid:90) dkkk (2π)3 v(k) ε(k) [1− Sii(kkk)] . Alternatively, this effective coupling can be expressed in terms of real-space quantities as (4.25) (4.26) (cid:12)(cid:12)(cid:12)(cid:12) = (cid:104)Z(cid:105)2 (cid:90) 3T (cid:90) ∞ 0 Γeff ii = ni 3T drrr ueff(r) [1− gii(rrr)] , where ueff is the inverse Fourier transform of (cid:104)Z(cid:105)2v(k)/ε(k). For a central effective potential, we can write (4.27) σ (r), σ (r) = (cid:90) ∞ ueff = (cid:104)Z(cid:105)2e2 dk . 2 π 0 sin(kr) kε(k) r If we further assume the pair-correlation function to be isotropic and introduce the dimensionless variable of integration x = r/ai, the effective coupling simplifies to Γeff ii = Γii dxxσ (x) [1− gii(x)] . (4.28) The evaluation of (4.28) requires gii(x) as an input; however, we can make several reasonable approxi- mations to estimate this relation analytically. For the screened Coulomb interaction, we have σ (x) = e−κx, 65 and if we approximate the pair-correlation function using that of Stewart and Pyatt [202], which smoothly interpolates between the (Debye-Hückel) DH and ion-sphere limits, we obtain (cid:34) 1− ˜κ3 + 1 (cid:32) (cid:17)1/3 1 + ζ (cid:21) − 1 Γii κ2 (cid:20)(cid:16) Γeff ii ≈ κ ˜κ ζ = (cid:33) (cid:35) κ3 + ˜κ3 ˜κ2(κ + ˜κ) e−ζ , ˜κ2 = κ2 + 3Γii. , In the weakly coupled limit (κ2, Γii (cid:28) 1), this formula yields the DH result of whereas in the strongly coupled limit (Γii (cid:29) 1), it yields the ion-sphere approximation It should also be noted that for the case of κ = 0, Equation (4.29) greatly simplifies to ΓDH ii = , κ + Γ2 ii κ2 + 3Γii (cid:113) (cid:16) 1− (1 + κ)e−κ(cid:17) (cid:26)(cid:104) (cid:105)2/3 (3Γii)3/2 + 1 κ−2. (cid:27) − 1 , Γis ii = Γii (cid:12)(cid:12)(cid:12)κ=0 Γeff ii = 1 6 (4.29) (4.30) (4.31) (4.32) (4.33) which has the weakly coupled result of Γeff ii ∼ Γii/2. Because these forms for the screened coupling parameter are based on thermodynamic arguments that begin with the general Coulomb Hamiltonian, we also report them with Γs ii/3)1/2 and the strongly coupled limit of Γeff ii to ensure that our estimates of the ii ∼ (Γ3 accessible regimes of PIPPs are well motivated. 4.4 Application of PIPPs Having established that DIH can be partly mitigated in a controlled manner to achieve a range of coupling strengths, we turn to some of the unique applications of PIPPs, and in particular, we explore (1) IPD, (2) transport phenomena, and (3) EOS properties over a wide range of conditions. 4.4.1 Ionization potential depression We begin with ionization potential depression (IPD), which has seen renewed interest recently in the context of WDM experiments [216, 42]. The energy shifts associated with IPD modify the ionization balance in a 66 Figure 4.3 Plasma properties versus gas pressure. A wide range of physical properties, in- cluding degeneracy (θ), screening (κ), species couplings (Γαγ), Coulomb logarithms (Λαγ) and equations of state from MD simulations (P) and using a mean-field approximation ((cid:101)P), are shown as a function of the initial gas pressure for a xenon plasma. Note that by varying the gas pressure, many quantities of interest can be varied substantially. For all of the cases shown, the electron tem- perature Te ∼ 6 eV, corresponding to (cid:104)Z(cid:105) = 2; these physical properties can be controlled further by varying Te, as well. plasma, and this balance affects important properties such as opacity, the EOS, and transport coefficients; these modifications to the ionization balance become more pronounced as the coupling increases. Despite the significance of IPD physics, there is a lack of well-established models that describe the phenomenon, as has been shown by recent WDM experiments [42]. The widely used Stewart-Pyatt (SP) model of IPD was found to increasingly disagree with measurements in these WDM experiments at higher ionization levels (see the appendix for a review of the SP model). Instead, the lesser-known Ecker-Kröll (EK) model, given by ∆E j = (Z j + 1)((cid:104)Z(cid:105) + 1)1/3e2/ai for the jth ionization state [55], offered better agreement with IPD measurements for the same conditions. As shown in the bottom-center panel of Fig. 4.3, the significant differences between the SP and EK predictions of IPD for PIPP conditions are similar to those observed in WDM experiments. Thus, the lower-energy-density environments of PIPPs (relative to WDM) and their controllability over a range of densities enable these plasmas to provide an optimal platform for validating different IPD models, without the complexities associated with WDM conditions. A notable feature in the IPD predictions is the sensitivity to pressure in the 50− 62 atm range, a prop- erty shared with the other applications discussed below. Re-examining Fig. 4.1, we see that the gas g(r) undergoes a monotonic-oscillatory transition, which is either a Kirkwood or Fisher-Widom transition [5]. 67 10210310410-1100101102103ΓαγΓiiΓieΓee10210310410-510-410-310-210-1100ΛαγΛiiΛieΛee1021031040.50.00.51.0P/neTePP0eP0i∆Pexee∆Pexie+∆Pexei∆Pexii102103104P(t = 0) (atm)1234567, θθ102103104P(t = 0) (atm)0510152025IPD (eV)∆E1 (SP)∆E2 (SP)∆E1 (EK)∆E2 (EK)102103104P(t = 0) (atm.)1.51.00.50.00.51.0fP/neTefP/10fP0efP0i∆fPexee∆fPexie+∆fPexei∆fPexii/10 Comparing the final g(r) for these cases (blue and green curves), we see that pressure variations in the 50− 62 atm range have a large impact on the final plasma g(r). Thus, establishing the Kirkwood/Fisher- Widom line for the gas of interest can yield considerable controllability. 4.4.2 Transport phenomena through Coulomb logarithms Turning next to transport phenomena, which are critical to a hydrodynamic description of plasmas, PIPPs can provide a platform for exploring quantities such as transport coefficients (diffusivities, viscosities, con- ductivities, etc.) that are relevant to conditions in ICF experiments [200, 90, 22, 21]. Most transport coef- ficients tend to be inversely proportional (or proportional) to a scattering cross section, which generally is not a transferable quantity among different transport processes [200]; however, we can explore the qualita- (cid:18) (cid:16) (cid:17)2(cid:19) tive effects of variations in PIPP conditions on these cross sections by employing a Coulomb logarithm (CL) as a proxy. Here, we consider CLs of the form [200] Λαγ = 1 , where the minimum impact parameter is given by bminαγ = |Zα Zγ|e2/(2Tαγ ), and the maximum impact param- eters are given by bmax e)1/2, with λi = (Ti/(4πni(cid:104)Z(cid:105)2e2))1/2 being the ionic screening length. This form is used to yield more physically rel- evant results in the strongly coupled regime, as the more common form ln is, in fact, negative e)−1)−1/2 and bmax ii = ((λ 2 ee = bmax ei = (λ 2 bmaxαγ /bminαγ i + a2 i )−1 + (λ 2 e + a2 2 ln 1 + (cid:16) e + a2 (cid:17) bmaxαγ /bminαγ in this regime. CLs are presented in the top-center panel of Fig. 4.3, where order-of-magnitude variations are seen in Λii. Note that all three CLs are near unity or smaller, representing the most theoretically challenging regime of transport in plasmas. 4.4.3 Multi-temperature equation of state model Finally, we examine PIPPs as a platform for investigating multi-temperature equation of state (EOS) models, which are essential to the hydrodynamic description of a plasma [72]. For example, large-scale experiments involving plasmas (e.g., ICF) require an accurate EOS for their design and interpretation [131, 22, 102, 34]. Because electrons and ions are at different temperatures in most laboratory plasmas, it is important to develop an accurate two-temperature EOS. Although theoretical efforts have explored this topic over many years [186, 27, 190, 64], much less experimental work on this subject has been conducted. To understand the role that PIPPs could play in such experimental studies, we calculate final system pressures using MD 68 simulations as the initial pressure of the gas phase is varied. For comparison, we employ Debye-Hückel (DH) theory, which is commonly used to describe a number of plasma properties [15, 16, 78] because of its analytic simplicity. In our EOS studies, we included electrons explicitly in the MD through the use of quantum-statistical potentials (QSPs). We use the diffractive contribution of Hansen-McDonald [98], while the Pauli exclusion term between electrons is derived from the low-density limit of the spin-averaged Lado potential [126, 109]. Each simulation used 104 ions and 2×104 electrons ((cid:104)Z(cid:105) = 2) with an initial condition given by the ionic coordinates after DIH. Pressures calculated from the MD simulations are shown in the top-right panel of Fig. 4.3 in units of the ideal-electron pressure. Here, the electron-ion contribution to the excess pressure is shown to vary the most, and this finding could be of interest for experiments focused on electron-ion interactions. However, the total pressure, which is primarily dominated by the ideal-electron pressure, is of most interest from a hydrodynamic perspective. For comparison, partial pressures predicted by the DH approximation are shown in the bottom-right panel of Fig. 4.3, denoted by (cid:101)P (see the appendix for explicit formulas). Note that we have reduced both (cid:101)P and ∆(cid:101)Pex ii by a factor of 10 to clearly show the variations in the other partial pressures. A striking observation is that the total pressure predicted by DH is negative in all cases because of the contribution of the ion-ion excess pressure. Moreover, the electron-ion excess pressure has a sign opposite that of its MD analog. These observations lead us to conclude that the DH approximation severely breaks down for the PIPP conditions studied here. Such disparities between the DH and MD predictions call for validation of the different EOS models with experiments. 4.5 Summary and conclusions In summary, we have proposed a paradigm for creating PIPPs that uses variable-pressure neutral gases along the continuum of densities between those of UCNPs and those of WDM. We have shown through MD sim- ulations that this scheme can mitigate DIH, thereby allowing for the creation of non-ideal plasmas with tunable physical parameters well into the strongly coupled regime. Some of the experimental challenges encountered with other approaches are reduced, and alternate diagnostics, such as near-visible Thomson scattering and terahertz spectroscopy [56, 103, 225, 1, 157], may be employed at these intermediate densi- 69 ties. We find that a wealth of plasma properties, including IPD, transport and the EOS, vary smoothly across the non-ideal plasma regime. The largest variations in these properties occur in the 50− 100 atm range, which has been used in previous experiments [15, 16]. This range of initial gas pressures corresponds to a change from monotonic to oscillatory behavior for the neutral Xe RDFs, which can be viewed as cross- ing the Fisher-Widom or Kirkwood line [5]; further analysis is needed to identify the exact nature of this transition. 70 CHAPTER 5 CONCLUSIONS In this chapter, I would like to summarize the conclusions of the three different aspects of my research work that I have described in this thesis. Beginning with my work on momentum-dependent potentials. A time-dependent, quantum-mechanical method was examined for large-scale simulations of non-equilibrium systems as an alternative to more expensive methods such as TD-DFT, TD-HF etc. and which relaxes most limitations associated with the relatively fast WPMD [93] and QSP [92] methods. The focus was on the use of a classical-like, Hamiltonian- based framework based on effective momentum-dependent interactions [220] that has desirable conservation properties and the computational scaling of standard classical molecular dynamics. For simplicity, the KW MDP form [220] was employed since they have been quite successful for many atomic properties [17, 47, 49, 50, 53, 217, 129, 127, 143, 51, 227]. Their strengths and weaknesses were examined for use in non-equilibrium dense plasma simulations using four criteria. The parameters of KW MDP were first trained based on HF calculations to give accurate ground state energies for neutral atoms with atomic number less than 20; the excellent agreement with HF calculations was demonstrated. The correlation between the parameters was found and the correlation was fitted to a parabolic relation, revealing a level of transferability of these parameters to previously untrained systems. Next, the properties of excited state orbits of electron-ion pair were computed and found to have dis- parate properties for the trajectories depending on the initial angular momentum, including unusual re- ciprocating patterns. However, an important feature of the KW MDP is the stabilization of the Coulomb catastrophe for the special case of zero initial angular momentum . Because plasma electrons persistently undergo ionization and recobination events, the ionization pro- cess was examined, with the ionization energy as the quality metric for the MDP. From the fixed, previously- determined ground state parameters the first- and second-ionization energies were predicted. The first ion- ization ionization energies were found to be in good agreement with experimental values, again suggesting good transferability. However, the second ionization energies were mixed, with half accurate and half with an error as large as a factor of two. The reason for this behavior in the second ionization energy is currently unknown. 71 Finally, the continuum states responsible for electronic transport processes, such as stopping power [6] and electrical and thermal conductivity [180] were examined. The ability of the KW MDP to reproduce the MTCS versus energy for electron-ion collisions was examined. A screened interaction was chosen because it more realistically represents the dense plasma environment and because the Rutherford MTCS has pathological properties (i.e., a large impact parameter divergence) in this context. Using a WKB approach [180], the classical and quantum MTCS were computed and compared with predictions from KW MDP trajectories. KW MDP cannot yield the correct MTCS despite training the parameters, suggesting that the functional form itself is responsible. Thus, although KW MDPs were successful for ground state properties, they cannot capture scattering properties important for plasma simulations. The results presented here were for isolated atoms and ions. To quantify the implications of the findings on many-body properties of a plasma, large-scale molecular dynamics simulation are needed; this is beyond the scope of the present work and will be explored in a future work. However, the results have suggested several areas of improvement. First, the KW forms did not include a Heisenberg interaction between elec- trons, which is unphysical. Second, the functional form of the KW appears to have emerged to recover a Bohr picture of the ground state, which yields poor properties for scattering states. Third, alternate training methods should be explored, including optimizing on several properties (e.g., ground state energy and cross sections) simultaneously. Such explorations are underway and will be the subject of a future work. Moving on to my work on identifying an efficient force calculation algorithm for screened Coulomb interactions that are medium-range. The Ewald decomposition was generalized to screened Coulomb inter- actions with arbitrary DRFs. The resulting effective interactions span from short to long range, depending on the physics contained in the DRF. This extension allows us to identify parameters for which the de- composition offers a computational advantage over efficient neighbor-list-only algorithms (e.g., LCL-only), thereby providing an operational definition of “medium range." An elegant way of performing the Ewald decomposition in Fourier-space was demonstrated and the decomposition was applied to two choices of the DRF, Yukawa and EGS interactions, that span a wide range of temperatures and densities of interest. However, this Ewald decomposition is quite general and can be applied to non-Coulomb interactions, such as those that employ electron-ion pseudopotentials [63]. For the Yukawa interaction, a standard PPPM [101, 201] algorithm was implemented. The implemen- tation was verified with an extensive analysis of force errors that are inherent to the PPPM algorithm. The 72 analysis involved deriving, for a random distribution of particles, analytical error estimates that agreed very well with the computed error. The error estimates provide insight into the effects of truncation and dis- cretization on the errors in the computed forces, and this information is needed when choosing optimal parameters for calculating forces efficiently. The error analysis was also extended to the plain Yukawa in- teraction, which allows for a cutoff radius because of its exponentially decaying form and, hence, can be implemented using the LCL algorithm. For the same force errors computed using the PPPM and LCL al- gorithms, the time required for one step of a force calculation was compared for a range of κ values and particle numbers. The PPPM algorithm was found to offer a significant computational advantage over the LCL algorithm for a medium-range Yukawa potential, which corresponds to screening lengths exceeding the ion-sphere radius (κ < 1). For screening lengths less than the ion-sphere radius (κ > 1), the computational advantage of PPPM over LCL was found to diminish. As screening lengths decreased further to values less than half of the ion-sphere radius (κ ≥ 2), the LCL algorithm became more advantageous than the PPPM algorithm. The method was then applied to the computation of the dynamic structure factor for a Yukawa system. By varying the screening length and the system size, cases were found for which the PPPM algorithm is computationally advantageous over the LCL algorithm. One such case is κ = 0.1 and 104 unit charges, corresponding to an RMS force error of approximately 10−6 (e2/a2 i ), which precludes the use of the LCL algorithm because the minimum image convention is not satisfied. The computed dynamic structure factor tends to be noisy and, therefore, must be smoothed, which is usually accomplished by averaging over the results of MD runs that vary in their initial conditions. An alternative and potentially more efficient way of smoothing S(k,ω) data using kernel smoothing has been presented. With Gaussian smoothing kernels, the dispersion relation obtained from the smoothed dynamic structure factor curves were compared with theoret- ical dispersion relations obtained using the QLCA, MNSE and VE-DDFT models. The differences between the dispersion relations obtained using MD and those obtained using theoretical models, highlight the need for model validation with accurate and computationally efficient MD simulations, which are possible using the PPPM algorithm for a wide range of temperatures and densities of interest. Moving on to the final aspect of my research work included in this thesis. A novel paradigm has been proposed for creating plasmas (PIPPs) that uses variable-pressure neutral gases along the continuum of densities between those of UCNPs and those of WDM. Through MD simulations it has been shown 73 that this scheme can mitigate DIH, thereby allowing for the creation of non-ideal plasmas with tunable physical parameters well into the strongly coupled regime. Some of the experimental challenges encountered with other approaches are reduced, and alternate diagnostics, such as near-visible Thomson scattering and terahertz spectroscopy [56, 103, 225, 1, 157], may be employed at these intermediate densities. It was found that a wealth of plasma properties, including IPD, transport and the EOS, vary smoothly across the non-ideal plasma regime. The largest variations in these properties occur in the 50− 100 atm range, which has been used in previous experiments [15, 16]. This range of initial gas pressures corresponds to a change from monotonic to oscillatory behavior for the neutral Xe RDFs, which can be viewed as crossing the Fisher- Widom or Kirkwood line [5]; further analysis is needed to identify the exact nature of this transition. I would like to end the Conclusions with pointers to some problems that could warrant further studies. Beginning with MDPs, it would be interesting to approach the problem of identifying better MDP functional forms through a machine learning approach. That is, given a set of quantities that are required to be captured by the MDP model, one can imagine employing a machine learning algorithm (for e.g. symbolic regression, neural networks etc.) to learning a MDP form that has the best prediction for the quantities of interest. With rapid developments currently occurring in the field of machine learning, such an approach could transform the MDP approach to a much more effective and versatile approach. Next, with respect to the development of efficient computational methods for computing forces from medium-range interactions, it would be very useful to extend the developed method for different boundary conditions such as those that can handle interactions with a surface or expansion into vacuum. In addition, parallelizing the algorithms would be a rewarding initiative because paralleliziation could offer orders of magnitude speedups as determined by the hardware and algorithm. Such speedups will play a crucial role is extending MD to extreme length- and time-scales. Finally, with the PIPP paradigm for controllable non-ideal plasmas, in order to identify the best regions for the paradigm, more MD studies with improved potentials would be needed. Further, it would be informative to explore other elements beyond Xe, including mixtures to examine the extent of DIH mitigation and associated coupling strengths that various other precorrelated neutral states can offer. The density, pressure and species regime examined for PIPPs are just pointers for possibly a whole new subfield within the domain of non-ideal plasmas. 74 APPENDICES 75 APPENDIX A APPENDIX TO CHAPTER 3 A.1 Ewald decomposition for Yukawa and EGS potentials In this section, we provide the details of Ewald decomposition for Coulomb interactions screened with Yukawa and EGS DRFs. In general, the Ewald decomposition of a Coulomb interaction linearly screened through an isotropic response takes the form 1− ζ e−k2/4α2(cid:21) (cid:20) φ (k) = v(k) ε(k) = φR (k) + φF (k), + ζ v(k) ε(k) e−k2/4α2 (A.1) (A.2) where v(k) = 4π/k2. The long-range component φF (k) can be evaluated in Fourier-space; however, the remaining short-range component φR (k) must be inverted to real-space as 1 (cid:90) (cid:90) ∞ (2π)3 2 πr 0 dkkk v(k) ε(k) sin(kr) ε(k) dk k φR (r) = = (cid:20) 1− ζ e−k2/4α2(cid:21) 1− ζ e−k2/4α2(cid:21) (cid:20) . eikkk·rrr A.1.1 Yukawa dielectric function We first consider the Yukawa dielectric function of the form ε−1 Y (k) = k2 k2s + k2 , where ks = λ−1 distribution. Using the known result s is the magnitude of the inverse screening vector associated with the background charge (A.3) (A.4) (A.5) (A.6) (cid:90) ∞ 0 dk k sin(kr) k2 + c2 e−b2k2 = π 4 eb2c2 (cid:16) [ f (−r)− f (r)] , bc + (cid:17) r 2b with f (r) = ecrerfc 76 and erfc(z) being the usual complimentary error function, we can evaluate the above expression as ek2s /4α2(cid:20) (cid:18) ks 2α − αr φ Y R (r) = ζ 1 e−ksr + r 2r − e−ksrerfc (cid:18) e−z2 erfc(z) ∼ 1− √πz (cid:18) 1− ζ ek2s /4α2(cid:19) e−ksr 1 2z2 + φ Y R (r) ∼ eksrerfc (cid:19)(cid:21) . (cid:19) 3 4z4 − . . . (cid:19) (cid:18) ks 2α + αr , Re{z} (cid:29) 1, (cid:18) r−2e−α2r2(cid:19) (A.7) (A.8) Lastly, the charge fraction ζ can be tuned to make this result as short range as possible. Noting that the large-r behavior of (A.7) is given by r + O . (A.9) Thus, ζ = e−k2s /4α2 real-space interaction is the optimal choice for the charge fraction, which results in the final form of the (cid:19) (cid:18) ks (cid:19)(cid:21) 2α + αr φ Y R (r) = = 1 1 e−ksr + r 2r (cid:18) − e−ksrerfc eksr 2r erfc αr + (cid:20) (cid:18) ks eksrerfc (cid:19) 2α − αr ks 2α + (cid:18) αr− ks 2α (cid:19) . (A.10) e−ksr 2r erfc A.1.2 EGS dielectric function We next turn to the next-order gradient correction to the Yukawa interaction using the EGS dielectric func- tion ε−1 EGS(k) = k2(k2 s + 1 4νk2) k4s + k2(k2s + 1 4νk2) , (A.11) where expressions for ν and further generalizations can be found in [[199]]; however, the analysis is appli- cable to any dielectric function that can be expressed in this form. To simplify the problem, this dielectric function can be written in terms of partial fractions as ε−1 EGS(k) = A− k2 k2 − + k2 − A+ k2 + + k2 , k2 (A.12) 77 where we have introduced the coefficients and new inverse lengths A± = s − 1 k2 4νk2 ± √1− νk2s , k2 ± = 2 ν (cid:16) 1±√1− ν (cid:17) k2 s . (A.13) We have now reduced the calculation to evaluating two Yukawa-type dielectric functions to yield φ EGS R (r) = A−φ−(r)− A+φ+(r) φ±(r) = ek2 ζ 1 e−k±r + r 2r − e−k±rerfc . ±/4α2(cid:20) (cid:18)k±2α − αr (cid:19)(cid:21) ek±rerfc (cid:18) −/4α2(cid:19) e−k−r 1− ζ ek2 +/4α2(cid:19) e−k+r (cid:18) 1− ζ ek2 r r − A+ large-r behavior for this interaction is given by φ EGS R (r) ∼ A− (cid:18)k±2α + αr (cid:19) (A.14) (A.15) , (A.16) (A.17) Note that this result is generic even for the case in which ν > 1, where the coefficients are complex. The and there is no longer a choice in ζ that eliminates the leading-order term. However, given that k− < k+ (for ν < 1), the optimal choice is ζ = e−k2 . −/4α2 A.1.3 Optimizing the charge fraction ζ In general, there will not be an analytic expression for the optimal ζ in terms of α and the physical parame- ters, and a numerical investigation will be necessary. An approximate expression can be obtained by noting that the large-r limit of φR (r) will be dominated by the low-k behavior of the integrand. By temporarily inserting the approximate dielectric function ε−1 app(k) = k2 0 + k2 , k2 k2 0 = lim k→0 (cid:16) (cid:17) k2ε(k) , (A.18) the resulting optimal choice of charge fraction becomes ζapp = e−k2s /4α2 . For example, k0 = ks for the EGS dielectric function; this becomes increasingly valid for small ν, when the system is relatively hot and/or dense. In the particular case of the EGS dielectric function, there is an alternative approach to accelerate the convergence of the real-space interaction. By introducing two screening clouds in the Ewald decomposition, 78 the real- and Fourier-space interactions respectively become (cid:20) 1− ζ1e−k2/4α2 (cid:20) ζ1e−k2/4α2 1 − ζ2e−k2/4α2 (cid:21) 2 1 + ζ2e−k2/4α2 2 . (cid:21) φR (k) = φF (k) = v(k) ε(k) v(k) ε(k) The above analysis can then be repeated to obtain (cid:19) + α1r 2 1 + (cid:18) k±2α1 (cid:19) (cid:20) φ EGS R (r) = A−φ−(r)− A+φ+(r) (cid:18) k±2α1 − α1r (cid:19)(cid:21) ek2 ±/4α2 ek±rerfc φ±(r) = (cid:18) k±2α2 (cid:20) (cid:19)(cid:21) (cid:18) k±2α2 − α2r ek±rerfc ζ1 1 e−k±r + r 2r − e−k±rerfc ek2 ±/4α2 ζ2 2r − e−k±rerfc (cid:18) −/4α2 1 − ζ2ek2 −/4α2 1− ζ1ek2 (cid:18) +/4α2 +/4α2 1− ζ1ek2 1 − ζ2ek2 2 φ EGS R (r) ∼ A− − A+ + α2r 2 , (cid:19) e−k−r (cid:19) e−k+r r r which now has the limiting behavior (A.19) (A.20) (A.21) (A.22) . (A.23) (A.24) (A.25) (A.26) (A.27) . This leading-order term can hence be eliminated by solving the system of equations 1 = ζ1ek2 1 = ζ1ek2 −/4α2 +/4α2 1 + ζ2ek2 1 + ζ2ek2 −/4α2 2 +/4α2 2 to obtain the optimal charge fractions ζ1 = ζ2 = ek2 +/4α2 1 +k2 +/4α2 ek2 −/4α2 +/4α2 1 +k2 2 − ek2 2 − ek2 1 − ek2 2 − ek2 −/4α2 2 −/4α2 2 +k2 +/4α2 1 −/4α2 2 +k2 ek2 −/4α2 ek2 −/4α2 +/4α2 1 +/4α2 1 79 Note that this approach results in two Ewald parameters to optimize, which could prove to be advantageous. A.2 Error analysis In this section, we provide details of the derivation of the RMS error estimate for forces computed using the PPPM algorithm. We begin with the analysis for the real-space portion of the algorithm. As introduced in Section 3.3, for a system of point charges with an arbitrary pair-wise short-range interaction, the estimate of the RMS force error is given by (3.29). For the real-space part of the Ewald decomposition for Yukawa interaction, (3.29) takes the form ∆FY R = Q √NV (cid:20)(cid:90) Vc(r≤rc<∞) (cid:21)1/2 , d3r |fffY R (rrr)|2 where fffY R (rrr) is given by fffY R (rrr) = −∇rrrφY R (r) = − ∂φY R (r) ∂ r rrr r (A.28) (A.29) because of the radial nature of φY respect to κ and express φY R (r) as a sum of two terms: R (r). To simplify the analysis, we make use of the symmetry in φY R (r) with φY R (r) = VY R (r,κ) +VY R (r,−κ) , (cid:16) (cid:17)(cid:105) αr− κ 2α . ∂VY R (r,−κ) ∂ r (cid:104) e−κrerfc ∂VY R (r,κ) ∂ r (A.30) (A.31) (A.32) , eκre− αr+ κ 2α + κeκrerfc (cid:16) αr + (cid:17) κ 2α (cid:16) αr + κ 2α (cid:17) (A.33) can be approxi- (A.34) where This implies that where ∂VY R (r,κ) ∂ r − VY R (r,κ) = 1 2r ∂φY R (r) ∂ r = − −  2α (cid:16) √π 1 = − 2r eκr 2r2 erfc + (cid:16) (cid:17) αr + . κ 2α − (cid:17)2 For a sufficiently large distance (r (cid:29) 1), the complementary error function erfc mated by the first term of the asymptotic expansion [118] given by (cid:16) erfc αr + (cid:17) κ 2α (cid:16) 2α2(cid:16) e− αr+ κ 2α r + κ 2α2 (cid:17)2 (cid:17) . 2α √π (cid:39) 80 Using this approximation in (A.33) gives ∂VY R (r,κ) ∂ r − eκr 2 (cid:39) − eκr 2α√π +  2α √π (cid:16) r2(cid:16) e− (cid:16) αr+ κ 2α e− r (cid:17)2 (cid:17) . αr+ κ 2α r + κ 2α2 (cid:17)2 + κ α√π (cid:16) (cid:16) e− r αr+ κ 2α r + κ 2α2 (cid:17)2 (cid:17)  (A.35) For r (cid:29) 1, the second and third terms on the right-hand side can be neglected compared to the first term, resulting in − ∂VY R (r,κ) ∂ r = − which implies that − − Substituting (A.37) in (A.28) gives (cid:16) ∆FY R (cid:17)2 (cid:39) = as , . e−κ2/4α2 e−α2r2 r (cid:19) (cid:18) α √π (cid:19) ∂VY ∂φY (cid:18) 2α R (r,κ) (cid:39) − ∂ r ∂VY R (r,−κ) ∂ r R (r) ∂ r (cid:39) − (cid:90) V (r≤rc<∞) √π Q2 NV 16α2Q2 NV . This, in turn, yields r e−κ2/4α2 e−α2r2 (cid:32) e−κ2/2α2(cid:90) ∞ 4α2 π (cid:33) e−κ2/2α2 e−2α2r2 r2 e−2α2r2 d3r dr . rc r (cid:90) ∞ rc e−2α2r2 r dr e−2α2r2c 4α2rc , (cid:39) (A.36) (A.37) (A.38) (A.39) (A.40) Using the first term of the asymptotic expansion [118], the integral on the right-hand side is approximated resulting in the simplified expression for ∆FY R given by ∆FY R (cid:39) 2Q √NV e−κ2/4α2 e−α2r2c √rc . We note that the approximation, which involves truncation of the asymptotic expansion that led to the above error estimate, is not valid for κ/α → ∞ as α → 0. Hence, this error estimate will not tend to the error estimate for the plain Yukawa interaction that is derived below. Similar error estimates can be derived for the plain Yukawa interaction because it decays exponentially. For the plain Yukawa interaction φY (r) = e−κr r , 81 (A.41) we have e−κr r2 (1 + κr) Using (3.29), the error estimate for the Yukawa interaction is given by fffY (rrr) = ∂φY (r) rrr r = ∂ r rrr r . (A.42) (A.43) (cid:16) ∆FY(cid:17)2 ∆FY(cid:17)2 (cid:16) (cid:90) ∞ rc 4πQ2 NV = e−2κr r2 (1 + κr)2 ; (cid:19)2 + κ (cid:18) 1 rc 2πQ2 NV (cid:39) e−2κrc ∆FY Q √NV (cid:39) √2πκ e−κrc . the integral in this equation can be approximated using the first term of an asymptotic expansion [118] to yield κ rc (cid:28) κ, we could further approximate (A.44) to obtain For 1 . (A.44) (A.45) We now proceed to simplifying the estimate of the RMS error in the Fourier-space force computed using the PPPM algorithm. As introduced in Section 3.3, the RMS error estimate for the Fourier-space force is given by where ∆FF (cid:39) Q √NV 2/3 χF , (cid:104) ∑mmm ˆU2 (cid:20) kn+mkn+mkn+m ∑mmm ˆU2  (cid:21)2 ˆGkn+mkn+mkn+mknknkn ·kn+mkn+mkn+m kn+mkn+mkn+m |knknkn|2 F V 2/3 = ∑ χ2 kkk(cid:54)=0 ˆGkkk|kkk|2 −∑ nnn (cid:105)2  . (A.46) (A.47) If the arbitrary Green function (denoted by ˆGknknkn) rapidly decays in Fourier-space such that alias contributions to the Green function can be neglected (mmm denotes the alias indices), then (A.47) can be simplified further, beginning with the approximation given by (cid:20) ∑ mmm ˆU2 kn+mkn+mkn+m ˆGkn+mkn+mkn+mknknkn ·kn+mkn+mkn+m This implies that  F V 2/3 (cid:39) ∑ χ2 kkk(cid:54)=0 ˆGkkk|kkk|2 −∑ nnn 82 (cid:21)2 (cid:104) ˆU2 (cid:39) = ˆG2 ˆGknknknknknkn ·knknkn n ˆU4 knknkn . knknkn knknknk4 n ˆU4 k4 knknkn (cid:21)2 (cid:20) ˆG2 knknkn ∑mmm ˆU2 kn+mkn+mkn+m k2n (cid:105)2  , (A.48) (A.49) which can be rewritten as 1− F V 2/3 (cid:39) ∑ χ2 nnn ˆG2 knknknk2 n  , (cid:19)2 (A.50) (cid:18) ˆU4 knknkn ∑mmm ˆU2 kn+mkn+mkn+m with the assumption that the maximum length of nnn (denoted by nnnmax) is sufficiently large such that ˆG2 knknkn 0 for nnn = nnnmax. For an isotropic system, following the derivation in [[61]], we find that k2 n (cid:39) 1− (cid:18) ˆU4 kkkn ∑mmm ˆU2 kkkn+m (cid:19)2  (cid:39) 3 (cid:18)knzh (cid:19)2p p−1 2 C(p) m ∑ m=0 (cid:18)knzh (cid:19)2m 2 , (A.51) m are listed in Table I of [[61]]. Therefore, (A.50) is approximated to a form given where the coefficients C(p) by F V 2/3 (cid:39) 3∑ χ2 (cid:18)h nnn = 3 2 (cid:19)2m(cid:35) (cid:18)knzh 2 ˆG2 knknknk2 n k2(m+p) nz (cid:34)(cid:18)knzh (cid:19)2p p−1 ˆG2 knknknk2 n 2 ∑ m=0 C(p) m 2 C(p) m (cid:19)2p p−1 ∑ (cid:19)2m (cid:20) (cid:18)h m=0 ∑ nnn (cid:19)3(cid:90) (cid:18) L (cid:19)3(cid:90) (cid:18) L d3k , 2π ∑ nnn → (cid:21) . (A.52) (A.53) (A.54) When the box dimensions are large (L (cid:29) 1, for a cube of edge length L), n k2(m+p) nz → 2π d3k ˆG2 kkk k2 k2(m+p) nz = I . (A.55) If the Green function is symmetric, then the integral I can be reduced further to the following 1D integral: d3k ˆG2 kkk k2 k2(m+p) nz dφ sinφ (cosφ )2(m+p) (cid:21) (cid:20)(cid:90) ∞ 0 (cid:21) , dk ˆG2 k k2(m+p+2) which implies that I = = where ∑ nnn (cid:18) L 2π V (2π)2 ˆG2 knknknk2 (cid:19)3(cid:90) (cid:20)(cid:90) π (cid:90) π 0 dφ sinφ (cosφ )2(m+p) = 2 1 + 2(p + m) . 0 The integral I can then be expressed as I = V 2 (2π)2 [1 + 2(p + m)] β (m, p) , 83 (A.56) (A.57) (A.58) (A.59) . 4π k2 + κ2 e−(k2+κ2)/4α2 (cid:20) (cid:21)2 4π e−k2/2α2 ˆGk = (cid:90) ∞ 0 dk The integral β (m, p) for this case takes the form βY (m, p) = which results in (3.38). k2(m+p+2) , (A.64) where (cid:90) ∞ 0 β (m, p) = Therefore, for an arbitrary Green function, dk ˆG2 k k2(m+p+2) . (A.60) χ2 F V−1/3 = AF 3 (cid:39) (2π)2 (cid:19)2(m+p) (cid:18)h 2 p−1 ∑ m=0 C(p) m 2 [1 + 2(p + m)] β (m, p) , (A.61) which yields the following approximated estimate of the RMS error for an arbitrary Green function: ∆FF = Q √NV A1/2 F . The Green function for the Ewald decomposition for Yukawa interaction is given by (A.62) (A.63) k2 + κ2 84 APPENDIX B APPENDIX TO CHAPTER 4 For completeness and improved reproducibility, we give further details of several of the models used in our study. In the first section, we review the non-ideal Saha equations as a model for ionization in pressure- induced precorrelation plasmas (PIPPs). Next, we provide insights into the coupling-strength model used in the main text used to quantify the impact of precorrelation. Then, we provide an analysis of how this precorrelation impacts the initial-state microfield; this analysis reveals the dense-gas pair correlations in a slightly different, but equivalent, way. Finally, we collect results from the main text for the equation of state of a two-temperature system. B.1 Atomic Ionization Model B.1.1 Saha Equation Because we are considering compressed gases that are below solid density, we employ a non-ideal Saha equation rather than an average-atom model [200] to compute ionization levels. The framework of the Saha approach has the additional advantage that it naturally yields the distribution of charge states, not only the average. In its ideal form, the well-known Saha equations [182, 183] describing the transitions between N ionization states (in addition to the neutral state) are given by (cid:18) (cid:19) ≡ f j+1, (B.1) n j+1ne n j = 2g j+1 λ 3 eeg j exp E j Te − for j ∈ [0,N−1], where {g j} are the statistical weights of the energy states, {E j} are the ionization energies, λee = ¯h/(2πmeTe)1/2 is the thermal de Broglie wavelength, Te is the electron temperature in units of energy, n j is the number density of the jth ionization state, and ne is the electron density. Finally, conservation of charge, ne = ∑ j jn j, and conservation of particle number, ni = ∑ j n j, result in a closed system of equations; here, ni is the ion density. To solve this system, we follow the approach of Zaghloul et al. [223]. We can write (B.1) as a single equation in terms of the concentrations x j = n j/ni and the mean ionization (cid:104)Z(cid:105) = ∑ j jx j 85 as n−k i N+1 + (cid:104)Z(cid:105) N ∑ k=1 ((cid:104)Z(cid:105)− k)(cid:104)Z(cid:105) N−k k ∏ j=1 f j  = 0. (B.2) A root solver can be used to solve for (cid:104)Z(cid:105), and once (cid:104)Z(cid:105) is known, the individual species concentrations can be calculated from the following recursion relations: −1  N ∑ (cid:18) f j+1 (cid:19) k=1 ni(cid:104)Z(cid:105) k ∏ j=1 f j k nk i (cid:104)Z(cid:105)k x j. x0 = (cid:104)Z(cid:105) x j+1 = , (B.3) (B.4) For the specific case of Xe considered in the main text, we need to explore only the first few ionization states because of the limited temperature range of interest. Using the ionization energy values of Xe, we take N = 5 and set {E j} = {12,21,32,46,57} eV. We estimate the statistical weights by assuming that Xe, like other noble gases (beyond He), initially ionizes by losing electrons from the outermost p-subshell. Taking values from [223, 179], we have {g j} = {1,6,9,4,9,6}. While the gas in consideration is below solid density, it is still compressed to form Xe-Xe correlations, and this basic Saha model must be extended to its non-ideal form to capture many-body screening effects. B.1.2 Ionization-Potential Depression In a dense plasma, the environment surrounding each atom can significantly modify the energy required for ionizing transitions. These modifications are collectively referred to as ionization-potential depression (IPD) because they tend to have the effect of lowering ionization energies. To capture the effects of IPD, we employ the standard method of Stewart and Pyatt (SP) [202]. In the SP model, the electrostatic potential about a point ion is calculated using the following two assumptions: • Within the ion-sphere radius (r < ai), the ion density is zero, and the electron density is constant. (Here, ai is the ion-sphere radius defined as ai = (3/4πni)1/3.) • Beyond the ion-sphere radius (r > ai), the electrons and ions exhibit linear screening, with screening lengths λe and λi, respectively. 86 Figure B.1 Xe ionization. The left panel shows the mean ionization (cid:104)Z(cid:105) as a function of electron temperature for different Xe densities. The right panel shows the SP prediction of the IPD energies for the first (dashed line) and the second (solid line) ionization levels of Xe as a function of electron temperature. The colors in the right panel correspond to the colors in the left panel that differentiate the Xe densities for different initial pressures. Solving the associated Poisson equation and enforcing C 2-continuity in the potential results in an electro- static potential that can be radially expanded about an ion position as u(r) ∼ −Z je2/r + ∆E j. This energy shift associated with the jth ionization state is given by (cid:32)(cid:18) ai λtot (cid:19)3 (cid:33)2/3 + 1 − 1  , ∆E j = 3Z je2λ 2 tot 2a3 i where λ−2 tot = λ−2 e + λ−2 i . The ion screening length can be estimated using the DH theory to obtain (B.5) (B.6) λ 2 i = Ti . 4π(cid:104)Z(cid:105)2e2ni (cid:113) 3EF )2 T 2e + (2 4πe2ne λ 2 e ≈ To estimate the electronic screening length, we must additionally incorporate degeneracy via the Thomas- Fermi length [200], which can be approximated as , (B.7) where the Fermi energy is given by EF = ¯h2(3π2ne)2/3/2me. From this SP model, the energy levels are then shifted according to {E j} → {E j − ∆E j}. The (cid:104)Z(cid:105)-dependence in the energy shifts transforms (B.2) from a polynomial to a transcendental equation. 87 0246810Te (eV)0.00.51.01.52.02.53.03.54.0›Zfi1.835×1021/cc (P(t=0) = 50 atm)8.826×1021/cc (P(t=0) = 62 atm)1.039×1022/cc (P(t=0) = 143 atm)1.204×1022/cc (P(t=0) = 450 atm)1.421×1022/cc (P(t=0) = 1393 atm)0246810Te (eV)024681012141618IPD (eV)∆E1 (SP)∆E2 (SP) Fig. B.1 shows (cid:104)Z(cid:105) and the energy shifts for the first (dashed) and second (solid) ionization levels of Xe, for densities corresponding to the initial pressures P = {50, 62, 143, 1393} atm considered in the main text. The values of (cid:104)Z(cid:105) for pressures from 62 to 143 atm are similar for the range of Te considered. This similarity arises because the corresponding densities are similar (∼ 1022/cc); however, the 50-atm case differs from the others for Te (cid:38) 4 eV because of the lower densities for those cases (∼ 1021/cc). The corresponding energy-level shifts are shown in the right panel. The Yukawa potentials used in the molecular dynamics (MD) simulations employ these ionization states and (B.7). While it remains unclear which IPD model is most appropriate, we have employed this standard model both to self-consistently predict an accurate ionization state to use in the MD simulations and to illustrate the degree to which IPD physics can be explored with PIPPs. Our goal here is not to develop Saha models beyond this standard description. Future plasma experiments employing dense gases could provide more insights into IPD, and more refined models could then be developed. B.2 Two-temperature equation of state in the mean-field approximation Most laboratory plasmas have separate electron and ion temperatures, giving great importance to our un- derstanding of a two-temperature equation of state (EOS). Although theoretical efforts towards this end have continued over many years [186, 27, 190, 64, 168], much less has been done experimentally. To ob- tain analytical results, a mean-field approximation is often invoked in the DH limit. Here, we present a two-temperature generalization for comparison with MD simulations in the main document. In our EOS studies, we included electrons explicitly in the MD through the use of quantum-statistical potentials (QSPs) described in the main text. Typical species RDFs, denoted by gαγ (r), from MD simula- tions are shown in Fig. B.2. The total pressure of a multi-component system can be decomposed as Pex αγ P = ∑ α α + ∑ P0 (cid:90) α,γ nα nγ Pex αγ = − 6 drrr r duαγ dr gαγ (r), (B.8) (B.9) where P0 α are the ideal pressures. Substituting gαγ (r) = 1 + hαγ (r), where hαγ (r) are total correlation functions, the (individually infinite) uniform background terms cancel and the change in excess pressures 88 Figure B.2 Partial radial distribution functions versus pressure. The functions gii(r) (blue), gei(r) (green), and gee(r) (red) were computed with MD using quantum-statistical potentials for 3×104 total particles (104 ions and (cid:104)Z(cid:105)=2). Note the extreme behavior of gee(r), which exceeds unity near r = 0, revealing substantial ion-mediated effective interactions as electrons localize near (relatively) isolated ions. can instead be expressed as (cid:90) ∆Pex αγ = − nα nγ 6 drrr r duαγ dr hαγ (r). (B.10) There are some interesting features in the RDFs shown in Fig. B.2 that are worth a brief discussion. The gee(r) undergoes a Kirkwood-like transition in the 50-atm-pressure region, with its peak at the same location as the ion peak; there are strong correlations in the peak locations of gei(r) and gee(r) with those of gii(r), suggesting that the electrons are highly localized around the ions in these plasmas. Moreover, gee(0) > 1 for most of the initial gas pressures considered, indicating that an effective attraction between electrons (mediated by the ions) exists. Assuming Coulomb interactions, the change in excess pressures can be expressed as (cid:32) (cid:33) e2 r (cid:90) (cid:90) Zα Zγ nα nγ drrr hαγ (r) ∆Pex αγ = = = 1 6 1 6 1 3π Zα Zγ nα nγ Zα Zγ e2nα nγ dkkk (cid:90) ∞ (2π)3 v(k)hαγ (k) dk hαγ (k). 0 (B.11) (B.12) (B.13) These correlation functions are connected to the direct correlation functions ci j(r), through the multicom- 89 01234560.00.51.01.52.02.5gαγ(r)P(t = 0) = 50 atm.giigeigee01234560.00.51.01.52.02.5P(t = 0) = 62 atm.giigeigee0123456r (ae)0.00.51.01.52.02.5gαγ(r)P(t = 0) = 143 atm.giigeigee0123456r (ae)0.00.51.01.52.02.5P(t = 0) = 1393 atm.giigeigee ponent Ornstein-Zernike equations, which, in Fourier space, are given by hαγ (k) = cαγ (k) +∑ l nlcαl(k)hlγ (k). (B.14) In the DH approximation, cαγ (k) ≈ −βαγ uαγ (k), where uαγ (k) = 4πZα Zγ e2/k2, and βαγ = 1/Tαγ is the inverse of the mass-weighted cross temperature (in energy units) given by Tαγ = (mα Tγ + mγ Tα )/(mα + mγ ) [64]. For the system of study, the two components are ions (i) and electrons (e), where Zi = (cid:104)Z(cid:105), Ze = −1, and ne = (cid:104)Z(cid:105)ni. Within the DH approximation, the total correlation functions are then given in Fourier space by (B.15) , hii = (cid:104) (cid:104) neφ 2 niφ 2 (cid:105) (cid:105) ie − (1 + neφee)φii ie − (1 + niφii)φee 1 D 1 D hie = − D ≡ (1 + niφii)(1 + neφee)− nineφ 2 ie, φie, 1 D hee = , (B.16) (B.17) (B.18) where we have introduced the function φαγ (k) = βαγ uαγ (k). Defining B = βiβe − β 2 these relations into (B.13) to obtain ie, we can substitute (cid:32) (cid:33) ∆Pex ii = − ∆Pex ie = ∆Pex 2π(cid:104)Z(cid:105)4e4n2 i 3(k+ + k−) ei = − (cid:34) βe +(cid:104)Z(cid:105)βi ± 4πe2(cid:104)Z(cid:105)niB , , k−k+ βi + 2π(cid:104)Z(cid:105)3e4n2 i βie 3(k+ + k−) (cid:113) (βe −(cid:104)Z(cid:105)βi)2 + 4(cid:104)Z(cid:105)β 2 ie (B.19) (B.20) (B.21) (cid:35) . becomes negative for B < 0, which leads to where the roots k± are defined as It should be noted that k2 ± k2 ± = 2πe2(cid:104)Z(cid:105)ni is always real; however, the root k2 − a spurious result for the correlation functions and consequently the pressures. As we have defined Ti j, this spurious range is encountered for the temperature ratios (me/mi)2 < Te/Ti < 1, where we have assumed that me < mi. Given that Te (cid:29) Ti for our systems of interest, this regime will be avoided naturally. Of course, in the single-temperature limit, where the temperature ratio is unity, we have β = βe = βi = βie, and the correlation functions can be expressed simply as hαγ (k) = −Zα Zγβ v(k) 1 + ni(cid:104)Z(cid:105)((cid:104)Z(cid:105) + 1)β v(k) . 90 (B.22) Finally, the total pressure can thus be calculated as (cid:34) P = niTi +(cid:104)Z(cid:105)niTe+∆Pex tot, βi(cid:104)Z(cid:105)2 + 2(cid:104)Z(cid:105)βie + βe + 4π(cid:104)Z(cid:105)2((cid:104)Z(cid:105) + 1) ∆Pex tot = − 2π(cid:104)Z(cid:105)2e4n2 i 3(k+ + k−) e2niB k+k− (cid:35) . In the single-temperature limit, the total excess pressure reduces to ∆Pex tot = − λDH = T , (cid:104) 24πλ 3 DH 4π(cid:104)Z(cid:105)e2niβ (1 +(cid:104)Z(cid:105)) (cid:105)−1/2 , (B.23) (B.24) (B.25) (B.26) (cid:16) 4π(cid:104)Z(cid:105)2e2niβi (cid:17)−1/2 which is the standard DH result. Similarly, in the limit Te (cid:29) Ti, the leading-order excess pressure is given by the ion-ion correlations in the form ∆Pex . tot ∼ −Ti/(24πλ 3 i ), where λi = While the result in the main document is based on MD calculations of the pressure, we have included the DH result here for two reasons. First, the DH provides a point of comparison that answers the question: are PIPPs in an interesting EOS regime or well described by DH-level correlations? The conclusion of the main text is clear: the two methods predict total pressures with the opposite sign and disparate magnitudes, as shown in Fig. 2 of the main text. Second, the DH model has the advantage of ease of use and is therefore widely used to make estimates for experiments, and thus, the comparison is needed to quantify the error made in such estimates. 91 BIBLIOGRAPHY 92 BIBLIOGRAPHY [1] [2] [3] [4] [5] [6] [7] [8] [9] Amit Agrawal, Tatsunosuke Matsui, Wenqi Zhu, A Nahata, and ZV Vardeny. Terahertz spectroscopy of plasmonic fractals. Physical review letters, 102(11):113901, 2009. F Anderegg, DHE Dubin, TM O’Neil, and CF Driscoll. Measurement of correlation-enhanced colli- sion rates. Physical review letters, 102(18):185001, 2009. Xavier Andrade, Joseba Alberdi-Rodriguez, David A Strubbe, Micael JT Oliveira, Fernando Nogueira, Alberto Castro, Javier Muguerza, Agustin Arruabarrena, Steven G Louie, Alán Aspuru- Guzik, et al. Time-dependent density-functional theory in massively parallel computer architectures: the octopus project. Journal of Physics: Condensed Matter, 24(23):233202, 2012. Andrew J Archer, D Pini, Robert Evans, and L Reatto. Model colloidal fluid with competing interac- tions: Bulk and interfacial properties. The Journal of chemical physics, 126(1):014104, 2007. Andrew J Archer, D Pini, Robert Evans, and L Reatto. Model colloidal fluid with competing interac- tions: Bulk and interfacial properties. The Journal of chemical physics, 126(1):014104, 2007. NR Arista and P Sigmund. Stopping of ions based on semiclassical phase shifts. Physical Review A, 76(6):062902, 2007. NW Ashcroft and D Stroud. Theory of the thermodynamics of simple liquid metals. Solid State Physics, 33:1–81, 1978. Andrew D Baczewski, Luke Shulenburger, Michael P Desjarlais, Stephanie B Hansen, and Rudolph J Magyar. X-ray thomson scattering in warm dense matter without the chihara decomposition. Physical review letters, 116(11):115004, 2016. V Ballenegger, JJ Cerda, O Lenz, and Ch Holm. The optimal p3m algorithm for computing electro- static energies in periodic systems. The Journal of chemical physics, 128(3):034109, 2008. [10] Vincent Ballenegger, Jan Josep Cerdà, and Christian Holm. Removal of spurious self-interactions in particle–mesh methods. Computer Physics Communications, 182(9):1919–1923, 2011. [11] W Bang, BJ Albright, PA Bradley, DC Gautier, S Palaniyappan, EL Vold, MA Santiago Cordoba, CE Hamilton, and JC Fernández. Visualization of expanding warm dense gold and diamond heated rapidly by laser-generated ion beams. Scientific reports, 5, 2015. [12] Woosuk Bang, Brian James Albright, Paul Andrew Bradley, Erik Lehman Vold, Jonathan Carl Boettger, and Juan Carlos Fernández. Linear dependence of surface expansion speed on initial plasma temperature in warm dense matter. Scientific reports, 6:29441, 2016. [13] G Bannasch, J Castro, P McQuillen, T Pohl, and TC Killian. Velocity relaxation in a strongly coupled plasma. Physical review letters, 109(18):185008, 2012. [14] G Bannasch, TC Killian, and T Pohl. Strongly coupled plasmas via rydberg blockade of cold atoms. Physical review letters, 110(25):253003, 2013. [15] A Bataller, B Kappus, C Camara, and S Putterman. Collision time measurements in a sonoluminesc- ing microplasma with a large plasma parameter. Physical review letters, 113(2):024301, 2014. 93 [16] A Bataller, GR Plateau, B Kappus, and S Putterman. Blackbody emission from laser breakdown in high-pressure gases. Physical review letters, 113(7):075001, 2014. [17] WA Beck, L Wilets, and MA Alberg. Semiclassical description of antiproton capture on atomic helium. Physical Review A, 48(4):2779, 1993. [18] WA Beck, L Wilets, and MA Alberg. Semiclassical description of protonium formation in antiproton collisions with molecular hydrogen. Physical Review A, 74(5):052706, 2006. [19] Anatoly B Belonoshko, O LeBacq, Rajeev Ahuja, and Börje Johansson. Molecular dynamics study of phase transitions in xe. The Journal of chemical physics, 117(15):7233–7244, 2002. [20] JF Benage Jr, WR Shanahan, EG Sherwood, LA Jones, and RJ Trainor. Measurement of the electrical resistivity of a dense strongly coupled plasma. Physical Review E, 49(5):4391, 1994. [21] Lorin X Benedict, Michael P Surh, Liam G Stanton, Christian R Scullard, Alfredo A Correa, John I Castor, Frank R Graziani, Lee A Collins, Ondˇrej ˇCertík, Joel D Kress, et al. Molecular dynam- ics studies of electron-ion temperature equilibration in hydrogen plasmas within the coupled-mode regime. Physical Review E, 95(4):043202, 2017. [22] R Betti, VN Goncharov, RL McCrory, and CP Verdon. Growth rates of the ablative rayleigh–taylor instability in inertial confinement fusion. Physics of Plasmas, 5(5):1446–1454, 1998. [23] Lars Bildsten and David M Hall. Gravitational settling of 22ne in liquid white dwarf interiors. The Astrophysical Journal Letters, 549(2):L219, 2001. [24] José A Bittencourt. Fundamentals of plasma physics. Springer Science & Business Media, 2013. [25] David H Boal and James N Glosli. Computational model for nuclear reaction studies: Quasiparticle dynamics. Physical Review C, 38(6):2621, 1988. [26] David H Boal and James N Glosli. Quasiparticle model for nuclear dynamics studies: Ground-state properties. Physical Review C, 38(4):1870, 1988. [27] David B Boercker and Richard M More. Statistical mechanics of a two-temperature, classical plasma. Physical Review A, 33(3):1859, 1986. [28] Christoph Bostedt, Sébastien Boutet, David M Fritz, Zhirong Huang, Hae Ja Lee, Henrik T Lemke, Aymeric Robert, William F Schlotter, Joshua J Turner, and Garth J Williams. Linac coherent light source: the first five years. Reviews of Modern Physics, 88(1):015007, 2016. [29] Silvana Botti, Arno Schindlmayr, Rodolfo Del Sole, and Lucia Reining. Time-dependent density- functional theory for extended systems. Reports on Progress in Physics, 70(3):357, 2007. [30] Erik G Brandt and Olle Edholm. Dynamic structure factors from lipid membrane molecular dynamics simulations. Biophysical journal, 96(5):1828–1838, 2009. [31] Edward F Brown. Nuclear heating and melted layers in the inner crust of an accreting neutron star. The Astrophysical Journal, 531(2):988, 2000. [32] Edward F Brown, Lars Bildsten, and Philip Chang. Variability in the thermal emission from accreting neutron star transients. The Astrophysical Journal, 574(2):920, 2002. 94 [33] SG Brush, HL Sahlin, and E Teller. Monte carlo study of a one-component plasma. i. The Journal of Chemical Physics, 45(6):2102–2118, 1966. [34] Laurent Caillabet, Benoit Canaud, Gwenaël Salin, Stéphane Mazevet, and Paul Loubeyre. Change in inertial confinement fusion implosions upon using an ab initio multiphase dt equation of state. Physical review letters, 107(11):115004, 2011. [35] Federico Canova and Luca Poletto. Optical Technologies for Extreme-Ultraviolet and Soft X-ray Coherent Sources, volume 197. Springer, 2015. [36] Richard Car and Mark Parrinello. Unified approach for molecular dynamics and density-functional theory. Physical review letters, 55(22):2471, 1985. [37] Jose Castro, Patrick McQuillen, and TC Killian. Ion acoustic waves in ultracold neutral plasmas. Physical review letters, 105(6):065004, 2010. [38] DA Chapman and Dirk O Gericke. Analysis of thomson scattering from nonequilibrium plasmas. Physical review letters, 107(16):165004, 2011. [39] YC Chen, CE Simien, S Laha, P Gupta, YN Martinez, PG Mickelson, SB Nagel, and TC Killian. Electron screening and kinetic-energy oscillations in a strongly coupled plasma. Physical review letters, 93(26):265003, 2004. [40] Z Chen, B Holst, SE Kirkwood, V Sametoglu, M Reid, YY Tsui, V Recoules, and A Ng. Evolution of ac conductivity in nonequilibrium warm dense gold. Physical review letters, 110(13):135001, 2013. [41] BI Cho, T Ogitsu, K Engelhorn, AA Correa, Y Ping, JW Lee, LJ Bae, D Prendergast, RW Falcone, and PA Heimann. Measurement of electron-ion relaxation in warm dense copper. Scientific reports, 6:18843, 2016. [42] O Ciricosta, SM Vinko, B Barbrel, DS Rackstraw, TR Preston, T Burian, J Chalupsk`y, Byeoung Ick Cho, H-K Chung, GL Dakovski, et al. Measurements of continuum lowering in solid-density plasmas created from elements and compounds. Nature communications, 7:ncomms11713, 2016. [43] O Ciricosta, SM Vinko, H-K Chung, B-I Cho, CRD Brown, T Burian, J Chalupsk`y, K Engelhorn, RW Falcone, C Graves, et al. Direct measurements of the ionization potential depression in a dense plasma. Physical review letters, 109(6):065002, 2012. [44] Jean Clérouin, Gregory Robert, Philippe Arnault, Joel D Kress, and Lee A Collins. Behavior of the coupling parameter under isochoric heating in a high-z plasma. Physical Review E, 87(6):061101, 2013. [45] Kyle R Cochrane, Raymond W Lemke, Z Riford, and John H Carpenter. Magnetically launched flyer plate technique for probing electrical conductivity of compressed copper. Journal of Applied Physics, 119(10):105902, 2016. [46] [47] James S Cohen. Quasiclassical effective hamiltonian structure of atoms with z= 1 to 38. Physical Review A, 51(1):266, 1995. James S Cohen. Molecular effects on antiproton capture by h2 and the states of p¯ p formed. Physical Review A, 56(5):3583, 1997. 95 [48] [49] [50] [51] [52] [53] James S Cohen. Extension of quasiclassical effective hamiltonian structure of atoms through z= 94. Physical Review A, 57(6):4964, 1998. James S Cohen. Multielectron effects in capture of antiprotons and muons by helium and neon. Physical Review A, 62(2):022512, 2000. James S Cohen. Preliminary results for capture of negative muons and antiprotons by noble-gas atoms. Hyperfine interactions, 138(1):159–166, 2001. James S Cohen. Reexamination of over-the-barrier and tunneling ionization of the hydrogen atom in an intense field. Physical Review A, 64(4):043412, 2001. James S Cohen. Comment on ?laser-assisted formation of antihydrogen? 67(1):017401, 2003. Physical Review A, James S Cohen. Capture of negative exotic particles by atoms, ions and molecules. Reports on Progress in Physics, 67(10):1769, 2004. [54] National Research Council, Plasma Science Committee, et al. Frontiers in High Energy Density Physics: The X-Games of Contemporary Science. National Academies Press, 2003. [55] BJB Crowley. Continuum lowering–a new perspective. High Energy Density Physics, 13:84–102, 2014. [56] Maxime Ben Dahan, Ekkehard Peik, Jakob Reichel, Yvan Castin, and Christophe Salomon. Bloch oscillations of atoms in an optical potential. Physical Review Letters, 76(24):4508, 1996. [57] Tom Darden, Darrin York, and Lee Pedersen. Particle mesh ewald: An n log (n) method for ewald sums in large systems. The Journal of chemical physics, 98(12):10089–10092, 1993. [58] William Daughton, Michael S Murillo, and Lester Thode. Empirical bridge function for strongly coupled yukawa systems. Physical Review E, 61(2):2129, 2000. [59] Leonardo de Ferrariis and Néstor R Arista. Classical and quantum-mechanical treatments of the energy loss of charged particles in dilute plasmas. Physical Review A, 29(4):2145, 1984. [60] Simon W de Leeuw, John William Perram, and Edgar Roderick Smith. Simulation of electrostatic In Proceedings systems in periodic boundary conditions. i. lattice sums and dielectric constants. of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, volume 373, pages 27–56. The Royal Society, 1980. [61] Markus Deserno and Christian Holm. How to mesh up ewald sums. ii. an accurate error estimate for the particle–particle–particle-mesh algorithm. The Journal of chemical physics, 109(18):7694–7701, 1998. [62] AW DeSilva and JD Katsouros. Electrical conductivity of dense copper and aluminum plasmas. Physical Review E, 57(5):5945, 1998. [63] MWC Dharma-Wardana. Electron-ion and ion-ion potentials for modeling warm dense matter: Ap- plications to laser-heated or shock-compressed al and si. Physical Review E, 86(3):036407, 2012. 96 [64] MWC Dharma-Wardana and Michael S Murillo. Pair-distribution functions of two-temperature two- mass systems: Comparison of molecular dynamics, classical-map hypernetted chain, quantum monte carlo, and kohn-sham calculations for dense hydrogen. Physical Review E, 77(2):026401, 2008. [65] Abdourahmane Diaw and Michael Sean Murillo. Generalized hydrodynamics model for strongly coupled plasmas. Physical Review E, 92(1):013107, 2015. [66] Hong-Qiang Ding, Naoki Karasawa, and William A Goddard. The reduced cell multipole method for coulomb interactions in periodic systems with million-atom unit cells. Chemical Physics Letters, 196(1-2):6–10, 1992. [67] Paul AM Dirac. Note on exchange phenomena in the thomas atom. In Mathematical Proceedings of the Cambridge Philosophical Society, volume 26, pages 376–385. Cambridge University Press, 1930. [68] TP Doerr and Yi-Kuo Yu. Electrical interactions in the cell: Asymmetric screening in a watery “antiverse”. American journal of physics, 82(5):460–465, 2014. [69] Daniel HE Dubin. Measurement of screening enhancement to nuclear reaction rates using a strongly magnetized and strongly correlated non-neutral plasma. Physical review letters, 94(2):025002, 2005. [70] G Ecker and W Kröll. Lowering of the ionization energy for a plasma in thermodynamic equilibrium. The Physics of Fluids, 6(1):62–69, 1963. [71] Wolfgang Eckhardt, Alexander Heinecke, Reinhold Bader, Matthias Brehm, Nicolay Hammer, Her- bert Huber, Hans-Georg Kleinhenz, Jadran Vrabec, Hans Hasse, Martin Horsch, et al. 591 tflops multi-trillion particles simulation on supermuc. In International Supercomputing Conference, pages 1–12. Springer, 2013. [72] Shalom Eliezer. Fundamentals of equations of state. Allied Publishers, 2005. [73] Ulrich Essmann, Lalith Perera, Max L Berkowitz, Tom Darden, Hsing Lee, and Lee G Pedersen. A smooth particle mesh ewald method. The Journal of chemical physics, 103(19):8577–8593, 1995. [74] Paul P Ewald. Die berechnung optischer und elektrostatischer gitterpotentiale. Annalen der Physik, 369(3):253–287, 1921. [75] Hans Feldmeier and Jürgen Schnack. Molecular dynamics for fermions. Reviews of Modern Physics, 72(3):655, 2000. [76] Nicolas Ferré, Michael Filatov, and Miquel Huix-Rotllant. Density-functional methods for excited states, volume 368. Springer, 2015. [77] Daniel S Fisher, BI Halperin, and PM Platzman. Phonon-ripplon coupling and the two-dimensional electron solid on a liquid-helium surface. Physical Review Letters, 42(12):798, 1979. [78] Alex Fletcher, Sigrid Close, and Donovan Mathias. Simulating plasma production from hypervelocity impacts. Physics of Plasmas, 22(9):093504, 2015. [79] LB Fletcher, AL Kritcher, A Pak, T Ma, T Döppner, C Fortmann, L Divol, OS Jones, OL Landen, HA Scott, et al. Observations of continuum depression in warm dense matter with x-ray thomson scattering. Physical review letters, 112(14):145004, 2014. 97 [80] LB Fletcher, HJ Lee, T Döppner, E Galtier, B Nagler, P Heimann, C Fortmann, S LePape, T Ma, M Millot, et al. Ultrabright x-ray laser scattering for dynamic warm dense matter physics. Nature Photonics, 9(4):274–279, 2015. [81] Roger Fletcher. Practical methods of optimization. John Wiley & Sons, 2013. [82] [83] Jonathan J Fortney and Nadine Nettelmann. The interior structure, composition, and evolution ofá- giant planets. Space Science Reviews, 152(1-4):423–447, 2010. JA Frenje, PE Grabowski, CK Li, FH Séguin, AB Zylstra, M Gatu Johnson, RD Petrasso, V Yu Glebov, and TC Sangster. Measurements of ion stopping around the bragg peak in high-energy- density plasmas. Physical review letters, 115(20):205001, 2015. [84] Daan Frenkel and Berend Smit. Understanding molecular simulation: from algorithms to applica- tions, volume 1. Elsevier, 2001. [85] [86] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. The elements of statistical learning, vol- ume 1. Springer series in statistics New York, 2001. Ikuo Fukuda and Haruki Nakamura. Non-ewald methods: systems. Biophysical reviews, 4(3):161–170, 2012. theory and applications to molecular [87] Camelia Gabriel, Sami Gabriel, and E Corthout. The dielectric properties of biological tissues: I. literature survey. Physics in medicine and biology, 41(11):2231, 1996. [88] EJ Gamboa, CM Huntington, MR Trantham, PA Keiter, RP Drake, DS Montgomery, John F Benage, and Samuel A Letzring. Imaging x-ray thomson scattering spectrometer design and demonstration (invited) a. Review of Scientific Instruments, 83(10):10E108, 2012. [89] EJ Gamboa, DS Montgomery, IM Hall, and RP Drake. Imaging x-ray crystal spectrometer for laser- produced plasmas. Journal of Instrumentation, 6(04):P04004, 2011. [90] DO Gericke, MS Murillo, and M Schlanges. Dense plasma temperature equilibration in the binary collision approximation. Physical Review E, 65(3):036418, 2002. [91] SH Glenzer, LB Fletcher, E Galtier, B Nagler, R Alonso-Mori, B Barbrel, SB Brown, DA Chapman, Z Chen, CB Curry, et al. Matter under extreme conditions experiments at the linac coherent light source. Journal of Physics B: Atomic, Molecular and Optical Physics, 49(9):092001, 2016. [92] JN Glosli, FR Graziani, RM More, MS Murillo, FH Streitz, MP Surh, LX Benedict, S Hau-Riege, AB Langdon, and RA London. Molecular dynamics simulations of temperature equilibration in dense hydrogen. Physical Review E, 78(2):025401, 2008. [93] Paul E Grabowski, Andreas Markmann, Igor V Morozov, Ilya A Valuev, Christopher A Fichtl, David F Richards, Victor S Batista, Frank R Graziani, and Michael S Murillo. Wave packet spreading and localization in electron-nuclear scattering. Physical Review E, 87(6):063104, 2013. [94] Paul E Grabowski, Michael P Surh, David F Richards, Frank R Graziani, and Michael S Murillo. Molecular dynamics simulations of classical stopping power. Physical review letters, 111(21):215002, 2013. 98 [95] CC Grimes and G Adams. Evidence for a liquid-to-crystal phase transition in a classical, two- dimensional sheet of electrons. Physical Review Letters, 42(12):795, 1979. [96] A Grinenko, V Tz Gurovich, A Saypin, S Efimov, Ya E Krasik, and Vladimir Ivanovich Oreshkin. Strongly coupled copper plasma generated by underwater electrical wire explosion. Physical Review E, 72(6):066401, 2005. [97] P Gupta, S Laha, CE Simien, H Gao, J Castro, TC Killian, and T Pohl. Electron-temperature evolution in expanding ultracold neutral plasmas. Physical review letters, 99(7):075005, 2007. [98] JP Hansen and IR McDonald. Microscopic simulation of a hydrogen plasma. Physical Review Letters, 41(20):1379, 1978. [99] Zach Haralson and J Goree. Overestimation of viscosity by the green-kubo method in a dusty plasma experiment. Physical Review Letters, 118(19):195001, 2017. [100] N Heilmann, JB Peatross, and Scott D Bergeson. “ultracold” neutral plasmas at room temperature. Physical review letters, 109(3):035002, 2012. [101] Roger W Hockney and James W Eastwood. Computer simulation using particles. crc Press, 1988. [102] SX Hu, B Militzer, VN Goncharov, and S Skupsky. Strong coupling and degeneracy effects in inertial confinement fusion implosions. Physical review letters, 104(23):235003, 2010. [103] Rupert Huber, F Tauser, A Brodschelm, M Bichler, G Abstreiter, and A Leitenstorfer. How many-particle interactions develop after ultrafast excitation of an electron-hole plasma. Nature, 414(6861):286–290, 2001. [104] J Hughto, AS Schneider, CJ Horowitz, and DK Berry. Diffusion of neon in white dwarf stars. Physical Review E, 82(6):066401, 2010. [105] Setsuo Ichimaru. Strongly coupled plasmas: high-density classical plasmas and degenerate electron liquids. Reviews of Modern Physics, 54(4):1017, 1982. [106] Rolf E Isele-Holder, Wayne Mitchell, and Ahmed E Ismail. Development and application of a particle-particle particle-mesh ewald method for dispersion interactions. The Journal of chemical physics, 137(17):174107, 2012. [107] MJ Jensen, T Hasegawa, John J Bollinger, and DHE Dubin. Rapid heating of a strongly coupled plasma near the solid-liquid phase transition. Physical review letters, 94(2):025001, 2005. [108] Wei Jin, Priya Vashishta, Rajiv K Kalia, and José P Rino. Dynamic structure factor and vibrational properties of sio 2 glass. Physical Review B, 48(13):9359, 1993. [109] Christopher S Jones and Michael S Murillo. Analysis of semi-classical potentials for molecular dynamics and monte carlo simulations of warm dense matter. High energy density physics, 3(3):379– 394, 2007. [110] Kai Kadau, Timothy C Germann, Nicolas G Hadjiconstantinou, Peter S Lomdahl, Guy Dimonte, Brad Lee Holian, and Berni J Alder. Nanohydrodynamics simulations: an atomistic view of the rayleigh–taylor instability. Proceedings of the National Academy of Sciences of the United States of America, 101(16):5851–5855, 2004. 99 [111] Sergey A Khrapak, Boris Klumov, Lénaïc Couëdel, and Hubertus M Thomas. On the long-waves dispersion in yukawa systems. Physics of Plasmas, 23(2):023702, 2016. [112] TC Killian, S Kulin, SD Bergeson, Luis A Orozco, C Orzel, and SL Rolston. Creation of an ultracold neutral plasma. Physical Review Letters, 83(23):4776, 1999. [113] Thomas C Killian. Ultracold neutral plasmas. Science, 316(5825):705–708, 2007. [114] John G Kirkwood. Statistical mechanics of fluid mixtures. The Journal of Chemical Physics, 3(5):300–313, 1935. [115] CL Kirschbaum and L Wilets. Classical many-body model for atomic collisions incorporating the heisenberg and pauli principles. Physical Review A, 21(3):834, 1980. [116] John L Klepeis, Kresten Lindorff-Larsen, Ron O Dror, and David E Shaw. Long-timescale molec- ular dynamics simulations of protein structure and function. Current opinion in structural biology, 19(2):120–127, 2009. [117] MD Knudson, MP Desjarlais, and DH Dolan. Shock-wave exploration of the high-pressure phases of carbon. Science, 322(5909):1822–1825, 2008. [118] Jiri Kolafa and John W Perram. Cutoff errors in the ewald summation formulae for point charge systems. Molecular Simulation, 9(5):351–368, 1992. [119] Takahiro Komatsu, Noriyuki Yoshii, Shinichi Miura, and Susumu Okazaki. A large-scale molecular dynamics study of dynamic structure factor and dispersion relation of acoustic mode in liquid and supercritical water. Fluid phase equilibria, 226:345–350, 2004. [120] PM Kozlowski, BJB Crowley, Dirk O Gericke, Sean P Regan, and Gianluca Gregori. Theory of thomson scattering in inhomogeneous media. Scientific reports, 6:24283, 2016. [121] A Kramida, Yu Ralchenko, J Reader, et al. Nist atomic spectra database (ver. 5.2), national institute of standards and technology, gaithersburg, md, 2014. There is no corresponding record for this reference, 2014. [122] D Kraus, DA Chapman, AL Kritcher, RA Baggott, B Bachmann, GW Collins, SH Glenzer, JA Hawre- liak, DH Kalantar, OL Landen, et al. X-ray scattering measurements on imploding ch spheres at the national ignition facility. Physical Review E, 94(1):011202, 2016. [123] Georg Kresse and Jürgen Hafner. Ab initio molecular dynamics for liquid metals. Physical Review B, 47(1):558, 1993. [124] Konstantin N Kudin and Gustavo E Scuseria. A fast multipole method for periodic systems with arbitrary unit cell geometries. Chemical Physics Letters, 283(1):61–68, 1998. [125] S Kulin, TC Killian, SD Bergeson, and SL Rolston. Plasma oscillations and expansion of an ultracold neutral plasma. Physical review letters, 85(2):318, 2000. [126] F Lado. Effective potential description of the quantum ideal gases. The Journal of Chemical Physics, 47(12):5369–5375, 1967. 100 [127] KJ LaGattuta. Multiple ionization of helium clusters by long wavelength laser radiation. The Euro- pean Physical Journal D-Atomic, Molecular, Optical and Plasma Physics, 2(3):267–272, 1998. [128] OL Landen and RJ Winfield. Laser scattering from dense cesium plasmas. Physical review letters, 54(15):1660, 1985. [129] Peter B Lerner, Kenneth J LaGattuta, and James S Cohen. Ionization of helium by a short pulse of radiation: A fermi molecular-dynamics calculation. Physical Review A, 49(1):R12, 1994. [130] Richard L Liboff. Kinetic theory: classical, quantum, and relativistic descriptions. Springer Science & Business Media, 2003. [131] John Lindl. Development of the indirect-drive approach to inertial confinement fusion and the target physics basis for ignition and gain. Physics of plasmas, 2(11):3933–4024, 1995. [132] M Lyon, Scott D Bergeson, and MS Murillo. Limit of strong ion coupling due to electron shielding. Physical Review E, 87(3):033101, 2013. [133] Viet Hoang Man, Nguyen-Thi Van-Oanh, Philippe Derreumaux, Mai Suan Li, Christopher Roland, Celeste Sagui, and Phuong H Nguyen. Picosecond infrared laser-induced all-atom nonequilibrium molecular dynamics simulation of dissociation of viruses. Physical Chemistry Chemical Physics, 18(17):11951–11958, 2016. [134] Paul Maragakis, Kresten Lindorff-Larsen, Michael P Eastwood, Ron O Dror, John L Klepeis, Isaiah T Arkin, Morten Ø Jensen, Huafeng Xu, Nikola Trbovic, Richard A Friesner, et al. Microsecond molec- ular dynamics simulation shows effect of slow loop dynamics on backbone amide order parameters of proteins. The journal of physical chemistry. B, 112(19):6155, 2008. [135] Miguel AL Marques, Neepa T Maitra, Fernando MS Nogueira, Eberhard KU Gross, and Angel Ru- bio. Fundamentals of time-dependent density functional theory, volume 837. Springer Science & Business Media, 2012. [136] JA Maruhn, P-G Reinhard, PD Stevenson, and AS Umar. The tdhf code sky3d. Computer Physics Communications, 185(7):2195–2216, 2014. [137] Martial Mazars. Long ranged interactions in computer simulations and for quasi-2d systems. Physics Reports, 500(2):43–116, 2011. [138] EJ McGuire, JM Peek, and LC Pitchford. Proton stopping power of aluminum ions. Physical Review A, 26(3):1318, 1982. [139] A McKelvey, GE Kemp, PA Sterne, A Fernandez-Panella, R Shepherd, M Marinak, A Link, GW Collins, H Sio, J King, et al. Thermal conductivity measurements of proton-heated warm dense aluminum. Scientific Reports, 7, 2017. [140] Matthias Mecke, Jochen Winkelmann, and Johann Fischer. Molecular dynamics simulation of the liquid–vapor interface: The lennard-jones fluid. The Journal of chemical physics, 107(21):9264– 9270, 1997. [141] TB Mitchell, JJ Bollinger, DHE Dubin, X-P Huang, WM Itano, and RH Baughman. Direct obser- vations of structural phase transitions in planar crystallized ion plasmas. Science, 282(5392):1290– 1293, 1998. 101 [142] Normand Arthur Modine and Ryan M Hatcher. Representing the thermal state in time-dependent density functional theory. The Journal of chemical physics, 142(20):204111, 2015. [143] S Morita, N Matsuda, N Toshima, and K Hino. Ionization of stabilized helium atoms by proton and antiproton impacts. Physical Review A, 66(4):042719, 2002. [144] EI Moses, RN Boyd, BA Remington, CJ Keane, and R Al-Ayat. The national ignition facility: Ushering in a new age for high energy density science. Physics of Plasmas, 16(4):041006, 2009. [145] Michael S Murillo. Strongly coupled plasma physics and high energy-density matter. Physics of Plasmas, 11(5):2964–2971, 2004. [146] Michael S Murillo. Ultrafast dynamics of strongly coupled plasmas. Physical review letters, 96(16):165001, 2006. [147] Michael S Murillo. Viscosity estimates of liquid metals and warm dense matter using the yukawa reference system. High Energy Density Physics, 4(1):49–57, 2008. [148] Michael S Murillo. X-ray thomson scattering in warm dense matter at low frequencies. Physical Review E, 81(3):036403, 2010. [149] Michael S Murillo and MWC Dharma-Wardana. Temperature relaxation in hot dense hydrogen. Physical review letters, 100(20):205005, 2008. [150] MS Murillo. Using fermi statistics to create strongly coupled ion plasmas in atom traps. Physical review letters, 87(11):115003, 2001. [151] D Murphy and BM Sparkes. Disorder-induced heating of ultracold neutral plasmas created from atoms in partially filled optical lattices. Physical Review E, 94(2):021201, 2016. [152] Bob Nagler, Brice Arnold, Gary Bouchard, Richard F Boyce, Richard M Boyce, Alice Callen, Marc Campell, Ruben Curiel, Eric Galtier, Justin Garofoli, et al. The matter in extreme conditions instru- ment at the linac coherent light source. Journal of synchrotron radiation, 22(3):520–525, 2015. [153] Alexey Neelov and Christian Holm. Interlaced p3m algorithm with analytical and ik-differentiation. The Journal of chemical physics, 132(23):234103, 2010. [154] Franziska Nestler. Parameter tuning for the nfft based fast ewald summation. Frontiers in Physics, 4:28, 2016. [155] V Nosenko and J Goree. Shear flows and shear viscosity in a two-dimensional yukawa system (dusty plasma). Physical review letters, 93(15):155004, 2004. [156] V Nosenko, S Zhdanov, AV Ivlev, G Morfill, J Goree, and A Piel. Heat transport in a two-dimensional complex (dusty) plasma at melting conditions. Physical review letters, 100(2):025003, 2008. [157] P Nouvel, H Marinchio, J Torres, C Palermo, D Gasquet, L Chusseau, L Varani, P Shiktorov, E Starikov, and V Gružinskis. Terahertz spectroscopy of plasma waves in high electron mobility transistors. Journal of Applied Physics, 106(1):013717, 2009. [158] Thomas M Nymand and Per Linse. Ewald summation and reaction field methods for potentials with atomic charges, dipoles, and polarizabilities. The Journal of Chemical Physics, 112(14):6152–6160, 2000. 102 [159] Jae-Yong Oh, Yang-Hyun Koo, Jin-Sik Cheon, Byung-Ho Lee, and Dong-Seong Sohn. Molecular dynamics simulation of the pressure–volume–temperature data of xenon for a nuclear fuel. Journal of nuclear materials, 372(1):89–93, 2008. [160] Ari Ojanperä, Arkady V Krasheninnikov, and Martti Puska. Electronic stopping power from first- principles calculations with account for core electron excitations and projectile ionization. Physical Review B, 89(3):035120, 2014. [161] T Ott, M Bonitz, LG Stanton, and MS Murillo. Coupling strength in coulomb and yukawa one- component plasmas. Physics of Plasmas, 21(11):113704, 2014. [162] Fang Peng, Edward F Brown, and James W Truran. Sedimentation and type i x-ray bursts at low accretion rates. The Astrophysical Journal, 654(2):1022, 2007. [163] Juan R Perilla, Jodi A Hadden, Boon Chong Goh, Christopher G Mayne, and Klaus Schulten. All- atom molecular dynamics of virus capsids as drug targets. The journal of physical chemistry letters, 7(10):1836, 2016. [164] John W Perram, Henrik G Petersen, and Simon W De Leeuw. An algorithm for the simulation of condensed matter which grows as the 3/2 power of the number of particles. Molecular Physics, 65(4):875–893, 1988. [165] JB Pieper and J Goree. Dispersion of plasma dust acoustic waves in the strong-coupling regime. Physical review letters, 77(15):3137, 1996. [166] Bill Poirier. Bohmian mechanics without pilot waves. Chemical Physics, 370(1):4–14, 2010. [167] EL Pollock and Jim Glosli. Comments on p3m, fmm, and the ewald method for large periodic coulombic systems. Computer Physics Communications, 95(2-3):93–110, 1996. [168] Alexander Y Potekhin, Gilles Chabrier, and Forrest J Rogers. Equation of state of classical coulomb plasma mixtures. Physical Review E, 79(1):016411, 2009. [169] Aurora Pribram-Jones, Paul E Grabowski, and Kieron Burke. Thermal density functional theory: Time-dependent linear response and approximate functionals from the fluctuation-dissipation theo- rem. Physical review letters, 116(23):233001, 2016. [170] Ferry PRINS, Aaron J GOODMAN, and William A TISDALE. Reduced dielectric screening and enhanced energy transfer in single-and few-layer mos2. Nano letters, 14(11):6087–6091, 2014. [171] Yue Qi, Tahir Cagin, Yoshitaka Kimura, and William A Goddard III. Molecular-dynamics simulations of glass formation and crystallization in binary liquid metals: Cu-ag and cu-ni. Physical review B, 59(5):3527, 1999. [172] Edwin E Quashie, Bidhan C Saha, and Alfredo A Correa. Electronic band structure effects in the stopping of protons in copper. Physical Review B, 94(15):155403, 2016. [173] DS Rackstraw, O Ciricosta, SM Vinko, B Barbrel, T Burian, J Chalupsk`y, BI Cho, H-K Chung, GL Dakovski, K Engelhorn, et al. Saturable absorption of an x-ray free-electron-laser heated solid- density aluminum plasma. Physical review letters, 114(1):015003, 2015. 103 [174] J Randrup, C Dorso, and S Duarte. Quasi-classical treatment of the nucleon gas. Le Journal de Physique Colloques, 48(C2):C2–143, 1987. [175] Dennis C Rapaport. The art of molecular dynamics simulation. Cambridge university press, 2004. [176] F Remacle and RD Levine. On the classical limit for electronic structure and dynamics in the orbital approximation. The Journal of Chemical Physics, 113(11):4515–4523, 2000. [177] Yaakov Rosenfeld. Ewald method for simulating yukawa systems. Molecular Physics, 88(5):1357– 1363, 1996. [178] M Ross and AK McMahan. Condensed xenon at high pressure. Physical Review B, 21(4):1658, 1980. [179] Carl A Rouse. Ionization equilibrium equation of state. The Astrophysical Journal, 134:435, 1961. [180] Balazs F Rozsnyai. Electron scattering in hot/warm plasmas. High Energy Density Physics, 4(1):64– 72, 2008. [181] Erich Runge and Eberhard KU Gross. Density-functional theory for time-dependent systems. Physi- cal Review Letters, 52(12):997, 1984. [182] Megh Nad Saha. Liii. ionization in the solar chromosphere. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 40(238):472–488, 1920. [183] Megh Nad Saha. On a physical theory of stellar spectra. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 99(697):135–153, 1921. [184] T Saigo and S Hamaguchi. Shear viscosity of strongly coupled yukawa systems. Physics of plasmas, 9(4):1210–1216, 2002. [185] Gwenaël Salin and Jean-Michel Caillol. Ewald sums for yukawa potentials. The Journal of Chemical Physics, 113(23):10459–10463, 2000. [186] Edwin E Salpeter. Electron density fluctuations in a plasma. Physical Review, 120(5):1528, 1960. [187] Johannes Sarnthein, Alfredo Pasquarello, and Roberto Car. Origin of the high-frequency doublet in the vibrational spectrum of vitreous sio2. Science, 275(5308):1925–1927, 1997. [188] Jeremy Schiff and Bill Poirier. Communication: Quantum mechanics without wavefunctions, 2012. [189] Tamar Schlick and Charles S Peskin. Can classical equations simulate quantum-mechanical behavior? a molecular dynamics investigation of a diatomic molecule with a morse potential. Communications on Pure and Applied Mathematics, 42(8):1141–1163, 1989. [190] P Seuferling, J Vogel, and C Toepffer. Correlations in a two-temperature plasma. Physical Review A, 40(1):323, 1989. [191] Dmitrii V Shalashilin and Irene Burghardt. Gaussian-based techniques for quantum propagation from the time-dependent variational principle: Formulation in terms of trajectories of coupled classical and quantum variables. The Journal of chemical physics, 129(8):084104, 2008. 104 [192] David E Shaw, Ron O Dror, John K Salmon, JP Grossman, Kenneth M Mackenzie, Joseph A Bank, Cliff Young, Martin M Deneroff, Brannon Batson, Kevin J Bowers, et al. Millisecond-scale molecular dynamics simulations on anton. In High performance computing networking, storage and analysis, proceedings of the conference on, pages 1–11. IEEE, 2009. [193] Vincent K Shen, Raymond D Mountain, and Jeffrey R Errington. Comparative study of the effect of tail corrections on surface tension determined by molecular simulation. The Journal of Physical Chemistry B, 111(22):6198–6207, 2007. [194] Jiro Shimada, Hiroki Kaneko, and Toshikazu Takada. Efficient calculations of coulombic interactions in biomolecular simulations with periodic boundary conditions. Journal of Computational Chemistry, 14(7):867–878, 1993. [195] CE Simien, YC Chen, P Gupta, S Laha, YN Martinez, PG Mickelson, SB Nagel, and TC Killian. Using absorption imaging to study ion dynamics in an ultracold neutral plasma. Physical review letters, 92(14):143001, 2004. [196] Jeffrey S Simonoff. Smoothing methods in statistics. Springer Science & Business Media, 2012. [197] John C Snyder, Matthias Rupp, Katja Hansen, Klaus-Robert Müller, and Kieron Burke. Finding density functionals with machine learning. Physical review letters, 108(25):253002, 2012. [198] Philipp Sperling, EJ Gamboa, HJ Lee, HK Chung, E Galtier, Y Omarbakiyeva, H Reinholz, G Röpke, U Zastrau, J Hastings, et al. Free-electron x-ray laser measurements of collisional-damped plasmons in isochorically heated warm dense matter. Physical review letters, 115(11):115001, 2015. [199] LG Stanton and MS Murillo. Unified description of linear screening in dense plasmas. Physical Review E, 91(3):033104, 2015. [200] Liam G Stanton and Michael S Murillo. Ionic transport in high-energy-density matter. Physical Review E, 93(4):043203, 2016. [201] Harry A Stern and Keith G Calkins. On mesh-based ewald methods: Optimal parameters for two differentiation schemes. The Journal of chemical physics, 128(21):214106, 2008. [202] John C Stewart and Kedar D Pyatt Jr. Lowering of ionization potentials in plasmas. The Astrophysical Journal, 144:1203, 1966. [203] Trevor S Strickler, Thomas K Langin, Paul McQuillen, J Daligault, and Tom C Killian. Experimental measurement of self-diffusion in a strongly coupled plasma. Physical Review X, 6(2):021021, 2016. [204] Julius T Su and William A Goddard III. Excited electron dynamics modeling of warm dense matter. Physical review letters, 99(18):185003, 2007. [205] AP Sutton and J Chen. Long-range finnis–sinclair potentials. Philosophical Magazine Letters, 61(3):139–146, 1990. [206] AGR Thomas, M Tzoufras, APL Robinson, RJ Kingham, CP Ridgers, Mark Sherlock, and AR Bell. A review of vlasov–fokker–planck numerical modeling of inertial confinement fusion plasma. Jour- nal of Computational Physics, 231(3):1051–1079, 2012. 105 [207] H Thomas, GE Morfill, V Demmel, J Goree, B Feuerbacher, and D Möhlmann. Plasma crystal: Coulomb crystallization in a dusty plasma. Physical Review Letters, 73(5):652, 1994. [208] K Tiedtke, A Azima, N Von Bargen, L Bittner, S Bonfigt, S Düsterer, B Faatz, U Frühling, M Gensch, Ch Gerth, et al. The soft x-ray free-electron laser flash at desy: beamlines, diagnostics and end- stations. New journal of physics, 11(2):023029, 2009. [209] H Totsuji. Conduction of strongly coupled electrons through narrow channels on liquid helium sur- face: Simulation. Contributions to Plasma Physics, 52(1):74–77, 2012. [210] Abdulnour Y Toukmaji and John A Board. Ewald summation techniques in perspective: a survey. Computer physics communications, 95(2-3):73–92, 1996. [211] Ya-Yi Tsai, Jun-Yi Tsai, and I Lin. Generation of acoustic rogue waves in dusty plasmas through three-dimensional particle focusing by distorted waveforms. Nature Physics, 12(6):573–577, 2016. [212] Mark Tuckerman. Statistical mechanics: theory and molecular simulation. Oxford University Press, 2010. [213] John P Verboncoeur. Particle simulation of plasmas: review and advances. Plasma Physics and Controlled Fusion, 47(5A):A231, 2005. [214] SM Vinko, O Ciricosta, BI Cho, K Engelhorn, H-K Chung, CRD Brown, T Burian, J Chalupsk`y, RW Falcone, C Graves, et al. Creation and diagnosis of a solid-density plasma with an x-ray free- electron laser. Nature, 482(7383):59–62, 2012. [215] SM Vinko, O Ciricosta, TR Preston, DS Rackstraw, CRD Brown, T Burian, J Chalupsk`y, Investigation of femtosecond collisional ion- Byeoung Ick Cho, H-K Chung, K Engelhorn, et al. ization rates in a solid-density aluminium plasma. Nature communications, 6:ncomms7397, 2015. [216] SM Vinko, O Ciricosta, and JS Wark. Density functional theory calculations of continuum lowering in strongly coupled plasmas. Nature communications, 5:ncomms4533, 2014. [217] DA Wasson and SE Koonin. Molecular-dynamics simulations of atomic ionization by strong laser fields. Physical Review A, 39(11):5676, 1989. [218] Ulrich Welling and Guido Germano. Efficiency of linked cell algorithms. Computer Physics Com- munications, 182(3):611–615, 2011. [219] TG White, S Richardson, BJB Crowley, LK Pattison, JWO Harris, and G Gregori. Orbital-free density-functional theory simulations of the dynamic structure factor of warm dense aluminum. Phys- ical review letters, 111(17):175002, 2013. [220] L Wilets, EM Henley, M Kraft, and AD MacKellar. Classical many-body model for heavy-ion colli- sions incorporating the pauli principle. Nuclear Physics A, 282(2):341–350, 1977. [221] Bastian BL Witte, LB Fletcher, E Galtier, E Gamboa, HJ Lee, U Zastrau, R Redmer, SH Glenzer, and P Sperling. Warm dense matter demonstrating non-drude conductivity from observations of nonlinear plasmon damping. Physical Review Letters, 118(22):225001, 2017. 106 [222] D Wolf, P Keblinski, SR Phillpot, and J Eggebrecht. Exact method for the simulation of coulom- bic systems by spherically truncated, pairwise r- 1 summation. The Journal of chemical physics, 110(17):8254–8282, 1999. [223] Mofreh R Zaghloul, Mohamed A Bourham, and J Michael Doster. A simple formulation and solution strategy of the saha equation for ideal and nonideal plasmas. Journal of Physics D: Applied Physics, 33(8):977, 2000. [224] Ming-Liang Zhang and DA Drabold. An extension of the kubo-greenwood formula for use in molec- ular simulations. arXiv preprint arXiv:0904.0212, 2009. [225] W Zhang, Abul K Azad, and D Grischkowsky. Terahertz studies of carrier dynamics and dielectric response of n-type, freestanding epitaxial gan. Applied Physics Letters, 82(17):2841–2843, 2003. [226] J Zheng, QF Chen, YJ Gu, ZY Chen, and CJ Li. Thermodynamics, compressibility, and phase The Journal of chemical physics, diagram: Shock compression of supercritical fluid xenon. 141(12):124201, 2014. [227] Yueming Zhou, Cheng Huang, Qing Liao, and Peixiang Lu. Classical simulations including electron correlations for sequential double ionization. Physical review letters, 109(5):053004, 2012. [228] AB Zylstra, JA Frenje, PE Grabowski, CK Li, GW Collins, P Fitzsimmons, S Glenzer, F Graziani, SB Hansen, SX Hu, et al. Development of a wdm platform for charged-particle stopping experiments. In Journal of Physics: Conference Series, volume 717, page 012118. IOP Publishing, 2016. 107