Molecular Dynamics Simulations of Temperature Equilibration in Dense Hydrogen
J. N. Glosli, F. R. Graziani, R. M. More, M. S. Murillo, F. H. Streitz, M. P. Surh, L. X. Benedict, S. Hau-Riege, A. B. Langdon, R. A. London
MMolecular Dynamics Simulations of Temperature Equilibration in Dense Hydrogen
J. N. Glosli, F. R. Graziani, R. M. More, M. S. Murillo, F. H. Streitz, M. P. Surh, L. X. Benedict, S. Hau-Riege, A. B. Langdon, and R. A. London Lawrence Livermore National Laboratory Physics Division, MS D410, Los Alamos National Laboratory, Los Alamos, NM 87545 (Dated: October 26, 2018; Received)The temperature equilibration rate in dense hydrogen (for both T i > T e and T i < T e ) hasbeen calculated with molecular dynamics simulations for temperatures between 10 and 600 eV anddensities between 10 /cc to 10 /cc . Careful attention has been devoted to convergence of thesimulations, including the role of semiclassical potentials. We find that for Coulomb logarithms L (cid:38)
1, a model by Gericke-Murillo-Schlanges (GMS) [Gericke et al., PRE , 036418 (2002)] basedon a T-matrix method and the approach by Brown-Preston-Singleton [Brown et al., Phys. Rep. ,237 (2005)] agrees with the simulation data to within the error bars of the simulation. For smallerCoulomb logarithms, the GMS model is consistent with the simulation results. Landau-Spitzermodels are consistent with the simulation data for L > I. INTRODUCTION
The strong temperature dependence of thermonuclearreaction rates suggests that even small deviations fromequilibrium can yield differences in burn rates. Thus, thepursuit of ignition in the laboratory will benefit from ac-curate models of relaxation processes in hot, dense plas-mas. One of the greatest uncertainties in the nonequi-librium energy balance is the electron-ion temperaturerelaxation rate. Although there have been indirect mea-surements for cool dense matter [1], there is no experi-mental data in the regime of interest. Even worse, theo-retical descriptions of Coulomb collisions suffer from di-vergences that make detailed models difficult to develop.Here we take a complementary approach to hot, denseplasmas by using molecular dynamics (MD) techniques.We use this method to test recent theoretical models andcompare with standard results.The electron-proton coupling rate was first calculatedby Landau [2] and Spitzer (LS) [3] for classical plasmaswith weak collisions. They write the electron-proton tem-perature exchange rate (1 /τ pe ) in the form,1 τ pe = 8 √ πn p Z e m e m p c (cid:26) k B T e m e c + k B T p m p c (cid:27) − / L ≡ L J LS , (1)where J LS is the LS pre-factor, n e ( n p ) are the electron(ion) number densities, Z = 1 is the proton charge, T e ( T p ) are the electron (ion) temperatures, and k B is theBoltzmann constant. L is the so-called Coulomb loga-rithm containing details of the collision process. LS used L LS = ln( b max /b min ) (2)where b max and b min are impact parameter cutoffsneeded to remove divergences that arose from their treat-ment. b min is chosen to be a minimum impact pa-rameter consistent with plasma conditions, such as theclassical distance of closest approach ( b C = Ze /k B T ).At high temperatures, b min is often modified to includequantum diffraction effects by introducing the length scale of the electron thermal deBroglie wavelength Λ = (cid:112) π (cid:125) /m e k B T e . Typically b max is chosen to be a screen-ing length arising from collective plasma phenomena,such as the Debye length λ D = (cid:112) k B T e / πe n e .The presence of ad hoc cut-offs and other inconsisten-cies led researchers to derive kinetic equations withoutcut-offs [5, 6, 7, 8]. The essence of these theories is the in-clusion of strong scattering in the presence of dynamicalcollective (screening) behavior. Two such theories havebeen recently proposed: Gericke, Murillo and Schlanges(GMS) [8] and Brown, Preston and Singleton (BPS) [9].GMS applied these ideas to dense plasma temperatureequilibration, They investigated various approximationsin evaluating of L , including issues with trajectories andcutoffs, and provided four different evaluations of the re-laxation rate based on quantum kinetic theory. Fromtheir numerical work, GMS suggest an effective Coulomblogarithm [10] L GMS = 12 ln (cid:0) (cid:2) λ D + R ion (cid:3) / (cid:2) Λ / π + b C (cid:3)(cid:1) , (3)where R ion = (3 / πn p ) / is the ion sphere radius. Thisexpression was described by GMS as the best fit to theirfull T-matrix theory.BPS and Brown and Singleton [11] used dimensionalcontinuation to obtain an expression for the electron-ioncoupling rate accurate to second order in the plasmacoupling parameter. The method is applicable to bothdegenerate and non-degenerate electrons. For the non-degenerate case, they derive L BP S = log( λ D / Λ) + (log (16 π ) − γ − / , (4)where γ is the Euler constant.The most direct method of studying temperature equi-libration in the classical limit is with numerical simula-tion; strong, collective scattering at all length scales isthe forte of MD. Hansen and McDonald (HM) [12] ex-plored temperature equilibration in dense hydrogen us-ing MD, comparing their results against a LS model with a r X i v : . [ phy s i c s . p l a s m - ph ] F e b L LS = ln(2 πλ D / Λ). However, the HM simulations in-volved a very small number of particles ( N = 128) withpresumably large error bars. Here, we expand upon theircalculations to not only reassess the HM result, but alsocompare with the modern approaches of GMS and BPS. II. MOLECULAR DYNAMICS: SIMULATIONSAND RESULTS
MD simulations are applied to two-temperature sys-tems of charged particles in a cubic cell with periodicboundary conditions. The MD is performed with a fullyparallel code using a basic leapfrog method[13] with theCoulomb interaction evaluated by an Ewald summation[14, 15]. Because the classical Coulomb many-body prob-lem is unstable for attractive interactions, we employsemi-classical potentials that reduce the Coulomb inter-action on short length scales in order to prevent unphys-ical, deeply-bound states. We tested several forms ofthe diffractive [17, 18] and Pauli [19, 20] terms for thesepotentials. The resulting equilibration times typicallyvary by less than 15%, which is within the statistical er-ror of the MD data. The similarity is not unexpected,since most semi-classical potentials resemble one anotherabove 10 eV [22]. We report results using the semi-classical potential in HM [12], V ab ( r ) = Z a Z b e r (1 − exp ( − πr/ Λ ab ))+ k B T ln 2 exp (cid:0) − πr / Λ ab / ln 2 (cid:1) δ ae δ be (5)where Λ ab = (cid:112) π (cid:126) /µ ab k B T , µ ab is the reduced mass,and T = T e except when a and b are both protons when T = T p . The potentials are temperature-dependent, butwere held constant in most of our short simulations. Forlong simulations to equilibration, we allow the tempera-ture parameters to evolve with time, using a smoothedexponential average of the instantaneous MD value.Simulations were run long enough to extract a re-laxation time (typically 10% of τ ∗ ), with some strong-coupling cases continued to complete equilibration. Weobtain τ pe by fitting the temperature over a brief interval, dT e dt = T p − T e τ ep ; dT p dt = T e − T p τ pe . (6)We choose the timestep to conserve total energy overthe entire simulation (cid:0) ∆ E/E < − (cid:1) when using fixedpotentials. Typically, ∆ t ranges from 5 × − to 10 − fs.Any drift in the energy is tightly controlled, as artificialheating can distort the true relaxation rate. In long runs,the potentials change slowly as the temperature relaxes.Although energy is not conserved in these cases, the totalenergy change remains less than 3%. In practice, τ pe calculated from the time-dependent potential is within10-15% of the result for the constant potential.Convergence with respect to system size is tested byemploying various particles numbers N ranging from Case n i (1 /cc ) T e ( eV ) T i ( eV ) τ ∗ ( fs ) σ ( fs )A 10 . . . × . × B 10 . . . × . × C 10 . . . × . × D 10 . . . × . × E 10 . . . × . × F 10 . . . × . × G 10 . . . × . × H 10 . . . × . × I 10 . . . × . × J 1 . × . . . × .
3K 1 . × .
47 12 . . × . × L 10 . . . × . × M . . . × . × M . . . × . × M . . . × . × TABLE I: Density, initial electron, and ion temperature, re-laxation time and standard deviation of the MD simulations. N = 128 (the number that HM employed), to as manyas N = 64 , N = 1024.Statistical uncertainty for each case is estimated by com-puting the relaxation rate from equivalent samples (from8 to 64) of a microcanonical ensemble and then takingthe average and standard deviation. Sensitivity to initialconditions is studied using the ensemble of simulationsand/or by discarding a portion of the initial temperatureevolution.The nonequilibrium system is prepared using two sep-arate Langevin thermostats for protons and electrons.Initial configurations are sampled from a stationary dis-tribution obtained after 10 − timesteps. The ther-mostats are then removed, and the species allowed to un-dergo (microcanonical) collisional relaxation for approx-imately 10 timesteps.Equations 6 are valid for an ideal gas equation of statefor the plasma. For strongly coupled plasmas there is asignificant potential energy contribution that would in-validate this assumption. However, the error associatedwith using the temperature evolution equations is smallin the temperature-density regimes of interest here[16].Although the MD temperature relaxation is asymmet-ric in the strong-coupling cases, we find | dT e /dt | and | dT p /dt | differ by only about 10%. Thus, we only report1 /τ ∗ = 1 /τ pe + 1 /τ ep ≈ /τ pe .Table 1 lists the set of initial conditions for 15 dif-ferent systems. The ensemble average temperature re-laxation, τ ∗ (calculated from d ∆ T /dt = ∆
T /τ ∗ ), andthe standard deviation, σ , are in femtoseconds. A rangeof initial conditions were chosen to span the weakly- tostrongly-coupled and the degenerate to non-degenerateregimes. We include two sets of initial conditions con-sidered by HM (Cases J and K). In most cases, hydro-gen plasma is simulated using the true electron-protonmass ratio of 1:1836. In Case L, the cold electrons werereplaced with cold protons in order to shorten the re-quired simulation time. Cases M − involve a compari-son of electron-proton and positron-proton systems andwill be discussed below. Cases F and G have degenerateelectrons. Degeneracy effects are treated in neither theclassical MD simulations nor in the LS, GMS6, or BPSmodels. Hence, the models can be directly compared tothe simulations even for those cases when comparisonswith experiment would be questionable. III. COMPARISON WITH THEORY
Figure 1 shows MD results for Case K run to near-full relaxation using potentials that are implicitly time-dependent (temperature-dependent). We also displaypredictions for LS, GMS6, and BPS. The MD data inFigure 1 is most closely matched by BPS (although thisis partly fortuitous, as will been seen below) followed byGMS6. The LS model predicts the fastest relaxation, ex-ceeding MD by about a factor of two. This disagreementcontradicts the conclusion reached by HM. At the sametime, our τ ∗ for Cases J and K agree with those reportedby HM. We attribute the discrepancy to inconsistent def-initions of τ pe , τ LS and τ ∗ : τ LS is properly equal to τ pe ,which is 2 τ ∗ (not τ ∗ ) if τ pe ≡ τ ep . As previously notedby HM, however, ambiguities in the b min and b max maybe sufficient to accommodate this difference.To make comparisons of our MD results with theoret-ical predictions more transparent, we define an effectiveCoulomb logarithm as L MD ≡ J LS /τ ∗ . This result isthen compared with the theoretical prediction for L com-ing from LS, GMS6, and BPS. Fig. (2) shows simulationresults for L MD with error bars along with theoreticalpredictions for L GMS (solid) and L BP S (dashed) as afunction of initial electron temperature. Numerical re-sults and analytic expressions for L are arranged accord-ing to density; n = 10 , , and 10 (blue, red andblack respectively).In regions where it is expected to be applicable, wefind that LS systematically overestimates the effectiveCoulomb logarithm and thus predicts a relaxation ratethat is too fast relative to the MD results. For plasmaswith L >
1, the MD results are consistent with boththe GMS6 and BPS, suggesting that approaches beyondLS are indeed more predictive. As expected, BPS in-creasingly underestimates the relaxation rate for L < L over a surprisingly broad range of densityand temperature. Further discrimination between thesetheories in the region where where they are expected tobe most accurate (low density and high temperature) is time (fs) T e m p e r a t u r e ( e V ) MD (8 runs)MD (ave)LSGMS6BPS
FIG. 1: Electron (top curves) and proton (bottom curves)temperature relaxation is shown based on MD, GMS, LS, andBPS for Case L. The MD results are shown by points fromseveral simulations, with a line through the average. Notethat all approaches relax slower than LS.
10 100 1000
Temperature (eV) L GMS6BPSLSMD cc -1 cc -1 cc -1 FIG. 2: Theoretical (GMS6 [solid], BPS [dashed], and LS[dotted]) and MD calculations of L as a function of initial T e for densities 10 /cc, /cc, and 10 /cc (blue, red andblack respectively). Additional detail is in the text. not possible given the large uncertainties present in ourcurrent MD simulations. However, our results suggestthat validation of these theories could be accomplishedwith carefully controlled experiments [24] and larger (andlonger) simulations that further reduce statistical error.Finally, LS predicts identical equilibration rates forlike-charge and opposite charge systems. We tested thisby performing three sets of simulations at the same den-sity and temperatures (Case M − in Table I.) We simu-lated electrons-protons (M ) and positrons-protons (M )using Equation 6, and positrons-protons using a pure1 /r Coulomb potential (M ). The relaxation rates forall three cases agree to within our error bars, suggest-ing that energy transfer in these systems is occurringpredominately on length scales longer than the thermaldeBroglie wavelength. IV. CONCLUSIONS
We have performed MD simulations of the tempera-ture relaxation process in hot, dense hydrogen. We in-vestigated systems containing as large as 64,000 parti-cles, finding that N (cid:39) L (cid:38)
1, the sim-ulations are consistent with both GMS6 and BPS. In con-trast, the LS approach systematically overestimates therelaxation rate. In the limit of high temperature and lowdensity, all models are in agreement, however. Our MDresults suggest that LS is accurate for L >
4, rather thanthe usual restriction of L >
10, in agreement with pre-vious work [8, 21]. More modern approaches exemplifiedhere by GMS6 and BPS clearly extend the accessible pa-rameter space closer to
L ∼
1, with GMS6 providing areasonable description of the MD data even for for L < V. ACKNOWLEDGMENTS
This work performed under the auspices of the U.S.Department of Energy by Lawrence Livermore NationalLaboratory under Contract DE-AC52-07NA27344. Thiswork was funded by the Laboratory Directed Researchand Development Program at LLNL under project track-ing code 07-ERD-044. [1] P. Celliers, A. Ng, G. Xu, and A. Forsman, Phys. Lett ,2305 (1992); A. Ng, P. Celliers, G. Xu, and A. Forsman,Phys. Rev E , 4299 (1995).[2] L. D. Landau, Phys. Z. Sowjetunion , 154 (1936); Sov.Phys. JETP , 203 (1937).[3] L. Spitzer, Jr., Physics of Fully Ionized Gases , 2nd ed.(Interscience, New York, 1962).[4] Y. T. Lee and R. M. More, Phys. Fluids , 1273 (1984);J. D. Jackson, Classical Electrodynamics (Wiley, NewYork, 1975); R. Cauble and W. Rozmus, Phys. Rev E , 2974 (1995).[5] E. A. Frieman and D. C. Book, Phys. Fluids , 1700(1963).[6] H. A. Gould and H. E. DeWitt, Phys. Rev , 68 (1967).[7] M.W.C. Dharma-wardana and F. Perrot, Phys. Rev. E , 3705 (1998).[8] D. O. Gericke, M. S. Murillo, and M. Schlanges, Phys.Rev. E , 036418 (2002).[9] L. S. Brown, D. L. Preston, and R. L. Singleton, Jr.,Phys. Rep. , 237 (2005).[10] This is model 6 in Ref. [8]. Note that we have recastEquations 3 and 4 to use a consistent definition of thethermal deBroglie wavelength.[11] L. S. Brown and R. L. Singleton, Jr., Phys. Rev. E ,066404 (2007).[12] J. P. Hansen and I. R. McDonald, Phys. Lett. , 42 (1983).[13] F. H. Streitz, J. N. Glosli, and M. V. Patel, Phys. Lett.96, 225701 (2006).[14] P. Ewald, Ann. Phys.
253 (1921).[15] S.W. DeLeeuw, J.W.Perram and E.R. Smith, Proc. Roy.Soc. Lond.
A373 , 27 (1980).[16] D. O. Gericke and M. S. Murillo, Inertial Fusion Sciencesand Applications ; 1002 (2004).[17] C. S. Jones and M. S. Murillo, High Energy DensityPhysics , 379 (2007).[18] T. Dunn and A. A. Broyles, Phys. Rev. , 156 (1967).[19] H. Minoo, M. M. Gombert and C. Deutsch, Phys. Rev A , 1544 (1981).[20] G. Kelbg, Ann. Phys. , 219 (1963).[21] C.-K. Li and R. D. Petrasso, Phys. Rev. Lett. , 3059(1993).[22] A.V. Filinov, V.O. Golubnychiy, M. Bonitz, W. Ebelingand J.W. Dufty, Phys. Rev. E , 046411(2004).[23] For densities near 10 cm − , L BPS happens to cross L MD between 50 and 100 eV, mid-range of the relaxationin Fig. 1.[24] J.M. Taccetti, R.P. Shurter, J.P. Roberts, J.F. Benage, B.Graden, B. Haberle, M.S. Murillo, B. Vigil, F.J. Wysocki,J. Phys. A;39