Mean Force Kinetic Theory: a Convergent Kinetic Theory for Weakly and Strongly Coupled Plasmas
aa r X i v : . [ phy s i c s . p l a s m - ph ] A p r Mean Force Kinetic Theory: a Convergent Kinetic Theory for Weakly andStrongly Coupled Plasmas
Scott D. Baalrud and J´erˆome Daligault Department of Physics and Astronomy, University of Iowa, Iowa City, IA 52242 Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545 (Dated: 22 April 2019)
A new closure of the BBGKY hierarchy is developed, which results in a convergent kinetic equation thatprovides a rigorous extension of plasma kinetic theory into the regime of strong Coulomb coupling. The ap-proach is based on a single expansion parameter which enforces that the exact equilibrium limit is maintainedat all orders. Because the expansion parameter does not explicitly depend on the range or the strength of theinteraction potential, the resulting kinetic theory does not suffer from the typical divergences at short andlong length scales encountered when applying the standard kinetic equations to Coulomb interactions. Theapproach demonstrates that particles effectively interact via the potential of mean force and that the rangeof this force determines the size of the collision volume. When applied to a plasma, the collision operator isshown to be related to the effective potential theory [Baalrud and Daligault, Phys. Rev. Lett , 235001(2013)]. The relationship between this and previous kinetic theories is discussed.
I. INTRODUCTION
Derivations of standard plasma kinetic equationsconfront either infrared divergences (Boltzmann equa-tion) , ultraviolet divergences (Lenard-Balescu equa-tion) , or both (Landau equation) . These divergencesare resolved by invoking physical arguments, such as De-bye shielding to resolve the infrared divergence, or thedistance of closest approach in a binary collision to re-solve the ultraviolet divergence. The results obtained ineither approach agree to logarithmic accuracy as long asthe plasma is weakly coupled and near local equilibrium.They have been validated for many collisional transportprocesses satisfying these limits over the several decadessince they were developed. Despite the success of theseequations, the theory stands in an unsatisfactory posi-tion in that it relies on ad hoc arguments. As a result, itis difficult to attempt generalizations to other importantsituations, such as moderate or strong Coulomb coupling,or strong magnetization.Using a new self-consistent closure of the BBGKY hi-erarchy, we present a plasma kinetic equation that (i) isconvergent, i.e. it does not confront either infrared or ul-traviolet divergences, and that (ii) provides a rigorous ex-tension of plasma kinetic theory into the regime of strongCoulomb coupling. The approach demonstrates that par-ticles effectively interact via the potential of mean forceand that the range of this force determines the size of thecollision volume.In order to put our approach into perspective, we firstbriefly recall the origin of standard kinetic equations.Each can be derived from a perturbative closure of theBBGKY hierarchy obtained by first identifying a smalldimensionless parameter characteristic of the system, andby then using its smallness to truncate the hierarchy atthe level of two-particle correlations. Two dimension-less parameters naturally arise in terms of the averageparticle density n , and the range l and the strength φ of the interaction potential φ ( r ). Namely, the “concen- tration parameter” nl , which measures how many par-ticles, on the average, simultaneously interact with agiven particle, and the “strength parameter” φ /k B T ,which measures the interaction potential energy of twocolliding particles in units of the mean particle kineticenergy. It is easily seen that, strictly speaking, theseparameters are not well adapted to plasmas since theCoulomb potential has an infinite range, and its magni-tude becomes arbitrarily strong at short distance. Thecelebrated Boltzmann equation is obtained by assum-ing nl ≪ φ /k B T ∼ O (1); it applies to dilutegases interacting via short range potentials of arbitraryinteraction strength. The Landau equation is obtainedin the so-called weakly coupled limit characterized by φ /k B T ≪ nl ∼ O (1); it applies to weak poten-tials, i.e. potentials that are uniformly small for all in-terparticle distances r , a condition that is generally notrespected due to the typical strong repulsions at shortdistances. The Lenard-Balescu equation, which was de-veloped to deal with weakly coupled plasmas, is obtainedin the so-called weakly-coupled, long-range limit assum-ing φ /k B T ≪ nl φ /k B T ∼ O (1). Plasmas donot fall in any of these categories, which explains whyone needs to regularize ad-hoc the standard kinetic equa-tions before applying them to plasmas. A solution to theproblem can be obtained by writing the Coulomb poten-tial as a sum of two terms: a weak long-range term, plus astrong short-range term. This amounts to carefully join-ing the Lenard-Balescu and Boltzmann equations (see,e.g., Refs. 8–15).Until today, no one has derived a practical kinetic equa-tion beyond the previous limits. The generalization of theBoltzmann equation to higher densities remains an un-solved challenge (much progress on the kinetic theory ofdense gases and liquids has been made, but for systemsin thermal equilibrium). It was shown that the system-atic inclusion of many-body collisions through a densityexpansion similar to the virial expansion for computingthe equilibrium properties of dense gases is plagued byunphysical divergences. This is because N -body colli-sions and correlations cannot be treated separately fromthose of all higher orders ( N +1-body, N +2-body, etc.).This difficulty is symptomatic of many perturbative ap-proaches in physics and could in principle be dealt withusing so-called renormalization techniques. While suchrenormalization schemes have been proposed, their ex-ceeding complexity has thus far restricted their usefulnessto formal, but not practical, solutions.In this paper, we present a closure scheme for theBBGKY hierarchy that, like the standard equations, re-sults from the identification of a small expansion param-eter but, unlike the concentration and strength param-eters, does not explicitly depend on the characteristicsof the interaction potential like l and φ . The presentexpansion parameter [Eq. (11) below] is a measure ofthe perturbation of the distribution function about ther-mal equilibrium. The closure retains the correct equi-librium limit at all orders, and ensures that screeningis captured in the near-equilibrium limit, while at thesame time allowing inclusion of short-range interactions.It results in a convergent and tractable kinetic equationthat provides a rigorous extension of plasma kinetic the-ory into the regime of strong Coulomb coupling. Here,strong Coulomb coupling refers to plasmas in which theaverage Coulomb potential energy of interacting particles( Ze ) / πǫ a ( Ze is the electric charge, a = (3 / πn ) / isthe average interparticle spacing) exceeds their averagekinetic energy k B T ( T is the temperature), i.e. Γ > ∼ Ze ) / πǫ ak B T . (1)The resulting collision operator is the same as what isobtained from the Boltzmann equation if the interparti-cle force in binary collisions is taken to be the potentialof mean force, rather than the Coulomb potential. Inthe weakly coupled limit, the potential of mean force isthe Debye-H¨uckel potential and the theory gives thesame results as standard plasma kinetic equations. Thetransport coefficients resulting from the convergent ki-netic equation are those of the recent “effective poten-tial theory” (EPT), which was shown to successfully ex-tend the conventional plasma transport theory into thestrongly coupled regime.
Previous work has shownthat EPT can extend plasma kinetic theory well into theregime of strong Coulomb coupling.
In addition, in the resulting hydrodynamic equationsthe pressure and internal energy are composed of twoparts; a kinetic (or ideal) part and a potential (or ex-cess) part. The ideal part, which is obtained with allkinetic equations, is the standard ideal gas componentthat represents the transfer of momentum or energy dueto the flow of particles. The excess part, which is typ-ically absent from standard kinetic theories, representsthe transfer of momentum or energy between particlesby the particle interactions.The paper is organized as follows. In Section II, weintroduce the expansion parameter at the basis of the closure of the BBGKY hierarchy. In Section III, the con-vergent kinetic equation resulting from this closure is de-rived and its properties are discussed. To this end, weclosely follow the method used by Grad in his derivationof the Boltzmann equation. In Section IV, the conver-gent kinetic equation is discussed in comparison to thestandard plasma kinetic equations.
II. BASIC EXPANSION PARAMETER
Kinetic theories can be derived from the BBGKY hi-erarchy (cid:20) ∂∂t + n X i =1 (cid:18) L i + L ext i + n X j =1 j = i L Cij (cid:19)(cid:21) f ( n ) ( r n , v n , t ) (2)= − n X i =1 Z d Γ n +1 L Ci,n +1 f ( n +1) ( r n +1 , v n +1 , t ) , where f ( n ) ( r n , v n , t ) = N !( N − n )! Z d Γ ( N − n ) f [ N ] ( r N , v N , t )(3)defines the n th -order reduced distribution functions interms of the N -particle distribution function f [ N ] . Here, r N = ( r , r , . . . r N ), where r i denotes the spatial loca-tion of particle i , v N = ( v , v , . . . v N ) where v i denotesthe velocity of particle i , d Γ n ≡ d r n d v n is a short-hand notation for the 6-dimensional phase-space, and d Γ ( N − n ) ≡ d Γ n +1 . . . d Γ N . The BBGKY hierarchy fol-lows from integrating the Liouville equation (cid:20) ∂∂t + N X i =1 (cid:18) L i + L ext i + N X j =1 j = i L Cij (cid:19)(cid:21) f [ N ] ( r N , v N , t ) = 0 , (4)in order to obtain an evolution equation for each reduceddistribution function. In equation (2), L i = v i · ∂∂ r i + q i m i ( v i × B ) · ∂∂ v i (5)where the second term is associated with the Lorentzforce due to an external magnetic field, and L ext i = 1 m i F ext i · ∂∂ v i (6)where F ext i = −∇ i φ ext ( r i ) is an external force, associ-ated with an external potential. The external potentialis added here for later convenience to introduce spatialinhomogeneities in the thermal equilibrium state. Theterm L Cij ≡ m i F Cij · ∂∂ v i (7)is associated with the electrostatic Coulomb interactionsbetween particles F Cij = q i q j πǫ o ( r i − r j ) | r i − r j | = −∇ i φ ( r ij ) . (8)Any kinetic theory aims to describe the evolution ofthe first-order reduced distribution function f (1) ( r , v , t )( n = 1): (cid:18) ∂∂t + L + L ext1 (cid:19) f (1) (1) = − Z d Γ L C f (2) (1 , , (9)which depends on f (2) (1 ,
2) [(1 ,
2) is shorthand notationfor ( r , v , r , v , t )]. The task of a collision operatoris to provide an approximate expression for f (2) (1 , , t ).For example, this can be related to f (3) (1 , , , t ) via the n = 2 equation (cid:18) ∂∂t + L + L ext1 + L + L ext2 + L C + L C (cid:19) f (2) (1 , − Z d Γ (cid:18) L C + L C (cid:19) f (3) (1 , , , but approximations must be introduced to close the hi-erarchy.The most famous closure is to neglect the third orderdistribution function f (3) → nl ≪ f (2) subject to a constraint on the length scale overwhich binary collisions occur, as well as a lack of initialcorrelations in a binary scattering event. This methodwas shown to lead to the Boltzmann equation. The cal-culation is not trivial and several variations exist; belowwe shall rely on the method proposed by Grad. However,as recalled in the introduction, this closure is a poor ap-proximation for plasmas. The reason is that the physicsof screening is contained in f (3) . Neglecting screeningnot only misses an important physical process, but theresulting kinetic equation diverges because the two-bodyCoulomb force has an infinite range.The closure scheme proposed here enforces that thecorrect equilibrium (i.e., thermodynamic) limit is main-tained at all orders. This ensures that screening is cap-tured, while at the same time allowing inclusion of short-range interactions. This can be accomplished by takingthe basic expansion parameter to be∆ f ( n +1) ≡ f ( n +1) o (cid:18) f ( n +1) f ( n +1) o − f ( n ) f ( n ) o (cid:19) . (11)Equation (11) measures the difference of the n + 1 and n probability distributions, referenced to their equilib-rium values. It is a measure of the perturbation of non-equilibrium correlations about equilibrium; ∆ f ( n +1) → f ( n ) o ( r n , v n ) = ρ ( n ) ( r n ) f ( n )M ( v n ) (12)is the equilibrium reduced distribution function, f ( n )M ( v n ) = (cid:18) m πk B T (cid:19) n/ exp (cid:18) − n X i =1 m v i k B T (cid:19) (13)is the Maxwellian velocity distribution function, ρ ( n ) ( r n ) = N !( N − n )! 1 Z N Z d r ( N − n ) e − ( V ext + V N ) /k B T (14)is the n − particle density distribution function, V N ( r N ) = P Ni =1 P Nj>i φ ( r ij ) is the electrostatic potential en-ergy, V ext ( r N ) = P Ni =1 φ ext ( r i ) is the interaction en-ergy with the external potential φ ext , and Z N = R exp [ − ( V ext + V N ) /k B T ] d r N is the configurational in-tegral.The BBGKY hierarchy from Eq. (2) can be rearrangedin a way that expresses ∆ f ( n +1) as an expansion param-eter (cid:20) ∂∂t + n X i =1 ( L i + ¯ L ( n ) i ) (cid:21) f ( n ) = − n X i =1 Z d Γ n +1 L Ci,n +1 ∆ f ( n +1) , (15)where¯ L ( n ) i ≡ m i ¯ F ( n ) i ( r n ) · ∂∂ v i (16a)= L ext i + n X j =1 j = i L Cij + Z d Γ n +1 L Ci,n +1 f ( n +1) o f ( n ) o (16b)is an operator associated with force¯ F ( n ) i ( r n ) = F ext i + n X j =1 j = i F Cij (17)+ Z F Ci,n +1 ρ ( n +1) ( r n +1 ) ρ ( n ) ( r n ) d r n +1 . The force (17) has a simple physical interpretation. It isthe mean force acting on particle i obtained when keep-ing both it and a set of other particles ( j = 1 , . . . , n ,excluding j = i ) at fixed positions and averaging over allequilibrium configurations of the other N − n particles.As shown in Appendix A, this statistical force can beexpressed as the gradient of a potential,¯ F ( n ) i ( r n ) = −∇ i w ( n ) ( r n ) , (18)where the potential of mean force is w ( n ) ( r n ) = − k B T ln (cid:20) ρ ( n ) ( r n ) ρ n (cid:21) , (19)where ρ = lim V →∞ R V d r ρ (1) ( r ) / V is the average par-ticle density. In the absence of the external potential, φ ext ≡ w ( n ) ( r n ) = − k B T ln g ( n ) ( r n ), where g ( n ) ( r n ) = ρ ( n ) ( r n ) /ρ n is the n -particle distribution function (when φ ext = 0, g ( n ) ( r n ) = ρ ( n ) ( r n ) / Π ni =1 ρ (1) ( r i )).Although Eq. (15) is equivalent to Eq. (2), writing it inthis way affords certain pedagogical clarities. The right-hand side, at any order n , can now be interpreted asa “collision operator” in the sense that it vanishes atequilibrium and is small compared to the left side of theequation for slight perturbations from equilibrium. Thiscontrasts with the right side of Eq. (2), which does notvanish at equilibrium. Of course, solving Eq. (15) stillrequires a closure. The scheme suggested here is that thedynamical evolution of f ( n ) be closed at order n by taking∆ f ( n +1) →
0. However, a closure for the equilibrium dis-tribution ρ ( n ) ( r n ) (or g ( n ) ( r n ) in the homogeneous case)is still required in order to determine the potential ofmean force arising on the left side of the equation. De-termining ρ ( n ) is a more tractable problem because onecan rely on methods of equilibrium statistical mechan-ics. At equilibrium, the BBGKY hierarchy reduces tothe Yvon-Born-Green (YBG) limit ∇ ρ ( n ) ( r n ) − k B T (cid:18) F ext1 + n X j =2 F C ,j (cid:19) ρ ( n ) ( r n ) (20)= ρk B T Z d r n +1 F ,n +1 ρ ( n +1) ( r n +1 ) . Although this is also a hierarchical equation, accurate ap-proximations have been developed for most potentials ofinterest and for quite general conditions of density andtemperature. We also note that Eq. (19) follows di-rectly from Eq. (20). To proceed with the kinetic theoryderivation, ρ ( n ) will be considered a known quantity.In particular, in the following section a kinetic equa-tion is derived from the second order ( n = 2) term ofEq. (15) assuming that ∆ f (3) = 0. The derivation closelyfollows Grad’s method. Many methods have been usedto derive the Boltzmann equation from Eq. (2), andany of these could also be used for our purposes. Wechoose Grad’s method because it follows directly fromthe closure f (3) → f (3) → ρ (2) ( r , r ) will be as-sumed known. We will focus on the near-homogenouslimit that is relevant to fluid theory. Here, the two-particle density consists of a component that varies onlarge scales indicative of the fluid-scale gradients, andanother short spatial scale indicative of the scale of inter-actions: ρ (2) ≈ ρ ( r ) g (2) ( r , r ), where r = | r − r | . Theradial distribution function g (2) ( r ) can be provided bythe hypernetted chain approximation, which is knownto be accurate for plasmas at conditions spanning weakto strong coupling. However, the theory does not dependon the method used to obtain g (2) . At weak coupling, g (2) simply asymptotes to the Debye-H¨uckel limit andconventional plasma kinetic theory will result from thisunified framework. III. PLASMA KINETIC EQUATIONA. Modified Grad Method
The plasma kinetic equation (9) in terms of the expan-sion parameter of Eq. (11) is (cid:18) ∂∂t + L + ¯ L (1)1 (cid:19) f (1) (1) = − Z d Γ L C ∆ f (2) (1 , . (21)Following Grad’s method, space is divided into a volumewithin which collisions occur V σ , and a volume outsideof this in which the observable (truncated) distributionfunctions are defined f ( n ) σ ( r n , v n , t ) ≡ N !( N − n )! Z ∼ V σ d Γ ( N − n ) f [ N ] ( r N , v N , t ) . (22)Here, ∼ V σ denotes that the spatial integral excludes thesmall volume V σ . By integrating Liouville’s equation (4),the first-order BBGKY hierarchy equation for the “trun-cated reduced distribution function” is (cid:18) ∂∂t + L + ¯ L (1)1 (cid:19) f (1) σ (1) = − I S d v d S · u f (2) σ (1 , − Z ∼ V σ d r Z d v L C ∆ f (2) (1 ,
2) (23)where u ≡ v − v . Here, an extra term arises in com-parison to Eq. (21) that depends on the two-body distri-bution function evaluated at the surface of the collisionvolume, S . The method relies on the assumption thatthe collision volume, V σ , is sufficiently small that there isan appropriate limit in which f ( n ) σ → f ( n ) (i.e., that f ( n ) σ is the observable of interest).A key question is, what determines the scale of V σ ? Fora dilute gas, this is typically associated with the rangeof the short-range intermolecular forces. In these theo-ries, f (2) replaces ∆ f (2) in Eq. (23) and the last termis considered negligible because the intermolecular force(inside the L operator) is small in the region of spaceoutside of the collision volume. This term is dropped,and the collision operator is derived by solving the two-body interaction problem inside the collision volume todetermine f (2) on its surface. This leads to the Boltz-mann equation. Clearly, the same argument does not apply to plasmasbecause the Coulomb force is long-range and effectivelyextends over all of space. Indeed, this is why applica-tion of the Boltzmann equation to plasma results in adivergence. However, writing Eq. (23) in terms of ∆ f (2) introduces the notion of screening because the effectiverange of interaction is associated with the spatial cor-relation scale of ∆ f (2) itself. Of course ∆ f (2) vanishesat equilibrium at all spatial scales, so the relevant cor-relation arises only away from equilibrium. Consider-ing a small perturbation from equilibrium, ∆ f (2) will al-ways become small on large scales at which two particlesbecome decorrelated. Only at sufficiently small spatialscales (i.e., the collision scale | r − r | < V σ ) will ∆ f (2) be appreciable in magnitude. For instance, consider ahomogenous plasma. Starting at large scales where spa-tial correlations are small f (2) (1 , ≈ ρ f (1) ( v ) f (2) ( v ),the spatial component of ∆ f (2) will have the possibilityof contributing at the scale associated with spatial corre-lation of f (2) o ; i.e., ∆ f (2) ∼ ρ − ρ (2) ≈ ρ [1 − g ( r )]. Thisscale is directly determined by the potential of mean forcevia Eq. (19), and thus the screening length. For example,in a weakly coupled plasma the potential of mean forceis the Debye-H¨uckel potential w (2)DH ( r ) = φ ( r ) k B T e − r/λ D (24)where λ D is the Debye length. This corresponds directlyto the conventional notion that V σ is a “Debye interac-tion sphere” at weak coupling. However, the potentialof mean force is a more general concept that extends tostrongly coupled plasmas as well.Since V σ is defined as the excluded region of the volumeintegral, the last term in Eq. (23) is negligible and thekinetic equation is (cid:18) ∂∂t + L + ¯ L (1)1 (cid:19) f (1) (1) = − I S d v d S · u f (2) (1 , . (25)Solving for the collision operator in Eq. (25) requires solv-ing for f (2) on the surface of the collision volume. Asdiscussed in the previous section, this will be determinedfrom the second order equation [ n = 2 in Eq. (15)] using∆ f (3) = 0 as a closure (cid:18) ∂∂t + L + L + ¯ L (2)1 + ¯ L (2)2 (cid:19) f (2) (1 ,
2) = 0 . (26)As does the derivation of the Boltzmann equation,Eq. (26) describes the dynamical evolution of two inter-acting particles. However, retention of the higher or-der terms arising in the operators ¯ L ( n ) i shows that theparticles effectively interact via the mean force, ratherthan the bare Coulomb force F Cij . In this way, surround-ing particles “mediate” the binary interactions; a con-cept Rostoker referred to as “dielectric dressing”. Thus,although the dynamics are two-body, the interaction ismore general than that used in the Boltzmann equation.In this sense, it relaxes the binary collision assumptionto a certain degree by including the statistical influenceof all other N − B. Collision Operator
The vast majority of plasmas are magnetized in thesense that the Lorentz force may influence macroscopic dynamics, but it does not influence microscopic dy-namics within the collision volume. In a weakly cou-pled plasma, this condition is satisfied if the gyroradiusof most particles is signficantly larger than the Debyelength: r c ≫ λ D , where r c = v T /ω c is the thermal gy-roradius and ω c = | Ze | B/m is the gyrofrequency. Ina strongly coupled plasma, it requires r c ≫ a . Mostof plasma kinetic theory is based on this scale separa-tion, and this section concentrates on this situation. Inthis case, the Lorentz force term in L contributes to theplasma kinetic equation (21), but it is negligible in theoperators L and L in Eq. (26), which describes f (2) inside the collision volume.Following the method of Ref. 1, we choose a coordi-nate system such that a diametral plane intersecting thesphere is aligned perpendicular to the relative velocityvector u . Using polar coordinates in this plane, an areaelement is denoted dω = rdrdǫ . This transformationmaps the projection of S onto the disk 0 ≤ r ≤ σ and0 ≤ ǫ ≤ π . The disk ω is covered twice: once by thehemisphere S +2 : u · d S > S − : u · d S < S +2 represent particles moving awayfrom one another (post-collision particles), while thosethat map onto S − represent particles moving towardone another (pre-collision particles). Denoting the pointsthat map onto S +2 as r +2 and those that map onto S − as r − , we find u · d S = − udω on S +2 and u · d S = udω on S − . In terms of this transformation, Eq. (25) becomes (cid:18) ∂∂t + L + ¯ L (1)1 (cid:19) f (1) (1) (27)= Z dωd v u [ f (2) ( r , v , r +2 , v ) − f (2) ( r , v , r − , v )] . To obtain an explicit collision operator, we invoke fourassumptions that are similar to those underlying theBoltzmann equation: (1) truncation (finite size of V σ ),(2) binary collisions, (3) modified molecular chaos, and(4) slow variation of f (1) . Assumption (1) was alreadyapplied above, which justified the neglect of the last termin Eq. (23), and that f (1) σ and f (2) σ accurately approxi-mate f (1) and f (2) . We point out again that assumption(1) has a different origin in the context of the theory pre-sented here, in comparison to the Boltzmann equation.Here, the small size of V σ is justified only by the limitedrange of the effective interaction implied by ∆ f (2) . Incontrast, in terms of f (2) the neglected term in Eq. (23)would diverge for a plasma. Neglecting this is justifiedin derivations of the Boltzmann equation because theyfocus on dilute gases which have short-range potentials.Since the collision volume is small, the number ofparticle inside of it is very small compared to the to-tal number of particles. Furthermore, we assume thatwithin the collision volume particle interactions are pre-dominately binary [assumption (2)], but that they oc-cur via the mean force. This justifies ∆ f (3) = 0 asa closure, leading to Eq. (26) within V σ . The solu-tion of this equation is f (2) = const on a 2-particletrajectory. Thus, the points that map into the hemi-sphere S +2 can be equated with post-collision coordinates: f (2) ( r , v , r +2 , v , t ) → f (2) (ˆ r , ˆ v , ˆ r , ˆ v , ˆ t ).To solve for the particle trajectories inside the sphere,the Boltzmann equation typically makes the “molecularchaos” approximation (3) whereby the two-particle dis-tribution is assumed to be uncorrelated in both the ini-tial and final states: f (2) (1 ,
2) = f (1) (1) f (1) (2). This isa rather more subtle point associated with the introduc-tion of irreversibility into the theory, which is discussedin depth in Ref. 1. The procedure presented here at-tempts to develop a model that is consistent with theequilibrium state. As such, we modify this approxima-tion slightly to account for the statistical spatial corre-lation of particles at the surface of the collision volume: f (2) (1 , | r = σ ≈ χf (1) (1) f (1) (2). This concept, first in-troduced by Enskog in the context of a hard-sphere gas,accounts for the “excluded volume” associated with thefact that hard spheres cannot overlap. For hard spheres, χ = g ( r = σ ). For a soft potential, σ is not a singlevalue for all interactions, and as such the concept is un-derstood as a statistical analog. This is discussed inRef. 22 for single-component plasmas. In this model, σ is associated with the location where g ( r = σ ) = 0 . χ is determined by using this effective diameter inthe virial expansion of the Enskog equation of state; seealso Ref. 28 for an extension to binary mixtures. Thevalue 0 .
87 comes from a fit to molecular dynamics simu-lations of the OCP, but is expected to be a universal valuefor repulsive potentials.
From this model χ → < ∼
1, and χ ≈ . − . < ∼ Γ < ∼ f (1) . In this limit, we expectthe spatial coordinates ˆ r , r and ˆ r to differ from r bya small amount on the order of σ . Likewise, the shortcollision time implies that ˆ t differs from t by a negligibleamount. Putting the results of these four approximationsinto Eq. (27) leads to the plasma kinetic equation (cid:18) ∂∂t + L + ¯ L (1)1 (cid:19) f (1) ( v ) (28)= χ Z dωd v u [ f (1) (ˆ v ) f (1) (ˆ v ) − f (1) ( v ) f (1) ( v )] , where it is understood that each occurrence of f (1) isevaluated at the same local position r and time t . Equa-tion (28) was derived for a one-component plasma. Thegeneralization to a multicomponent system (cid:18) ∂∂t + L + ¯ L (1)1 (cid:19) f (1) s ( v ) (29)= X s ′ χ ss ′ Z dωd v u [ f (1) s (ˆ v ) f (1) s ′ (ˆ v ) − f (1) s ( v ) f (1) s ′ ( v )] , follows directly by tracking the derivation above and dis-tinguishing species s from s ′ . Since the mean force is central and conservative, thedifferential area of the disk can be related to a scatteringcross section, dω = σ ss ′ d Ω where d Ω = dφdθ sin θ , or animpact parameter, dω = bdbdφ , in the usual way. Theresults of the standard two-body scattering problem aredetermined by the scattering angle θ = π − b Z ∞ r o dr r − (cid:20) − b r − w (2) ss ′ ( r ) m ss ′ u (cid:21) − / . (30)Here, m ss ′ = m s m s ′ / ( m s + m s ′ ) is the reduced mass, and r o is the distance of closest approach, which is determinedby the largest root of the term in square brackets. Notethat this is the same as the scattering angle used in theBoltzmann equation, except that the potential of meanforce w (2) replaces the bare interaction potential φ ( r ). C. Comment on Strongly Magnetized Plasmas
If the external magnetic field is sufficiently strongthat the gyromotion of particles fits within the col-lision volume, the situation is much more compli-cated. O’Neil considered very strongly magnetized one-component plasmas (such that r c ≪ r o ), and derived akinetic equation by aligning the surface d S along themagnetic field (rather than u ). A collision operatorwas then derived based on an assumption that the adi-abatic invariant | u + | /B is preserved during the inter-action. Here u + ≡ u x + iu y are Cartesian componentsof u and the magnetic field is in the ˆ z direction. Thisadiabatic assumption is restrictive, and a more generaltheory is desirable.The method outlined above may provide a fruitfulstarting point for developing such a theory. The main ad-vantage is that it self-consistently accounts for screening.A matter of continuing investigation in strongly mag-netized plasmas is the relative importance of short andlong range collisions. Early theories modified Landau-Spitzer based approaches by changing the Debye lengthto the gyroradius in the Coulomb logarithm.
Oth-ers modified the Lenard-Balescu equation to accountfor magnetization in the plasma dielectric function.
However, recent molecular dynamics results have showntrends inconsistent with either prediction at strong mag-netization. By accounting for screening via the poten-tial of mean force, and for gyromotion via the 2-bodyinteraction inside the collision sphere, one may addressin a self-contained way the combined influence of shortand long-range interactions. The Bohr-van Leeuwen the-orem ensures that the magnetic field does not influenceany statistical property at equilibrium, so the potentialof mean force is expected to remain unchanged from thatpresented above.
D. Equation of State
In addition to providing a convergent collision opera-tor, Eq. (25) introduces a term in the convective deriva-tive ( ¯ L (1)1 ) that is absent from the Boltzmann equation.This term is associated with non-ideal components of theequation of state, i.e. the excess (or potential) compo-nent of the pressure and internal energy. Although thisis negligible in a weakly coupled plasma, it is the domi-nant component in moderately and strongly coupled plas-mas. The fact that these terms arise naturally shouldbe expected because the closure is designed to ensure thatthe exact equilibrium state consistent with the YBG hi-erarchy, Eq. (20), is maintained.To make the connection with pressure, it is first usefulto notice that the mean force acting on one particle canbe written as ¯ F (1)1 = − ∇ · P Φ ρ (1) ( r ) (31)where P Φ = − Z d r rr φ ′ ( r ) r Z dµ ρ (2) ( r − (1 − µ ) r , r + µ r )(32)is a second-rank tensor, r ≡ r − r and φ ′ ( r ) ≡ dφ/dr .A derivation of Eq. (31) from (17) is provided in Ap-pendix B. In the common limit that the plasma is suffi-ciently homogeneous that ρ (1) ( r ) is constant on the spa-tial scale of the two-body correlations, then ρ (2) ( r , r ) ≈ [ ρ (1) ( r )] g ( r , r ); see Pohl, et al for a related discus-sion. In this case, the divergence of the pressure tensor inEq. (31) reduces to the gradient of a scalar, ∇ p Φ , where p Φ = − ρ Z ∞ d r φ ′ ( r ) rg ( r ) (33)is the well-known equilibrium expression for the potentialcomponent of the scalar pressure. We note that Eq. (33) makes use of the assumptionof “local thermodynamic equilibrium” whereby there isa scale separation between large (fluid-scale) gradientsrepresented by r and the interaction scale representedby r = | r − r | . Thus, the force represented by ∇ p Φ arises due to a slight inhomogeneity of the backgrounddensity across the collision volume ( L ≫ σ ). This originof the excess pressure can be compared with Enskog’stheory of hard sphere gases. Enskog’s theory expands f (1) ( r ± σ s k ) about the local position r , where σ s is thesphere diameter and k is vector connecting the center oftwo spheres at the instant of contact. The excess pressure(and excess internal energy) arising from the lowest-orderterm in this expansion comes about due to any inhomo-geneity of f (1) ( r ) across the collision volume, which issimply the volume of the sphere in a hard-sphere gas.The higher order terms proportional to ∇ f (1) , etc., areassociated with the non-local nature of the collision itself(a hard sphere gas is a singular case in that σ s is both the maximum and minimum scale of an interaction). Equa-tion (33) represents a generalization of the first of thesecontributions (local collisions taking place in a finite sizedcollision volume) to the case of arbitrary potential. Itsexistence in the theory is enforced by the requirementthat the exact equilibrium, or local equilibrium, limit bemaintained.A closed fluid description follows from a Chapman-Enskog solution of the kinetic equation, but the basicfeatures of non-ideality can be illustrated directly fromthe conservation equations. In particular, the densitymoment of Eq. (29) leads to the usual continuity equation ∂ρ∂t + ∇ · ( ρ V ) = 0 (34)where V ≡ R d v v f (1) /ρ is the hydrodynamic velocityand ρ = ρ (1) ( r ). This is the same as what results fromthe Boltzmann equation. However, the momentum mo-ment leads to ρ m d V dt = −∇ · P + ρ q E ′ + J × B (35)in which the total pressure tensor now includes both ki-netic (ideal) and potential (excess) components P = P K + P Φ . (36)Here, d/dt = ∂/∂t + V · ∇ , P K = P s m s R d v v r v r f (1) s /ρ is the kinetic part of the pressure tensor, v r = v − V , ρ m = P s m s ρ s is the mass density, ρ q = P s q s ρ s thecharge density, J = P s q s ρ s ( V s − V ) the current density,and E ′ = E + V × B the electric field in the referenceframe of the fluid.The mean force will also contribute to the energy bal-ance. Taking the m s v r moment of the kinetic equation,and summing over species, gives ρ du K dt = −P K : ∇ V − ∇ · q K + J · E ′ − V · ( ∇ · P Φ ) , (37)where u K = P s m s R d v v r f (1) s is the internal kineticenergy, and q = P s m s R d v v r v r f (1) s is the kineticcomponent of the heat flux. Noting that V · ( ∇ · P Φ ) = ∇ · ( V · P Φ ) − P Φ : ∇ V , Eq. (37) can alternatively bewritten ρ du K dt = −P : ∇ V − ∇ · q + J · E ′ , (38)where q = q K + q Φ , and q Φ = V · P Φ (39)is a contribution to the heat flux from the mean forcethat is absent in the Boltzmann equation.An evolution equation for the total internal energymust also include the contribution from the average in-terparticle potential energy ρu Φ = 12 Z φ ( r ) f (2) (1 , t ) d v d v d r | r = r . (40)The evolution of u Φ is described by the φ ( r ) moment ofthe second order equation of the BBGKY hierarchy ρ du Φ dt = − ∇ r · Z φ ( r ) v r f (2) (1 , t ) d v d v d r | r = r − Z u · ( ∇ r φ ) f (2) (1 , t ) d v d v d r | r = r . (41)Explicit evaluation of this requires a time-dependent so-lution for f (2) . However, if we apply a Kirkwood-like su-perposition approximation at second order, f (2) (1 , ≈ ρ (2) ( r , r ) f (1) ( v ) f (2) ( v ), as in Sec. III B above,Eq (40) reduces to ρu Φ = R φ ( r ) ρ (2) ( r , r ) d r | r = r ,and Eq. (41) to du Φ /dt = 0. In this limit, the totalenergy evolution equation then has the form ρ dudt = −P : ∇ V − ∇ · q + J · E ′ , (42)where u = u K + u Φ is the total internal energy. Note thatbecause of the non-ideal terms, thermodynamic relationsmust be invoked to relate this to temperature. In thelimit that the plasma is locally homogeneous, the familiarexpression for the potential energy density u Φ = ρ Z d r φ ( r ) g ( r ) (43)is obtained. IV. COMPARISON WITH PREVIOUS MODELS
Common plasma theories can also be derived from clo-sures of the BBGKY hierarchy. Here, these closures arecompared with the one presented above. Specific atten-tion is paid to the approximations that lead to diver-gences by considering the equilibrium limit of the secondorder distribution function resulting from each closure.A common first step in the derivation of these equa-tions is to define a cluster expansion of the form f (1) (1) = f (1) (44a) f (2) (1 ,
2) = f (1) f (2) + g (1 ,
2) (44b) f (3) (1 , ,
3) = f (1) f (2) f (3) + f (1) g (2 ,
3) (44c)+ f (2) g (1 ,
3) + f (3) g (1 ,
2) + g (1 , , . The cluster expansion allows aspects of the triplet distri-bution function to be retained via products of f (1) and g . Applying this, Eq. (10) is (cid:18) ∂∂t + L + L (cid:19) g (1 ,
2) + (cid:0) L C + L C (cid:1) f (1) f (2) (45a)= − (cid:0) L C + L C (cid:1) g (1 ,
2) (45b) − Z d Γ (cid:2) L C g (2 , f (1) + L C g (1 , f (2) (cid:3) (45c) − Z d Γ (cid:0) L C + L C (cid:1) g (1 , , . (45d) Each of the theories discussed below are based on a bi-nary collision approximation in which the triplet corre-lation term Eq. (45d) is neglected, but where each keepsa different subset of terms (45a)–(45c) selected, as illus-trated below for the Landau and Boltzmann equations,by identifying the leading order in the expansion param-eters discussed in the introduction. Divergences arise in each, not because of the neglect oftriplet correlations, but rather because of neglecting oneor more of the terms (45b)–(45c). Although the collisionoperator vanishes at equilibrium due to the velocity de-pendence of the Maxwellian distribution, the divergencesare associated with the spatial component of the distri-bution function. Thus, they can be illustrated by consid-ering just the spatial dependence of the collision operatorat equilibrium C ∝ Z d r L C g eq2 (1 , ∝ Z dr [1 − g ( r )] (46)where g eq2 (1 ,
2) = ρ f M ( v ) f M ( v )[ g ( r ) − . (47) Landau Equation : As recalled in the introduction, theLandau equation is obtained in the weakly coupled limitcharacterized by λ = φ /k B T ≪
1. In terms of thisperturbative parameter, the Liouville operators satisfy L i = O ( λ ) and L ij = O ( λ ), while the distribution func-tions f = O ( λ ) (since it carries the complete normaliza-tion regardless of λ ) and it is assumed g n = O ( λ n − )(e.g., at least one direct interaction is required to cre-ate a two-body correlation from an uncorrelated state).Using these orderings in Eq. (45) and keeping the lowestorder in λ , Landau’s seminal kinetic equation can be de-rived from Eq. (45) by keeping only term (45a) to model g . Solving Eq. (45a) at equilibrium implies g ( r ) = 1 − φ ( r ) k B T . (48)Using this result in Eq. (46) shows that C L ∝ Z dr φ ( r ) k B T ∝ Z dr r , (49)which diverges logarithmically in both the large andsmall r limits. This indicates that the Landau closure ne-glects both short range physics [term (45b)] and screening[term (45c)]. Analogously, the same logarithmic diver-gences would be observed in the terms (45b) and (45c),but these are dropped from the analysis. The same resultis obtained for other plasma kinetic equations that areequivalent to the Landau equation, such as Rosenbluth’sFokker-Planck equation. Landau argued that the limitsof integration should range from the thermal distance ofclosest approach in a binary collision, r L = e /k B T , tothe Debye length λ D . Boltzmann Equation : The Boltzmann equation as-sumes λ = nl ≪
1. In terms of this perturbative pa-rameter, the Liouville operators satisfy L i = O ( λ ) and L ij = O ( λ ), and, from their dependence on the par-ticle density, the distribution functions f = O ( λ ) and g n = O ( λ n ). Using these orderings in Eq. (45) and keep-ing the lowest order in λ , the Boltzmann equation can bederived from Eq. (45) by keeping terms (45a) and (45b),while dropping terms (45c) and (45d) (this is equivalentto taking a closure f (3) = 0 and applying the method ofSec. III). Substituting Eq. (47) into the equation result-ing from this approximation shows that it implies g ( r ) = e − φ ( r ) /k B T (50)at equilibrium. Using this result in Eq. (46) shows that C B ∝ Z dr (cid:20) exp (cid:18) − φ ( r ) k B T (cid:19) − (cid:21) . (51)This expression converges in the close interaction limit( r → λ D ad hoc . Lenard-Balescu Equation : The Lenard-Balescu equa-tion can be derived from Eq. (45) by keeping terms(45a) and (45c) [while dropping terms (45b) and (45d)]. Substituting Eq. (47) into the resultant approximationshows that this implies g ( r ) = 1 − φ ( r ) k B T e − r/λ D . (52)Using this result in Eq. (46) shows that C LB ∝ Z dr φ ( r ) k B T e − r/λ D . (53)This expression converges in the far interaction limit( r → ∞ ), but since the integrand is proportional to 1 /r for small r , it diverges logarithmically in the close in-teraction limit. The neglected term (45b) contains thephysics of the close interaction. This limit is typically re-solved by truncating the spatial integral at the thermaldistance of closest approach in a binary collision, r L . Frieman-Book Equation : Self-consistently convergentkinetic equations have previously been developed. Frie-man and Book derived such an equation by matchingsolutions of Eq. (45) in asymptotic limits of close inter-action ( r < ∼ r L ), far interaction ( r > ∼ λ D ) and interme-diate scale ( r ∼ n − / ). The solution in each limit isobtained by neglecting the terms described above for theBoltzmann equation, Lenard-Balescu equation and Lan-dau equation, respectively. The resultant expression for g eq2 consists of the sum of the Boltzmann equation plusthe Lenard-Balescu equation minus the Landau equation.The radial distribution function then has the form g ( r ) = e − φ/k B T − φk B T e − r/λ D + φk B T (54) (see also Refs. 46 and 47 for further discussion). Usingthis result in Eqs. (46) and (47) shows that C FB ∝ Z dr (cid:18) − e − φ/k B T − φk B T e − r/λ D + φk B T (cid:19) . (55)This expression converges in both the close and far in-teraction limits. This shows that the divergences inthe other theories are associated with neglecting one ofthe terms (45b)-(45c), rather than the triplet correlation(45d). Relation with the theory of Sec. III : Since this proce-dure preserves the exact equilibrium limit, the relation-ship implied by the second order equation is simply theexact statistical statement g ( r ) = exp (cid:20) − w (2) ( r ) k B T (cid:21) . (56)However, a closure of Eq. (20) is still required to deter-mine g ( r ). Since this is an equilibrium quantity, the toolsof equilibrium statistical mechanics are available for thistask, and accurate approximations have been developed.In principle, any such approximation for g ( r ), or evenexperimental data, can serve as the input to this theory.One example of an accurate closure for Coulomb sys-tems is the combination of the hypernetted chain andOrnstein-Zernike equations g ( r ) = exp[ − φ ( r ) /k B T + h ( r ) − c ( r )] (57a)ˆ h ( k ) = ˆ c ( k )[1 + n ˆ h ( k )] (57b)where h ( r ) ≡ g ( r ) −
1, and ˆ h ( k ) denotes the Fouriertransform of h ( r ). This closed set of equations is knownto accurately describe g ( r ) in plasmas spanning fromasymptotically weak coupling well into the strong cou-pling regime. Further extensions have been provided,via models for the bridge function that extend this evenfurther to near solid-state plasma conditions. Thus,very accurate methods are available to determine g ( r ) atequilibrium at essentially any conditions.At weak coupling, Eq. (57) reduces to the Debye-H¨uckel limit g ( r ) = exp (cid:20) − φ ( r ) k B T e − r/λ D (cid:21) . (58)O’Neil and Rostoker showed that this is valid to orderΓ ln Γ for Γ ≪ None of the approximate theoriesdescribed above capture Eq. (58). However, the conver-gent Frieman-Book result is consistent with it to orderΓ ; see Refs. 46 and 47 for further discussion.Other previous kinetic theories, including Liboff andPaquette, have modified the Boltzmann equation inan ad-hoc manner by modeling binary collisions as oc-curring via the Debye-H¨uckel potential (rather than theCoulomb potential). This leads to a convergent kineticequation that accurately models weakly coupled plasmas.0Indeed, if φ ( r ) is replaced by the Debye-H¨uckel poten-tial in Eq. (50) the correct result for g at equilibrium[Eq. (58)] is obtained.In a sense, these theories, as well as Lenard-Balescutheory, can be understood as variants of the “dressedtest particle” concept. Differences lie in the details ofhow interactions are modeled. Lenard-Balescu type ap-proaches model interactions via the correlation of lin-ear fluctuations. This loses the “particle” concept, andthus does not include close collisions, but it does ac-count for dynamics in the dielectric “dressing”, which canbe important for super-thermal particles. For instance,the dynamic aspect of the dielectric response results inorder-unity contributions to fast particle stopping in aplasma. In contrast, the modified Boltzmann based ap-proaches retain the particle concept, and physics of closecollisions. However, the dielectric response in this caseis input in an ad hoc manner, via the screened Coulombpotential, which is accurate for sub-thermal particles and thus expected to apply to near equilibrium trans-port processes.The approach of Sec. III essentially formalizes the lat-ter line of reasoning by enforcing the viewpoint that theexact equilibrium limit should be maintained in the clo-sure to the BBGKY hierarchy. Since this is the Debye-H¨uckel result at weak coupling, the theory reduces tothe result of Liboff in that limit. It also extends thisline of reasoning by revealing that the appropriate in-teraction potential is the potential of mean force, whichdiffers substantially from the Debye-H¨uckel potential atstrong Coulomb coupling. Recent work has shown thatmodeling the potential of mean force via Eq. (57) en-ables an extension of traditional plasma transport the-ory well into the strongly coupled regime (typically forΓ < ∼ Similar results have also been obtained innon-equilibrium two-temperature plasmas.
V. CONCLUSIONS AND OUTLOOK
The closure presented in this paper makes two ad-vances. First, it provides a formal expansion param-eter [Eq. (11)] for the BBGKY hierarchy that enablesa self-contained derivation of a plasma kinetic equationthat includes both static screening and close interactions.Second, it provides a means to extend traditional plasmatheory to strong Coulomb coupling, not only in the col-lision operator, but also in the equation of state. It doesso by ensuring that the equilibrium limit is preserved inthe spatial correlation at second order of the hierarchy.This formalizes and extends the EPT collision oper-ator obtained previously from physical arguments.
It also provides a conceptual basis that may prove use-ful for extending the theory to other regimes of plasmaphysics. For example, kinetic theories have been devel-oped to include three-body interactions.
A similaranalysis could be applied while using ∆ f (4) = 0 to closethe hierarchy, thereby including the static screening re- sponse in effective three-body dynamics. It may providea basis for addressing strongly magnetized plasmas inwhich the Lorentz force acts at the collision scale. Fi-nally, the approach was motivated by plasmas, whichemphasize the need to account for many-body screen-ing self-consistently. However, the theory may also beapplied to other types of systems, particularly in the ki-netic theory of dense gases. Appendix A: Derivation of Equation (18)
Here, it is shown that the mean force given by Eqs. (18)and (19) is equivalent to (17). First, note that for i ≤ n , −∇ i (cid:2) V ext ( r N ) + V N ( r N ) (cid:3) = F ext i + n X j =1 j = i F Cij + N X j = n +1 F Cij so that the gradient of the n − particle density ρ ( n ) ( r n )defined by Eq. (14) can be written k B T ∇ i ρ ( n ) ( r n )= N !( N − n )! 1 Z N Z dr N − n e − ( V ext + V N ) /k B T × F ext i ( r i ) + n X j =1 j = i F Cij + N X j = n +1 F Cij = F ext i ( r i ) + n X j =1 j = i F Cij ρ ( n ) ( r n )+ N !( N − n )! 1 Z N Z dr N − n N X j = n +1 F Cij e − ( V ext + V N ) /k B T . Using, Z dr N − n N X j = n +1 F Cij e − ( V ext + V N ) /k B T = ( N − n ) Z dr N − n F Ci,n +1 e − ( V ext + V N ) /k B T = ( N − n ) Z N ( N − ( n + 1))! N ! Z dr N − n F Ci,n +1 ρ ( n +1) = Z N ( N − n )! N ! Z dr N − n F Ci,n +1 ρ ( n +1) , we find k B T ∇ i ρ ( n ) ( r n ) /ρ ( n ) ( r n )= F ext i ( r i ) + n X j =1 j = i F Cij + Z dr N − n F Ci,n +1 ρ ( n +1) ρ ( n ) , which completes the proof.1 Appendix B: Derivation of Equation (31)
Equation (31) follows from the n = 1 component ofEq. (17) ¯ F (1)1 = − Z d r [ ∇ φ ( r )] ρ (2) ( r , r ) ρ (1) ( r ) . (B1)This derivation follows an analogous method presentedin Ref. 3. First, note that Z d r ∇ φ ( r ) ρ (2) ( r , r ) = − Z d r φ ′ ( r ) r r ρ (2) ( r , r + r )(B2a)= − Z d r φ ′ ( r ) r r ρ (2) ( r + r , r )(B2b)where the second line follows from the property that ρ (2) is constant under the interchange 1 ↔
2. Combiningthese two equivalent expressions, and applying the sub-stitution r → − r as the integration variable in Eq. (B2b),gives Z d r ∇ φ ( r ) ρ (2) ( r , r ) = (B3) − Z d r φ ′ ( r ) r r [ ρ (2) ( r , r + r ) − ρ (2) ( r − r , r )] . Next, observing that ρ (2) ( r , r + r ) − ρ (2) ( r − r , r )= Z dµ ∂∂µ ρ (2) ( r − (1 − µ ) r , r + µ r ) (B4)= Z dµ r · ∇ ρ (2) ( r − (1 − µ ) r , r + µ r ) (B5)shows that Eq. (B3) can be written as the divergence ofa tensor Z d r [ ∇ φ ( r )] ρ (2) ( r , r ) = (B6) − ∇ · Z d r rr φ ′ ( r ) r Z dµρ (2) ( r − (1 − µ ) r , r + µ r ) . Putting Eq. (B6) into (B1) completes the derivation ofEq. (31).
ACKNOWLEDGMENTS
The authors thank Nathaniel Shaffer and Louis Josefor helpful discussions. This work was supported by theU. S. Air Force Office of Scientific Research under AwardNo. FA9550-16-1-0221; and by the U. S. Department ofEnergy, Office of Fusion Energy Sciences, under AwardNo. de-sc0016159. The work of J. D. was supported bythe US Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratoryis operated by Triad National Security, LLC, for the Na-tional Nuclear Security Administration of U. S. Depart-ment of Energy (Contract No. 89233218CNA000001). H. Grad, in
Handbuch der Physik , vol 12 (Springer-Verlag,Berlin, 1958). S. Harris,
An Introduction to the Boltzmann Equation (Dover,New York, 1971). J. H. Ferziger and H. G. Kaper,
Mathematical theory of transportprocesses in gases , (Elsevier, New York, 1972). A. Lenard, Ann. Phys. , 390 (1960). R. Balescu, Phys. Fluids , 52 (1960). R. L. Guernsey, thesis, University of Michigan (1960); Phys. Flu-ids , 1600 (1964). L. Landau, Physik. Z. Sowjetunion , 154 (1936). W. Thompson and J. Hubbard, Rev. Mod. Phys. , 714 (1960). J. Hubbard, Proc. Roy. Soc. (London)
A261 , 371 (1961). D. Baldwin, Phys. Fluids , 1523 (1962). E. A. Frieman and D. L. Book, Phys. Fluids , 1700 (1963). T. Kihara and O. Aono, J. Phys. Soc. Japan , 837 (1963);T. Kihara, ibid . , 108 (1964). J. Weinstock, Phys. Rev. , A673 (1964). O. Aono, J. Phys. Soc. Japan , 1250 (1965); Phys. Fluids ,341 (1968). H. A. Gould and H. E. DeWitt, Phys. Rev. , 68 (1967). J. R. Dorfman and E. G. D. Cohen, J. Math. Phys. , 282 (1967). J. Daligault, J. Stat. Phys. , 1189 (2011). T. L. Hill,
An Introduction to Statistical Thermodynamics (Addison-Wesley, Reading, 1960) p. 313. P. Debye and E. H¨uckel, Phys. Z. , 185 (1923). S. D. Baalrud and J. Daligault, Phys. Plasmas , 055707 (2014). S. D. Baalrud and J. Daligault, Phys. Rev. Lett. , 235001(2013). S. D. Baalrud and J. Daligault, Phys. Rev. E , 063107 (2015). J. Daligault, K. O. Rasmussen, and S. D. Baalrud, Phys. Rev. E , 033105 (2014). M. V. Beznogov, and D. G. Yakovlev, Phys. Rev. E , 033102(2014). T. Haxhimali, R. E. Rudd, W. H. Cabot, and F. R. Graziani,Phys. Rev. E , 023104 (2014). T. S. Strickler, T. K. Langin, P. McQuillen, J. Daligault, andT. C. Killian, Phys. Rev. X , 021021 (2016). J. Daligault, S. D. Baalrud, C. E. Starrett, D. Saumon, andT. Sjostrom, Phys. Rev. Lett. , 075002 (2016). N. R. Shaffer, S. D. Baalrud and J. Daligault, Phys. Rev. E ,013206 (2017). J.-P. Hansen and I. R. McDonald,
Theory of Simple Liquids, 3rdEdition (Academic Press, Oxford, 2006). E. G. D. Cohen, Physica , 163 (1961). N. Rostoker, Phys. Fluids , 479 (1964). S. D. Baalrud and J. Daligault, Phys. Rev. E , 043202 (2017). D. Enskog, Kungl. Svenska Vet.-Ak. Handl. , 1 (1922). H. J. M. Hanley, R. D. McCarty and E. G. D. Cohen, Physica , 322 (1972). T. M. O’Neil, Phys. Fluids , 2128 (1983). D. H. E. Dubin, Phys. Plasmas , 052108 (2014). V. P. Silin, Sov. Phys. JETP , 1281 (1963). D. Montgomery, G. Joyce, and L. Turner, Phys. Fluids , 2201(1974). N. Rostoker, Phys. Fluids , 922 (1960). M.H.A. Hassan, Plasma Physics , 1043 (1977). T. Pohl, T. Pattard, and J. M. Rost, Phys. Rev. A , 033416(2004). J. G. Kirkwood, J. Chem. Phys. , 180 (1946). R. Balescu,
Equilibrium and Nonequilibrium Statistical Mechan-ics (Wiley-Interscience, 1975). M. N. Rosenbluth, W. M. MacDonald, and D. L. Judd,Phys. Rev. , 1 (1957). D. Nicholson,
Introduction to Plasma Theory (John Wiley &Sons, New York, 1983). F. Shure, Phys. Rev. Lett. , 353 (1964). H. E. DeWitt, Phys. Rev. , A466 (1965). S. Ichimaru,
Statistical Plasma Physics, Vol. 1: Basic Principles (Addison-Wesley, Boston, 1992). S. Tanaka and S. Ichimaru, Phys. Rev. A , 4163 (1986). T. O’Neil and N. Rostoker, Phys. Fluids , 1109 (1965). R. L. Liboff, Phys. Fluids , 40 (1959). C. Paquette, C. Pelletier, G. Fontaine and G. Michaud, Astro-phys. Journal Suppl. Series , 177 (1986). G. Zwicknagel, C. Toepffer, P.-G. Reinhard, Phys. Reports ,117 (1999). P. Seuferling, J. Vogel, and C. Toepffer, Phys. Rev. A , 1859(1986). N. R. Shaffer, S. K. Tiwari, and S. D. Baalrud, Phys. Plasmas , 092703 (2017). N. R. Shaffer, and S. D. Baalrud, Phys. Plasmas (2019). S. T. Choh and G. E. Uhlenbeck,