Accurate Determination of the Shear Viscosity of the One-Component Plasma
AAccurate Determination of the Shear Viscosity of the One-Component Plasma
J´erˆome Daligault ∗ and Kim Ø. Rasmussen Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
Scott D. Baalrud
Department of Physics and Astronomy, University of Iowa, Iowa City, Iowa 52242, USA (Dated: May 10, 2019)The shear viscosity coefficient of the one-component plasma is calculated with unprecedentedaccuracy using equilibrium molecular dynamics simulations and the Green-Kubo relation. Numericaland statistical uncertainties and their mitigation for improving accuracy are analyzed. In the weaklycoupled regime, our the results agree with the Landau-Spitzer prediction. In the moderately andstrongly coupled regimes, our results are found in good agreement with recent results obtained forthe Yukawa one-component plasma using non-equilibrium molecular dynamics. A practical formulais provided for evaluating the viscosity coefficient across coupling regimes, from the weakly-coupledregime up to solidification threshold. The results are used to test theoretical predictions of theviscosity coefficients found in the literature.
PACS numbers: 52.27.Gr,52.25.Fi,52.27.Lw
I. INTRODUCTION
Like the hard-sphere model in the theory of simple liq-uids, the classical one-component plasma (OCP) is a ref-erence model in the study of strongly coupled Coulombsystems and, in particular, of ions in strongly coupledplasmas [1]. By definition, the OCP consists of a sys-tem of identical ions of charge Ze , mass m and numberdensity n in an infinite three-dimensional space. Parti-cle dynamics is governed by the laws of classical, non-relativistic mechanics. The interaction energy betweentwo ions separated by the distance r is modeled by aYukawa potential v ( r ) = q e − r/λ sc r , where λ sc ≥ q = ( Ze ) / π(cid:15) . In the limit λ sc → + ∞ ,particles interact via the bare Coulomb interaction andthe ions must be immersed in a uniform, neutralizingbackground for well-posedness of the model.The equilibrium properties of the OCP depend ononly two dimensionless parameters: the screening pa-rameter κ = a/λ sc and the Coulomb coupling parameterΓ = q /ak B T , where a = (3 / πn ) / is the Wigner-Seitzradius and T is the temperature. The Coulomb couplingparameter measures the degree of non-ideality of the sys-tem, i.e. the degree to which many-body interactions af-fect the properties of the ensemble of ions. Given a valuefor κ , the OCP shows transitions from a nearly collision-less, gaseous regime for Γ << m (e.g., Γ m (cid:39)
175 at κ = 0, Γ m = 440 at κ = 2). The gas-like to liquid-likecrossover manifests itself in several ways in the micro- ∗ Electronic address: [email protected] scopic properties of the OCP. Most noticeably, the coef-ficient of reduced shear viscosity η ∗ = ηmna ω p , (1)where ω p is the plasma frequency defined below, exhibitsa minimum at intermediate values around Γ min ∼ Γ m / η increases monotonically with density along any isotherm,whereas along any isochore, η exhibits a minimum as afunction of temperature. In a fluid, transport of momen-tum occurs not only by the bodily movement of parti-cles, but also by the direct transmission of intermolecularforces, which results from a competition between kineticand interaction effects. At small coupling Γ (cid:28) Γ min ,the former mechanism is predominant and, like in a gas,the OCP viscosity increases with increasing temperature.At large coupling Γ (cid:29) Γ min , the latter mechanism ispredominant and, like in a liquid, the OCP viscosity de-creases with increasing temperature. Strong interparticleinteractions give rise to the cage-effect [4], whereby eachparticle finds itself trapped for some period of time in thecage formed by its immediate neighbors, rebounding untilit overcomes the energy barrier and diffuses to a neigh-boring cage. At intermediate coupling Γ ∼ Γ min , the twomomentum transport mechanisms contribute with simi-lar magnitude, resulting in a shallow minimum in theviscosity coefficient.Despite the apparent simplicity of the OCP model, ac-curate determination of the viscosity coefficient of theOCP by molecular dynamics (MD) simulations is diffi-cult. This is exemplified by the significant discrepanciesin the results obtained over the years by different authors.A compilation of these has been provided in [5]. Remark-ably, important differences are found not only betweenresults obtained using different MD techniques, but alsobetween those obtained using the same technique. Forexample, previous results obtained using equilibrium MD a r X i v : . [ phy s i c s . p l a s m - ph ] J u l −1 Coupling Strength ( Γ ) S hea r V i sc o s i t y ( η * ) Vieillefosse/HansenBernu/VieillefosseWallenborn/BausSalin/CaillolBastea
FIG. 1: Shear-viscosity coefficient of the Coulomb OCP( κ = 0) obtained by different authors [6–9] using equilibriummolecular dynamics and the Green-Kubo relation. Resultsobtained using non-equilibrium molecular dynamics are com-piled in [5]. are shown in Figure 1. Among these, the results of Bastea[6] are believed to be most accurate. Recently, Donk´oand Hartmann have presented arguably the most accu-rate results for moderately and strongly coupled OCP’sat κ = 1 , , non-equilibrium MDsimulation methods, namely the M¨uller-Plathe reverseMD approach and the Evans-Morriss homogeneous shearalgorithm. In the present paper, we use equilibrium
MDbased on the evaluation of the Green-Kubo relation tovalidate the non-equilibrium MD results of Donk´o andHartmann and the equilibrium MD results of Bastea forthe Coulomb OCP ( κ = 0).There are a couple of reasons why this is important.Donk´o and Hartmann chose non-equilibrium MD meth-ods, claiming that they are generally more efficient thanequilibrium calculations. Indeed, we shall see that thedetermination of viscosity from the Green-Kubo relationis made difficult by the large statistical imprecision inthe calculation of the shear-stress autocorrelation func-tion. This arises mainly due to the fact that simulationaverages are taken over finite-length runs. The noise canbe satisfactorily reduced at the price of very long simula-tions, but this requires much longer run times than havepreviously been reported. Despite the computationalcost, equilibrium MD has advantages. It provides in-formation about the microphysical ion dynamics, in par-ticular the time-correlation function of the shear stress.It works equally well for all κ and Γ values, includingthe Coulomb OCP. The same simulation can be used toconsistently calculate all other transport and static prop-erties of the plasma. Furthermore, it provides an inde-pendent method that allows us to confirm Donk´o andHartmann’s results.In the past, much effort has been devoted to develop a Particle number 5000 (Γ ≥ .
5) ; 50000 (Γ < . . /ω p ; 0 . /ω p (Γ < . . /ω p Equilibration length 1000 /ω p Numerical method Ewald sums with P M algorithm [10]Ewald parameter aα N = 5000) 1 (N=50000)Short-range cutoff r C /a (N=5000) 64 (N=50000)rms forces < − TABLE I: Main parameters of the equilibrium MD simula-tions used in this work to calculate the viscosity coefficients. theory that extends the traditional plasma regime validat small Γ to the moderate and strongly coupled regimes.An accurate determination of the shear-viscosity is de-sirable to test existing and future theories. In this pa-per, we test the conventional result of Landau-Spitzer,the kinetic theories of Wallenborn-Baus [8], Viellifosse-Hansen [3] and of Tanaka-Ichimaru [21], and the recenteffective potential theory of Baalrud-Daligault [18].This paper is organized as follows. Section II describesthe equilibrium simulations used to determine the shearviscosity coefficient. We present a detailed study of thestatistical convergence necessary to ensure quality of thefinal results. The simulation methods and parametersare explicitly given to help anyone who wishes to repro-duce our results. The results are described in section III.Finally, in section IV, we compare the results with theabove mentioned theories.
II. SIMULATION METHODS AND RESULTS
In the following, ω p = (4 πnq /m ) / denotes theplasma frequency. In our MD simulation 1 /ω p and theWigner-Seitz radius a are used as unit of time and length,respectively. A. Basic definitions
Here we describe a typical equilibrium MD simulationof an OCP at given Γ and κ . N particles are placed in acubic box of volume V = L and periodic conditions areimposed on all boundaries. Particle trajectories are de-termined by solving Newton’s equations of motion withthe velocity Verlet integrator [11]. The force on an ionthat results from its interaction with the ions in the sim-ulation box and with those in the periodically replicatedcells is calculated using the Ewald summation technique.This is essential for small κ sc := a/λ sc values becausethe range of the interaction is larger than the simulationbox in this case. A formulation of the Ewald sum ap-proach for Yukawa potentials can be found in [9]. Theinteraction energy v ( r ) between two particles at distance r is represented by a sum of a short-range (sr) and along-range (lr) component, v ( r ) = q φ sr ( r ) + q φ lr ( r ) (2)where φ sr ( r ) = 12 r (cid:104) erfc (cid:16) α r + κ sc α (cid:17) e κ sc r + erfc (cid:16) α r − κ sc α (cid:17) e − κ sc r (cid:105) and φ lr ( r ) =4 πV (cid:88) n ∈ Z e − ( k + κ sc ) / (4 α ) k + κ sc e i k · r , k = 2 πL n (3)where α > M) method, which combineshigh-resolution of close encounters (the sr term is calcu-lated using nearest neighbor techniques) and rapid, long-range force calculations (the lr forces are computed ona mesh using three-dimensional fast Fourier transforms)[10].Table I lists the main numerical parameters used in thepresent study to compute the viscosity coefficients. Ourtimestep δt = 10 − /ω p is chosen small enough to ensureexcellent energy conservation for all Γ and κ values. Inall simulations, N = 5000 for Γ > .
5, while N = 50000was chosen for Γ < . B. Shear viscosity coefficient
The shear viscosity coefficient, η , was computed usingthe Green-Kubo relation that expresses η as the timeintegral of the equilibrium autocorrelation function of theoff-diagonal components σ xy of the shear stress tensor σ ↔ [12], η = 16 V k B T (cid:88) x =1 3 (cid:88) y (cid:54) = x =1 (cid:90) ∞ J xy ( t ) dt . (4)where J xy ( t ) is the shear-stress autocorrelation function J xy ( t ) = (cid:10) σ xy ( t ) σ xy (0) (cid:11) eq . (5)In Eq.(4), the brackets (cid:10) . . . (cid:11) eq the denote equilibrium(thermal) average at temperature T . As shown Ref. [9],it follows from the Ewald decomposition of the interac-tion potential that the components of the shear stresstensor can be conveniently split into a kinetic compo-nent, a short-range interaction component, and a long-range interaction component σ ↔ ( t ) = σ ↔ kin ( t ) + σ ↔ sr ( t ) + σ ↔ lr ( t ) , (6)where σ ↔ kin ( t ) = N (cid:88) i =1 m v i ( t ) v i ( t ) (7) σ ↔ sr ( t ) = − q N (cid:88) i (cid:54) = j =1 (cid:88) n ∈ Z ( r ij ( t ) + n L )( r ij ( t ) + n L ) || r ij ( t ) + n L || φ (cid:48) sr ( || r ij ( t ) + n L || ) (8) σ ↔ lr ( t ) = q V N (cid:88) i (cid:54) = j =1 (cid:88) n ∈ Z (cid:32) ˜ φ lr ( k ) ↔ U + d ˜ φ lr ( k ) dk kk k (cid:33) e i k · r ij ( t ) , k = 2 πL n (9)in which r i ( t ) and v i ( t ) are the instantaneous postionand velocity of particle i at time t , r ij = r i − r j , ↔ U is theunit dyad tensor, and ˜ φ lr ( k ) = 4 π e − ( k κ sc ) / (4 α k + κ sc . C. Typical MD simulation
The simulations were performed as follows. Initial par-ticle positions were assigned randomly in the simulationbox, with a small region surrounding each particle ex-cluded to avoid initial explosion. Initial particle veloci-ties were assigned randomly from a Maxwell-Boltzmann distribution at the desired temperature. The simula-tion time consisted of an equilibration phase of length t eq = N eq δt followed by the main MD run of length t run = N run δt , for a total of N eq + N run time steps. Dur-ing the equilibration phase, velocity scaling (also knownas the Berendsen thermostat [11]) was used at everytimestep to maintain the desired temperature. Veloc-ity scaling was turned off after the equilibration phase,at which point the simulations transitioned to the mainMD run phase in which particle positions and velocitieswere recorded at every timestep.The viscosity coefficient was evaluated from the Green-Kubo relation, Eq. (4), as follows. First, σ xy ( t ) was t t ω p J xy ( t ) ω p η ( t ) t t t ω p J xy ( t ) ω p η ( t ) t t t ω p J xy ( t ) ω p η ( t ) t FIG. 2: Illustration of the calculation with equilibrium MD of the Green-Kubo expression (4) for the shear viscosity coefficient.Shown are results for the OCP with κ = 0 and for three different values of the coupling parameter Γ: Γ = 0 . (cid:80) x (cid:54) = y (cid:104) σ xy ( t ) σ xy (0) (cid:105) obtained using Eq.(13) are shown in the left-hand column. Corresponding (dimensionless) cumulated sums η ( t ) defined by Eq.(15) are shown in the right-hand column. In each figure, the results corresponding to five different simulationlengths are shown (in 1 /ω p units); the other simulation parameters are listed in table I. In the left-hand side figures, the curveshave been shifted vertically for clarity (the horizontal dashed lines corresponds to the shifted zero of the vertical axis.) computed using Eqs. (7)-(9) at each time step t = nδt ,0 ≤ n ≤ N run using the positions and velocities from theMD simulation. Second, the shear stress autocorrelationfunction (5) was computed by first replacing the thermal average by a time average J xy ( τ ) = lim t →∞ ¯ J xy ( t, τ ) (10) (cid:39) ¯ J xy ( t run , τ ) (11)where¯ J xy ( t, τ ) := 1 t − τ (cid:90) t − τ σ xy ( s + τ ) σ xy ( s ) ds (12)and then discretizing in time¯ J xy ( t run = N run δt, τ = nδt ) (13) (cid:39) N run + 1 − n N run − n (cid:88) m =0 σ xy (( m + n ) δt ) σ xy ( mδt ) . Third, the cumulated sum η ( τ = nδt ) := 16 V k B T (cid:88) x =1 3 (cid:88) y (cid:54) = x =1 (cid:90) τ (cid:10) σ xy ( t ) σ xy (0) (cid:11) eq dt (14) (cid:39) n (cid:88) m =0 δt V k B T (cid:88) x (cid:54) = y ¯ J xy ( N run δt, mδt ) (15)was calculated. Ideally, according to the Green-Kuborelation (4), the viscosity coefficient is given by η =lim t run →∞ η ( t run ) (neglecting other systematic errorsdue to, e.g., size effects, N , force accuracy, etc.) Oneexpects the sum (14) to converge towards η after a timelonger than the correlation time scale of the correlationfunction. Beyond that time, the correlation function van-ishes and the cumulated sum reaches a plateau valueequal to the viscosity coefficient. In practice, as we shallsee, the convergence to a plateau is quite slow and the ac-curate determination of the viscosity is impossible unlessone performs very long simulations.Figure 2 illustrates the method and its convergence fora Coulomb OCP ( κ = 0) at three values of the couplingparameter across the fluid regime: Γ = 0 . ,
30 and 120.The figures on the left-hand side show the shear-stress au-tocorrelation function calculated using Eq. (13), those onthe right-hand side show the cumulated sum Eq. (15). Ineach case, the results of simulations for four different sim-ulation lengths are shown, namely t run = t , t , t , t and 16 t with t = 5242 . /ω p , all the other numericalparameters being identical (see table I). The final resultsreported in Sec. III were all obtained using t run = 16 t ;note also that the shortest runs shown here are in factlonger than the simulation duration ∼ − /ω p pre-viously reported in the literature [6, 9]. For all Γ values,the autocorrelation function decays toward zero on a timescale t c ∼ (cid:48) s /ω p , which is much smaller than the simu-lation time scale since t c (cid:28) t . Nevertheless, the lengthof the simulation has a significant effect on the noise and,in turn, on the convergence of the Green-Kubo calcula-tion. Thus, in all cases, the correlation function obtainedwith the shortest MD run t run = t never fully vanishesas time increases: it decays toward zero and then alterna-tively stays above and below the horizontal axis. Whenintegrated over time, this leads to a significant noise inthe cumulated sum, which never quite reaches a plateauvalue. J xy ( ) t run FIG. 3: Illustration of the convergence of the sum rule (16)with the length of the simulation t run for Γ = 1 (upper panel),Γ = 20 (middle panel), and Γ = 175 (lower panel). In all cases κ = 0. Dashed lines represent Eq. (16) evaluated with the g ( r )obtained from the MD simulation, while the dots representsresults of direct calculations of J xy (0) from the stress tensor. D. Convergence Study
In this section, we report on the analysis of the speed ofconvergence of the calculation that we have undertakento select the numerical parameters of table I used to cal-culate the viscosity coefficients reported in section III.The goal is to empirically answer the question: how largeshould the simulation length t run be in order to attainthe desired accuracy.
1. Initial time behavior
We start with a calculation of the initial value of thecorrelation function: J xy (0) = (cid:104) σ xy ( t ) σ xy ( t ) (cid:105) . Accuracyin the determination of J xy (0) is important since a shiftin its value would certainly correspond to a shift of theentire time-evolution, which would lead to error in the cu-mulated sum and viscosity. We find that although J xy (0)is less subject to statistical noise than J xy ( t ) for t >
0, itis sufficient to cause concern. Remarkably, a number ofproperties (or “sum rules”) concerning J xy (0) are knownand can be used to monitor the converge of its numeri-cal determination. In particular, the following exact sumrule can be shown J xy (0) = N ( k B T ) (16)+ 2 πN nk B T (cid:90) ∞ drr ( g ( r ) − δ κ, ) [4 φ (cid:48) ( r )+ rφ (cid:48)(cid:48) ( r )] −2 −1 N/(3 Γ ) ~ Γ −1.2 Coupling Strength ( Γ ) B u l k M odu l u s FIG. 4: (color online) Bulk modulus of the Coulomb OCPas a function of Γ. The vertical axis shows the dimensionlessquantity ˜ J xy (0) = J xy (0)[ m ( aω p ) ] = N (3Γ) Gnk B T . Alternatively, using unitless quantities ˜ σ = σ/ ( m ( aω p ) )and ˜ v ( r ) = v ( r/a ) /k B T ,˜ J xy (0) = J xy (0)[ m ( aω p ) ] (17)= N (18)+ N (cid:90) ∞ drr [ g ( r ) − δ κ, ] [4˜ v (cid:48) ( r ) + r ˜ v (cid:48)(cid:48) ( r )] . Figure 3 shows a plot of the evolution of J xy (0) withthe length of the simulation for a Coulomb OCP atΓ = 1 ,
20 and 175. In all cases, the initial value of thecorrelation function equals the expected value for simula-tion lengths greater than t ∗ ∼ /ω p . The inaccuracygrows rapidly when the simulation length lies below thisvalue. The convergence of the correlation function attimes t ≥ J xy (0) as afunction of the coupling strength Γ. Note that this quan-tity is simply related to the isothermal bulk modulus G = − V dPdV = n dPdn , such that J xy (0) = V k B T G. (19)In the weakly coupled regime, the equation of state isdominated by the ideal gas, kinetic pressure P = nk B T ,i.e. G = nk B T and J xy (0) = N ( k B T ) ; accordingly˜ J xy (0) (cid:39) N . In the strongly coupled regime, the inter-action energy dominates the pressure and the MD datashow that ˜ J xy (0) scales like ∼ N Γ . .A more detailed study of the convergence of the initialvalue can be obtained by considering the different com-ponents of the sum rule (16) obtained by substituting thedecomposition (6) in Eq.(5). The first term in Eq. (16) K i n − K i n ( P o t − P o t ) s r −0.0200.02 t run P o t − K i n FIG. 5: Illustration of the convergence of the detailed com-ponents of the sum rule (16) with the length of the simula-tion t run for Γ = 2 and κ = 0. Dashed lines represent theexact value, the dots represents results of direct calculationsfrom the components of the stress tensor. Top panel: kinetic-kinetic term (20). Middle panel: sr potential - potential term(23). bottom panel: kinetic-potential term (25) corresponds to the kinetic-kinetic contribution J kinxy (0) = (cid:10) σ kin xy (0) σ kin xy (0) (cid:11) eq = N ( k B T ) . (20)As shown in the appendix, the direct evaluation of J kinxy (0) using Eq. (13) in an MD simulation amounts tocalculating J kinxy (0) (cid:39) (cid:34) t run (cid:90) t run N (cid:88) i =1 mv x,i ( t ) dt (cid:35) × (cid:34) t run (cid:90) t run N (cid:88) i =1 mv y,i ( t ) dt (cid:35) + cross term . (21)where the full expression for the cross term is given inEq.(A1). In the large t run limit, the first term, whichis related to the product of the averaged instantaneouskinetic energy (cid:80) Ni =1 12 mv x,i ( t ) , is expected to convergeto the exact result (20), while the remaining cross termis expected to vanish. In a MD simulation, the instanta-neous kinetic energy fluctuates around the target velocityand the first term rapidly converges with t run to its limitvalue. On the contrary, the cross term, which involvescontributions that are quartic in the velocities, does notconverge as fast as the kinetic energy to its limit value,which is quadratic in the velocities. This is illustrated inFig. 5 (top panel): the value of (21) converges to the ex-pected value at large enough t run beyond which the crossterms are negligibly small compared to the first term.We now discuss the term involving the interaction only, J potxy (0) := (cid:10) (cid:2) σ sr xy (0) + σ lr xy (0) (cid:3) (cid:2) σ sr xy (0) + σ lr xy (0) (cid:3) (cid:11) eq (22)= 2 πN nk B T (cid:90) ∞ drr g ( r ) [4 φ (cid:48) ( r ) + rφ (cid:48)(cid:48) ( r )]The later can actually be further broken into two com-ponents; for instance [13], J srxy (0) := (cid:10) σ sr xy (0) (cid:2) σ sr xy (0) + σ lr xy (0) (cid:3) (cid:11) eq (23)= 2 πN nk B T (cid:90) ∞ drr g ( r ) [4 φ (cid:48) sr ( r ) + rφ (cid:48)(cid:48) sr ( r )] . (24)Figure 5 shows the convergence of J srxy (0) with the simu-lation length t run toward the exact value Eq.(24). Againwe find that long t run must be used to ensure conver-gence.Finally, the kinetic-potential term J kin − potxy (0) := (cid:10) σ kin xy (0) (cid:2) σ sr xy (0) + σ lr xy (0) (cid:3) (cid:11) eq (25)= (cid:10) (cid:2) σ sr xy (0) + σ lr xy (0) (cid:3) σ kin xy (0) (cid:11) eq = 0 . (26)In practice, the term is negligibly small for long enoughsimulation length t run .In conclusion, the initial value of the correlation func-tion converges relatively slowly towards its expectedvalue; simulations longer than t ∗ are necessary to repro-duce the expected value. The slow convergence is foundto be caused by cross terms that vanish in the ideal limitbut are finite in practice.
2. Finite time correlation function
The initial time correlation function determines a lowerbound, t ∗ , for the simulation length needed to calculatethe viscosity coefficient. However, Fig. 2 shows that thisis far too short to obtain the accurate correlation func-tions necessary to evaluate the viscosity coefficient. Theintermediate time dynamics of J xy ( t ) does not convergeas fast as the short-time dynamics, which leads to largevariations in the evaluation of the viscosity coefficient.As a consequence, the cumulated sum does not reach aplateau value at time t ∗ . It is noteworthy that the time t ∗ is actually larger than the simulation lengths used inprevious studies [9]. Figure 2 reveals that satisfying con-vergence can be achieved for simulation times on the or-der of ∼ t ∗ , corresponding to over 8 million time stepsfor δt = 0 . /ω p .In order to understand this behavior, we employ thestatistical error analysis of Zwanzig and Ailawadi [11, 14].Zwanzig and Ailawadi gave an error estimate for the de-viation between the shear-stress autocorrelation functionat time t obtained with an MD simulation of finite length t run , and its exact value ¯ J xy ( ∞ , t )∆( t ) = ¯ J xy ( t run , t ) − ¯ J xy ( ∞ , t ); (27) see Eq. (12). To this end, they assumed that σ xy ( t ) is aGaussian random variable (average denoted by (cid:104) . . . (cid:105) be-low), which was shown to give the correct order of magni-tude of error estimates [11, 15]. Under this assumption,they arrived at the result (cid:104) ∆( t )∆( t ) (cid:105) (cid:39) τ c t run ¯ J xy ( ∞ , (28)with 0 ≤ t , t ≤ τ c (cid:28) t run where τ c = 2 (cid:82) ∞ dt J xy ( t ) J xy (0) (29)measures the relaxation time within which the exact cor-relation function ¯ J xy ( ∞ , t ) decays to zero from its initialvalue. Applying Eq. (28) with t = t = t shows thatthe absolute error (cid:104) ∆( t ) (cid:105) in J xy ( t ) is independent of t .Therefore, the relative error (cid:68)(cid:2) ¯ J xy ( t run , t ) − ¯ J xy ( ∞ , t ) (cid:3) (cid:69) ¯ J xy ( ∞ , t ) (cid:39) τ c t run (cid:20) ¯ J xy ( ∞ , J xy ( ∞ , t ) (cid:21) (30)increases rapidly as ¯ J xy ( ∞ , t ) goes to 0 [11]. Equation(30) shows that this increase in the relative error can belessened by increasing the simulation length t run . This isindeed consistent with the results in Fig. 2 for t run = t and t run = 2 t , although the initial value of the corre-lation function has converged to its expected value, thevalues at later times t ≤ τ c are not converged. Thisnoise is reduced when the simulation length is increasedto 16 t . III. FINAL RESULTS
The shear viscosity coefficients obtained using themethod described in section II, along with the numer-ical parameters collected in table I, are shown in figure 6for κ = 0 and κ = 2. Also shown are the data of Donk´oand Hartmann, obtained from non-equilibrium MD for κ = 2 (see table I in [5]) and the data of Bastea obtainedfrom equilibrium MD for κ = 0 (we plot the fitting for-mula (11) of [6]). For convenience, the numerical valuesare given in table II for κ = 0 and in table III for κ = 2.We highlight the following important features of thepresent results.(1) In Fig. 6, the data at κ = 2 are compared with theresults of Donk´o and Hartmann obtained using two in-dependent non-equilibrium molecular dynamics calcula-tions [5]. We find very good agreement between all threeindependent calculations. For κ = 0, our data are in verygood agreement with Bastea’s fit, which was obtained byinterpolating MD data over 0 . ≤ Γ ≤
100 [6]; see Fig. 7.(2) As shown in Fig. 7, for all Γ <
10, the viscosity coef-ficient of the coulomb OCP ( κ = 0) is well approximated Γ η/ ( mna ω p ) Γ η/ ( mna ω p ) Γ η/ ( mna ω p )0.1 75.2 20 0.084 75 0.1520.5 3.6907 21 0.084 80 0.1600.7 2.1546 23 0.085 85 0.1681 1.1831 25 0.085 90 0.1762 0.4440 27 0.086 95 0.1843 0.2755 30 0.088 100 0.1914 0.1928 32 0.089 105 0.1995 0.1713 35 0.092 110 0.2077 0.1345 40 0.097 115 0.21410 0.101 45 0.102 120 0.22212 0.0953 50 0.110 130 0.23615 0.0864 55 0.119 140 0.25117 0.083 60 0.128 175 0.320918 0.0830 65 0.136 200 0.440019 0.0810 70 0.144TABLE II: Shear viscosity coefficient η of the one-componentplasma with κ = 0 at various coupling parameters Γ as ob-tained with the molecular dynamics simulations described inthe main text. Data are shown in units of mna ω p . by η = η δ ln (cid:16) C λ D r L (cid:17) (31)where η = 54 (cid:114) mπ ( k B T ) / q , λ D = (cid:112) πq n/k B T is theDebye length, and r L = q /k B T is the so-called distanceof closest approach. Here δ = 0 .
466 and C = 1 .
493 arenumerical parameters determined by interpolating thenumerical data. The model (31) represents a straight-forward modification of the traditional Landau-Spitzer(LS) formula [16] η LS = η (cid:16) λ D r L (cid:17) (32)derived for weakly-coupled plasmas. Indeed, in theweakly-coupled limit, Eq.(31) reduces to δη / ln (cid:16) C λ D r L (cid:17) .In the LS theory, the Coulomb logarithm ln (cid:16) λ D r L (cid:17) arisesfrom the long-range nature of the Coulomb force. Itis usually expressed in terms of the Debye length λ D (which represents the largest impact parameter beyondwhich interactions are screened out), and of the distance r L (which characterizes the smallest impact parameter).Our MD simulations reveal that, while the LS theoryprovides the right scaling at Γ <<
1, the model must becorrected through the coefficients C and δ to match thedata. The coefficient C is a correction to the somewhatarbitrary parameters λ D and R L , which can be predictedby more advanced theories ([18] and literature therein).The prefactor δ is a correction to the fact that LS cor-responds to a single Sonine polynomial approximation Γ η/ ( mna ω p ) Γ η/ ( mna ω p ) Γ η/ ( mna ω p )2 0.8638 102 0.0654 242 0.111712 0.1170 112 0.0665 262 0.117032 0.0619 122 0.0736 282 0.124242 0.0584 132 0.0728 302 0.129652 0.05572 142 0.0742 322 0.131662 0.0562 162 0.0840 342 0.142672 0.05882 182 0.0906 362 0.147882 0.0575 202 0.0955 382 0.155092 0.0637 222 0.1003 402 0.1571TABLE III: Shear viscosity coefficient η of the one-componentplasma with κ = 2 at various coupling parameters Γ as ob-tained with the molecular dynamics simulations described inthe main text. Data are shown in units of mna ω p . in the Chapman-Enskog solution of the plasma kineticequation. Figure 7 shows that the modified LS result δη / ln (cid:16) C λ D r L (cid:17) breaks down at Γ ∼ .
1, while the simplemodification (31) extends its validity to the moderatelycoupled regime up to Γ ∼
10. Remarkably, the same ex-tension of the LS theory was found to work as well forother transport processes, including the electron-ion tem-perature relaxation rate [17] and the diffusion coefficientsin mixtures [19].(3) The curve η ∗ (Γ) presents a shallow minimum thatis located in the range 18 ≤ Γ ≤
20. A more precisedetermination of the minimum is not possible with theaccuracy of the present data.(4) At high Γ, the viscosity η and the self-diffusion coef-ficient D satisfy the Stokes-Einstein relation πak B T Dη = 0 . ±
2% for Γ ≥ , (33)for κ = 0. For a detailed discussion on the Stokes-Einstein and its physical interpretation, see [4].(5) Finally, we provide a practical fit that reproducesthe viscosity coefficient across coupling regimes, from theweakly coupled regime to the solid-liquid transition, inthe form η ∗ (Γ) = ηmna ω p = a Γ / ln (cid:0) b Γ / (cid:1) a Γ + a Γ + a Γ b Γ + b Γ + b Γ + b Γ (34)In Eq.(34), we enforce the model (31) valid for Γ < η ∗ = 0 . / Γ + 0 . / Γ . +0 . . ≤ Γ ≤ −1 −2 −1 Coupling Strength ( Γ ) S hea r V i sc o s i t y ( η * ) κ =0 κ =2Donko/Hartmann, κ =2 FIG. 6: (color online) Shear viscosity coefficient across thefluid phase of the one-component plasma at κ = 0 (blue opendots) and at κ = 2 (black full dots) obtained with the equilib-rium MD simulations described in the text. The lines betweenthe dots are included to guide the eyes. At κ = 2, the red dotsshow the non-equilibrium MD results of Donk´o and Hartmann[5]. at moderate and strong coupling, it fails to reproducethe traditional Landau-Spitzer behavior in the weaklycoupled regime. IV. COMPARISON TO THEORETICALMODELS
In the previous section, we compared the MD resultswith the seminal theory of Landau-Spitzer. In this sec-tion, we test the validity of theories that have been devel-oped to predict the viscosity coefficients of the CoulombOCP ( κ = 0) in the moderately and strongly cou-pled regime, namely the theory of Vieillefosse-Hansen,the kinetic theories of Wallenborn-Baus and of Tanaka-Ichimaru, and the recent effective potential theory ofBaalrud-Daligault.The predictions of these theories are compared withour new MD results in Fig. 8. In the following, we brieflyrecall some basic facts about the various theories anddiscuss their validity with regard to the comparison withthe MD results. a b a a a b b b b −3 −2 −1 −1 Coupling Strength ( Γ ) S hea r V i sc o s i t y ( η * ) MDEq.(31)Eq.(34)a/[ Γ ln(b/ Γ )]Bastea FIG. 7: (color online) Comparison of the MD data (dots)with the fitting formula Eq.(34). The red line shows thefull expression (34), the green dashed line shows the LSlimit η ∗ (Γ) = a Γ / ln (cid:16) b Γ3 / (cid:17) , the blue line shows Eq.(31), i.e. η ∗ (Γ) = a Γ / ln (cid:16) b Γ3 / (cid:17) , the black dashed line shows theformula of Bastea [6]. A. The Vieillefosse-Hansen theory
Vieillefosse and Hansen [3] applied the framework ofthe generalized hydrodynamics formalism. Briefly, theknown short-time expansion of the transverse-current au-tocorrelation function C ⊥ ( k, t ) up to fourth-order in time t was used to build a Gaussian approximation of the −2 −1 −2 −1 Coupling Strength ( Γ ) Sh e a r V i s c o s i t y ( η ∗ ) MD dataVieillefosse-HansenWallenborn-BausTanaka-IchimaruLandau-SpitzerBaalrud-Daligault
FIG. 8: (color online) Comparison of the MD results for theshear viscosity coefficient with the corrected Landau-Spitzerprediction discussed in Sec.III, the theory of Vieillefosse-Hansen, the kinetic theories of Wallenborn-Baus and ofTanaka-Ichimaru, and the effective potential theory ofBaalrud-Daligault. See main text for a detailed comparison. C ⊥ . The coefficients ofthe Gaussian approximation depended on the first threefrequency sum-rules of C ⊥ that can be exactly writtenin terms of the pair distribution function g ( r ) and ofthe ternary distribution function g ( r, r (cid:48) ). Using the su-perposition approximation to express g in term of g ,the theory of Vieillefosse-Hansen depends on the pairdistribution g only. Figure 8 displays the results re-ported in Table III of the original paper [3]. Remark-ably the predicted viscosity exhibits a minimum as afunction of Γ around Γ = 20. At Γ = 20, Vieille-fosse and Hansen give for the reduced viscosity coefficient η ∗ = 0 . ± . η ∗ = 0 .
084 reported in Table II. However, thisgood agreement may be fortuitous since, as seen in Fig. 8,the Vieillefosse-Hansen model greatly underestimate theviscosity at all other values of Γ.
B. The Wallenborn-Baus theory
Wallenborn and Baus applied the framework of renor-malized equilibrium kinetic theory, a general kinetic the-ory of phase-space correlation functions, to derive an an-alytical model for the shear-viscosity coefficient [8]. Inthis framework, the shear-viscosity coefficient can be ex-actly expressed in terms of the only unknown of the the-ory, the so-called generalized memory function. They de-rived a sophisticated approximation for the latter that,by construction, attempts to account for (i.e. renormal-ize) the correlated motion of ions. Their approximationreduces to the Lenard-Balescu collision operator when allthe quantities involved in the memory function (e.g., thedirect correlation function) are approximated by theirweakly-coupled limiting values. They then used theirapproximate memory function to calculate the shear-viscosity coefficient across coupling regimes. The val-ues of the shear-viscosity coefficient given in the originalpaper [8] are displayed in Fig. 8. At weak coupling, theWallenborn-Baus theory agrees with the MD data, whichis consistent with the fact that the theory reduces to theLenard-Balescu result with corrections due to short-rangecorrelations, which determine the correction factor C inthe Coulomb logarithm (see Sec. III). This theory doespredict a minimum of the reduced viscosity coefficientwith a value η ∗ = 0 .
007 in fair agreement with the simu-lations, but at a coupling strength Γ ≈
8, which is belowthe MD value of (cid:39)
C. The Tanaka-Ichimaru theory
Tanaka and Ichimaru obtained a model for the shearviscosity coefficient by applying the framework of non-equilibrium kinetic theory, i.e. a theory for the tempo-ral evolution of the non-equilibrium single-particle phase-space distribution functions f ( r , p , t ). Using quasi-lineartheory, they postulate an expression for the collision op- erator by introducing the notion of static local field cor-rection G ( k ), a quantity that accounts for static correla-tions between particles. Their collision operator is C I ( f, f ) = πm (cid:90) d k (2 π ) k · ∂∂ p (cid:90) d p (cid:48) v ( k )[1 − G ( k )] | (cid:15) ( k , k · p /m ) | × δ [ k · ( p − p (cid:48) )] k · (cid:18) f ( p (cid:48) ) ∂f∂ p − f ( p (cid:48) ) ∂f∂ p (cid:48) (cid:19) , (35)where (cid:15) ( k , ω ) = 1 − v ( k )[1 − G ( k )] χ (0) ( k , ω ) is theplasma dielectric function, v ( k ) = πe k , and χ (0) ( k , ω ) = − (cid:82) d p k · ∂F/∂ p ω − k · v the density response function of the idealgas, and F the Boltzmann distribution function at tem-perature T and density n . In traditional weakly coupledplasma physics, correlations are neglected, i.e. G ( k ) is setto zero, and Eq. (35) reduces to the Lenard-Balescu colli-sion operator. By applying the Chapman-Enskog methodto lowest order in the Sonine polynomial expansion, thefollowing expression for the viscosity coefficient can beobtained [20] η T I = η T I . (36)Here, the generalized Coulomb logarithmΞ T I = 2 √ π (cid:90) ∞ dk [1 − G ( k )] k (cid:90) ∞ dz e − z | (cid:15) ( k, kv T z ) | (37)arises, where v T = (cid:112) k B T /m . This can be comparedwith Eq. (32).Tanaka and Ichimaru have presented results for η T I using Eq. (37) with a local-field correction obtained bysolving the hypernetted chain (HNC) equations with thebridge function correction of Ichimaru [20, 21]. The HNCequation gives access to the direct correlation function c ( k ), which provides G ( k ) = 1 + k B Tv ( k ) c ( k ). Reference [20]provides results for 0 . ≤ Γ ≤
20. We have evaluated η T I using the same HNC equations, including Ichimaru’sbridge function, for a wider range of values; see Fig. 8.This method agrees well with the MD data for Γ < ∼ (cid:39) .
4, Ξ
T I crosses from positive to negative values,leading to a divergence of η T I .The Tanaka-Ichimaru theory reduces to traditionalplasma physics results in the weakly coupled limit. Thesimplest (Landau-Spitzer) plasma limit can be obtainedby setting G ( k ) ≡ (cid:15) ( k, ω ) ≡ T I reduces to the traditionalCoulomb logarithm ln Λ = ln (cid:16) λ D r L (cid:17) when the usual cut-offs λ D and r L (see Sec. III) are used to regularize the k -integral in Eq. (37). The Lenard-Balescu result is ob-tained by setting G ( k ) ≡ (cid:15) . In this case, the k -integral converges at k = 0but a cutoff is necessary to regularize the remaining diver-gence at k = ∞ . This case was worked out by Braun [22],who expressed the result as a correction to the Landau-Spitzer viscosity coefficient as η LB = η LS . / ln (cid:16) λ D r L (cid:17) . (38)1 −1 Wavenumber ( ka ) L o c a l F i e l d C o rr ec t i o n [ − G ( k ) ] exp( − Γ ka )HNC0.1 0.01 0.001 Γ = 1 FIG. 9: (color online) Comparison between the local fieldcorrection obtained from the hypernetted chain (HNC) ap-proximation, and Eq. (39) for weakly coupled OCP.
Alternative to these cutoffs, simple results for the lo-cal field correction can be obtained in the weakly coupledlimit that allow analytic evaluation of the convergent in-tegral in Eq. (37). Figure 9 shows that1 − G ( k ) = exp( − Γ ka ) (39)provides a good approximation for the OCP local fieldcorrection in the weakly coupled limit. If we also takethe static dielectric function ˜ ε = 1 + 3Γ[1 − G ( k )] / ( ka ) and note that the local field correction is negligible inthis for weakly coupled plasmas [˜ ε (cid:39) / ( ka ) ], wefind Ξ T I (cid:39) (cid:90) ∞ d ¯ k ¯ k exp( − Γ¯ k )¯ k + 3Γ (40)= 12 (cid:20) E ( i/ Λ) e i/ Λ + E ( − i/ Λ) e − i/ Λ (cid:21) in which Λ = 1 / ( √ / ) = λ D /r L is the OCP plasmaparameter and E is the exponential integral. Expandingfor Λ (cid:29) T I → ln Λ − γ + O (Λ − ) (41)where γ is Euler’s constant. Note that Eq. (41) is thesame result, including the order unity correction, as hasbeen obtained from other methods, including using thescreened Coulomb potential in the effective potential the-ory ([18] and references therein). D. The effective potential theory
Recently, we proposed another approach for extendingtraditional plasma transport theories into the strong cou-pling regime [18, 23]. Like traditional plasma theories, this is based on a binary scattering approximation, butwhere physics associated with many body correlations isincluded through the use of an effective interaction po-tential. This effective interaction potential was relatedto the potential of mean force, which is the interactionpotential between two particles taking all surroundingparticles to be at fixed positions. Like the other theo-ries previously discussed, this also requires only the pair-distribution as input. Figure 8 shows that this approachis accurate across coupling regimes up to approximatelythe minimum in the viscosity coefficient.Breakdown of the effective potential theory arises atsufficiently strong coupling that the potential compo-nent of the viscosity dominates. This is expected be-cause transport theories based on binary collisions onlyaccount for changes in the particle momenta, so they canat most describe the kinetic contribution. This is shownin detail in Fig. 10. This figure shows the kinetic-kineticand potential-potential terms of the viscosity computedfrom MD using components of J xy ( t ) based on σ ↔ kin and σ ↔ pot = σ ↔ sr + σ ↔ lr . We found that the cross terms (kinetic-potential and potential-kinetic) were negligible across thedomain.For the theoretical evaluation, the viscosity was com-puted from the Chapman-Enskog relation η ∗ = 5 √ π √ / Ξ (2 , (42)where Ξ (2 , was obtained using the method of [18, 23]inputing a pair distribution function calculated from theHNC approximation (no bridge function was included forthe HNC computations used here). Figure 10 shows thatthis theory accurately tracks the kinetic-kinetic term,but contains no information about the potential-potentialterm. This is similar to how binary collision operatorspredict only the ideal gas component of the equation ofstate, whereas an additional term dependent on the pairdistribution is required to describe the potential contri-bution at strong coupling. The effective potential theorybreaks down at sufficiently strong coupling even for trans-port coefficients that do not have potential components,such as diffusion or temperature relaxation rates [18, 23],but the inaccuracy beyond this threshold is not as severefor these coefficients. V. SUMMARY
We have carried out a detailed study of the calcula-tion of the shear viscosity coefficient of one-componentplasmas with equilibrium MD simulation in order to in-dependently validate the non-equilibrium MD results of[5] for κ > −1 −2 −1 Coupling Strength ( Γ ) Sh e a r V i s c o s i t y ( η ∗ ) MD: KineticMD: PotentialMD: TotalBaalrud-DaligaultLandau-Spitzer
FIG. 10: Contributions to the viscosity coefficient computedfrom MD for κ = 0 (for finite κ values, see ?? ): kinetic(squares), potential (triangles), and total (circles). Alsoshown is the prediction of the effective potential theory (dia-monds) and Landau-Spitzer theory (green line). Acknowledgments
This work was carried out under the auspices of theNational Nuclear Security Administration of the U.S. De-partment of Energy (DOE) at Los Alamos National Lab-oratory under Contract No. DE-AC52-06NA25396. Thework of J.D. and K.Ø.R. was supported by the DOEOffice of Fusion Sciences. The work of S.D.B was sup-ported in part by the University of Iowa and in part byLos Alamos National Laboratory. [1] M. Baus and J.-P. Hansen,
Phys. Rep. , 1 (1980).[2] T. Saigo and S. Hamaguchi, Phys. Plasmas , 1210(2002).[3] P. Vieillefosse and J.P. Hansen, Phys. Rev. A , 1106(1975).[4] J. Daligault, Phys. Rev. Lett. , 065003 (2006). TheMD results for the viscosity coefficients published in thispaper are incorrect because of an unintentional mistake ofthe author in implementing the formula for the viscosityin his code.[5] Z. Donk´o and P. Hartmann, Phys. Rev. E , 026408(2008).[6] S. Bastea, Phys. Rev. E , 056405 (2005).[7] B. Bernu and P. Vieillefosse, Phys. Rev. A , 2345(1978).[8] J. Wallenborn and M. Baus, Phys. Rev. A , 1737(1978).[9] G. Salin and J.-M. Caillol, Phys. Plasmas , 1220(2003).[10] R. Hockney and J. Eastwood, Computer Simulation usingParticles (IOP Publishing, 1988).[11] D. Frenkel and B. Smit,
Understanding Molecular Dy-namics (Academic Press, 2002).[12] J.P. Hansen and I.R. McDonald, theory of Simple Liquids (Academic, London, 1986).[13] J. Daligault, unpublished.[14] R. Zwanzig and N.K. Ailawadi,
Phys. Rev. , 280(1969).[15] I. Bitsanis, M. Tirrell and H. Ted Davis,
Phys. Rev. A , 958 (1987).[16] L. Spitzer, Jr., Physics of Fully Ionized Gases, 2nd Ed. (Interscience, New York, 1962).[17] G. Dimonte and J. Daligault, Phys. Rev. Lett. ,135001 (2008).[18] S.D. Baalrud and J. Daligault,
Phys. Rev. Lett. ,235001 (2013).[19] J. Daligault, Phys. Rev. Lett. 108, 225004 (2012).[20] S. Ichimaru,
Statistical Plasma Physics, Vol. I: BasicPrinciples , Addison-Wesley Publ. Company (1992).[21] S. Tanaka and S. Ichimaru,
Phys. Rev. A , 4163 (1986).[22] E. Braun, Phys. of Plasmas , 731 (1967).[23] S.D. Baalrud and J. Daligault, Phys. Plasmas , 055707(2014).[24] T. Saigo and S. Hamaguchi, Phys. Plasmas , 1210(2002). Appendix A: Kinetic-kinetic term
The inital value of the kinetic-kinetic contribution tothe shear stress correlation function is (cid:10) σ kin xy (0) σ kin xy (0) (cid:11) eq = lim t →∞ ¯ J kinxy ( t )where3¯ J kinxy ( t ) = 1 t (cid:90) t N (cid:88) i =1 mv x,i ( s ) (cid:88) j =1 mv y,j ( s ) ds + 1 t (cid:90) t ds N (cid:88) i (cid:54) = j =1 m v x,i ( s ) v x,j ( s ) v y,i ( s ) v y,j ( s ) ds = (cid:34) t (cid:90) t N (cid:88) i =1 mv x,i ( s ) ds (cid:35) (cid:34) t (cid:90) t N (cid:88) i =1 mv y,i ( s ) ds (cid:35) + cross terms( t )and cross terms( t ) = 1 t (cid:90) t (cid:34) N (cid:88) i =1 mv x,i ( s ) − t (cid:90) t N (cid:88) i =1 mv x,i ( s ) ds (cid:35) × (cid:34) N (cid:88) i =1 mv y,i ( s ) − t (cid:90) t N (cid:88) i =1 mv y,i ( s ) ds (cid:35) ds + 1 t (cid:90) t N (cid:88) i (cid:54) = j =1 m v x,i ( s ) v x,j ( s ) v y,i ( s ) v y,j ( s ) ds (A1)In the limit t → ∞ , (cid:34) t (cid:90) t N (cid:88) i =1 mv x,i ( s ) ds (cid:35) (cid:34) t (cid:90) t N (cid:88) i =1 mv y,i ( s ) ds (cid:35) = ( N k B T )( N k B T ) and cross terms( tt