Mean Force Kinetic Theory: a Convergent Kinetic Theory for Weakly and Strongly Coupled Plasmas

TL;DR

This paper introduces a convergent kinetic theory for plasmas that remains valid from weak to strong Coulomb coupling, avoiding divergences of standard theories by using a novel closure based on the potential of mean force.

Contribution

A new closure of the BBGKY hierarchy is developed, providing a convergent kinetic equation applicable to strongly coupled plasmas, extending traditional plasma kinetic theory.

Findings

The kinetic theory remains valid at strong Coulomb coupling.

Particles interact via the potential of mean force, defining the collision volume.

The collision operator relates to the effective potential theory.

Abstract

A new closure of the BBGKY hierarchy is developed, which results in a convergent kinetic equation that provides a rigorous extension of plasma kinetic theory into the regime of strong Coulomb coupling. The approach is based on a single expansion parameter which enforces that the exact equilibrium limit is maintained at all orders. Because the expansion parameter does not explicitly depend on the range or the strength of the interaction potential, the resulting kinetic theory does not suffer from the typical divergences at short and long length scales encountered when applying the standard kinetic equations to Coulomb interactions. The approach demonstrates that particles effectively interact via the potential of mean force and that the range of this force determines the size of the collision volume. When applied to a plasma, the collision operator is shown to be related to the effective potential theory [Baalrud and Daligault, Phys. Rev. Lett 110, 235001 (2013)]. The relationship between this and previous kinetic theories is discussed.

Equations

$$\Gamma =\frac{(Ze{)}^2/4\pi {ϵ}_{0}a}{k_{\text{B}}T}.$$

$$\displaystyle\biggl{[}\frac{\partial}{\partial t}$$

$$=-\sum _{i=1}^n\int d{\mathbf{\Gamma }}_{n+1}{\mathcal{L}}_{i,n+1}^C{f}^{(n+1)}({\mathbf{r}}^{n+1},{\mathbf{v}}^{n+1},t),$$

$$f^{(n)}({\mathbf{r}}^n,{\mathbf{v}}^n,t)=\frac{N!}{(N-n)!}\int d{\mathbf{\Gamma }}^{(N-n)}{f}^{[N]}({\mathbf{r}}^N,{\mathbf{v}}^N,t)$$

$$\biggl{[}\frac{\partial}{\partial t}+\sum_{i=1}^{N}\biggl{(}\mathcal{L}_{i}+\mathcal{L}^{\textrm{ext}}_{i}+\sum_{\stackrel{{\scriptstyle j=1}}{{j\neq i}}}^{N}\mathcal{L}_{ij}^{C}\biggr{)}\biggr{]}f^{[N]}(\mathbf{r}^{N},\mathbf{v}^{N},t)=0,$$

$${\mathcal{L}}_i={\mathbf{v}}_i\cdot \frac{\partial }{\partial {\mathbf{r}}_i}+\frac{q_i}{m_i}({\mathbf{v}}_i\times \mathbf{B})\cdot \frac{\partial }{\partial {\mathbf{v}}_i}$$

$${\mathcal{L}}_i^{\text{ext}}=\frac{1}{m_i}{\mathbf{F}}_i^{\text{ext}}\cdot \frac{\partial }{\partial {\mathbf{v}}_i}$$

$${\mathcal{L}}_{ij}^C\equiv \frac{1}{m_i}{\mathbf{F}}_{ij}^C\cdot \frac{\partial }{\partial {\mathbf{v}}_i}$$

$${\mathbf{F}}_{ij}^C=\frac{q_i{q}_j}{4\pi {ϵ}_o}\frac{({\mathbf{r}}_i-{\mathbf{r}}_j)}{|{\mathbf{r}}_i-{\mathbf{r}}_j{|}^3}=-{\nabla }_i\varphi (r_{ij}).$$

$$\biggl{(}\frac{\partial}{\partial t}+\mathcal{L}_{1}+\mathcal{L}^{\textrm{ext}}_{1}\biggr{)}f^{(1)}(1)=-\int d\mathbf{\Gamma}_{2}\,\mathcal{L}_{12}^{C}f^{(2)}(1,2),$$

$$\displaystyle\biggl{(}\frac{\partial}{\partial t}$$

$$\displaystyle=-\int d\mathbf{\Gamma}_{3}\biggl{(}\mathcal{L}_{13}^{C}+\mathcal{L}_{23}^{C}\biggr{)}f^{(3)}(1,2,3),$$

$$\Delta f^{(n+1)}\equiv f_{o}^{(n+1)}\biggl{(}\frac{f^{(n+1)}}{f_{o}^{(n+1)}}-\frac{f^{(n)}}{f_{o}^{(n)}}\biggr{)}.$$

$$f_o^{(n)}({\mathbf{r}}^n,{\mathbf{v}}^n)={\rho }^{(n)}({\mathbf{r}}^n)f_{\text{M}}^{(n)}({\mathbf{v}}^n)$$

$$f_{\textrm{M}}^{(n)}(\mathbf{v}^{n})=\biggl{(}\frac{m}{2\pi k_{\textrm{B}}T}\biggr{)}^{3n/2}\exp\biggl{(}-\sum_{i=1}^{n}\frac{m\mathbf{v}_{i}^{2}}{2k_{\textrm{B}}T}\biggr{)}$$

$${\rho }^{(n)}({\mathbf{r}}^n)=\frac{N!}{(N-n)!}\frac{1}{{\mathcal{Z}}_N}\int d{\mathbf{r}}^{(N-n)}{e}^{-(V_{\text{ext}}+V_N)/k_{\text{B}}T}$$

$$\biggl{[}\frac{\partial}{\partial t}+\sum_{i=1}^{n}(\mathcal{L}_{i}+\bar{\mathcal{L}}_{i}^{(n)})\biggr{]}f^{(n)}=-\sum_{i=1}^{n}\int d\mathbf{\Gamma}_{n+1}\mathcal{L}_{i,n+1}^{C}\Delta f^{(n+1)},$$

$${\bar{\mathcal{L}}}_i^{(n)}$$

$$={\mathcal{L}}_i^{\text{ext}}+\sum _{\stackrel{j=1}{j\ne i}}^n{\mathcal{L}}_{ij}^C+\int d{\mathbf{\Gamma }}_{n+1}{\mathcal{L}}_{i,n+1}^C\frac{f_o^{(n+1)}}{f_o^{(n)}}$$

$${\bar{\mathbf{F}}}_i^{(n)}({\mathbf{r}}^n)$$

$${\bar{\mathbf{F}}}_i^{(n)}({\mathbf{r}}^n)=-{\nabla }_i{w}^{(n)}({\mathbf{r}}^n),$$

$$w^{(n)}({\mathbf{r}}^n)=-k_{\text{B}}T\ln \left[\frac{{\rho }^{(n)}({\mathbf{r}}^n)}{{\rho }^n}\right],$$

$${\nabla }_1{\rho }^{(n)}({\mathbf{r}}^n)$$

$$=\frac{\rho }{k_{\text{B}}T}\int d{\mathbf{r}}_{n+1}{\mathbf{F}}_{1,n+1}{\rho }^{(n+1)}({\mathbf{r}}^{n+1}).$$

$$\biggl{(}\frac{\partial}{\partial t}+\mathcal{L}_{1}+\bar{\mathcal{L}}_{1}^{(1)}\biggr{)}f^{(1)}(1)=-\int d\mathbf{\Gamma}_{2}\mathcal{L}_{12}^{C}\Delta f^{(2)}(1,2).$$

$$f_{\sigma }^{(n)}({\mathbf{r}}^n,{\mathbf{v}}^n,t)\equiv \frac{N!}{(N-n)!}{\int }_{\sim V_{\sigma }}d{\mathbf{\Gamma }}^{(N-n)}{f}^{[N]}({\mathbf{r}}^N,{\mathbf{v}}^N,t).$$

$$\displaystyle\biggl{(}\frac{\partial}{\partial t}+\mathcal{L}_{1}+\bar{\mathcal{L}}_{1}^{(1)}$$

$$-{\int }_{\sim V_{\sigma }}d{\mathbf{r}}_2\int d{\mathbf{v}}_2{\mathcal{L}}_{12}^C\Delta f^{(2)}(1,2)$$

$$w_{\text{DH}}^{(2)}(r)=\frac{\varphi (r)}{k_{B}T}{e}^{-r/{\lambda }_D}$$

$$\biggl{(}\frac{\partial}{\partial t}+\mathcal{L}_{1}+\bar{\mathcal{L}}_{1}^{(1)}\biggr{)}f^{(1)}(1)=-\oint_{S_{2}}d\mathbf{v}_{2}d\mathbf{S}_{2}\cdot\mathbf{u}f^{(2)}(1,2).$$

$$\biggl{(}\frac{\partial}{\partial t}+\mathcal{L}_{1}+\mathcal{L}_{2}+\bar{\mathcal{L}}_{1}^{(2)}+\bar{\mathcal{L}}_{2}^{(2)}\biggr{)}f^{(2)}(1,2)=0.$$

$$\displaystyle\biggl{(}\frac{\partial}{\partial t}+\mathcal{L}_{1}+\bar{\mathcal{L}}_{1}^{(1)}\biggr{)}f^{(1)}(1)$$

$$=\int d\omega d{\mathbf{v}}_{2}u[f^{(2)}({\mathbf{r}}_1,{\mathbf{v}}_1,{\mathbf{r}}_2^{+},{\mathbf{v}}_2)-f^{(2)}({\mathbf{r}}_1,{\mathbf{v}}_1,{\mathbf{r}}_2^{-},{\mathbf{v}}_2)].$$

Full text

I Introduction

Derivations of standard plasma kinetic equations confront either infrared divergences (Boltzmann equation) grad:58 ; harr:71 ; ferz:72 , ultraviolet divergences (Lenard-Balescu equation) lena:60 ; bale:60 ; guer:60 , or both (Landau equation) land:36 . These divergences are resolved by invoking physical arguments, such as Debye shielding to resolve the infrared divergence, or the distance of closest approach in a binary collision to resolve the ultraviolet divergence. The results obtained in either approach agree to logarithmic accuracy as long as the plasma is weakly coupled and near local equilibrium. They have been validated for many collisional transport processes satisfying these limits over the several decades since they were developed. Despite the success of these equations, the theory stands in an unsatisfactory position in that it relies on ad hoc arguments. As a result, it is difficult to attempt generalizations to other important situations, such as moderate or strong Coulomb coupling, or strong magnetization.

Using a new self-consistent closure of the BBGKY hierarchy, we present a plasma kinetic equation that (i) is convergent, i.e. it does not confront either infrared or ultraviolet divergences, and that (ii) provides a rigorous extension of plasma kinetic theory into the regime of strong Coulomb coupling. The approach demonstrates that particles effectively interact via the potential of mean force and that the range of this force determines the size of the collision volume.

In order to put our approach into perspective, we first briefly recall the origin of standard kinetic equations. Each can be derived from a perturbative closure of the BBGKY hierarchy obtained by first identifying a small dimensionless parameter characteristic of the system, and by then using its smallness to truncate the hierarchy at the level of two-particle correlations. Two dimensionless parameters naturally arise in terms of the average particle density $n$, and the range $l$ and the strength $\phi_{0}$ of the interaction potential $\phi(r)$. Namely, the “concentration parameter” $nl^{3}$, which measures how many particles, on the average, simultaneously interact with a given particle, and the “strength parameter” $\phi_{0}/k_{\textrm{B}}T$, which measures the interaction potential energy of two colliding particles in units of the mean particle kinetic energy. It is easily seen that, strictly speaking, these parameters are not well adapted to plasmas since the Coulomb potential has an infinite range, and its magnitude becomes arbitrarily strong at short distance. The celebrated Boltzmann equation is obtained by assuming $nl^{3}\ll 1$ and $\phi_{0}/k_{\textrm{B}}T\sim O(1)$; it applies to dilute gases interacting via short range potentials of arbitrary interaction strength. grad:58 ; harr:71 ; ferz:72 The Landau equation is obtained in the so-called weakly coupled limit characterized by $\phi_{0}/k_{\textrm{B}}T\ll 1$ and $nl^{3}\sim O(1)$; it applies to weak potentials, i.e. potentials that are uniformly small for all interparticle distances $r$, a condition that is generally not respected due to the typical strong repulsions at short distances. land:36 The Lenard-Balescu equation, which was developed to deal with weakly coupled plasmas, is obtained in the so-called weakly-coupled, long-range limit assuming $\phi_{0}/k_{\textrm{B}}T\ll 1$ and $nl^{3}\phi_{0}/k_{\textrm{B}}T\sim O(1)$. lena:60 ; bale:60 ; guer:60 Plasmas do not fall in any of these categories, which explains why one needs to regularize ad-hoc the standard kinetic equations before applying them to plasmas. A solution to the problem can be obtained by writing the Coulomb potential as a sum of two terms: a weak long-range term, plus a strong short-range term. This amounts to carefully joining the Lenard-Balescu and Boltzmann equations (see, e.g., Refs. thom:60, ; hubb:61, ; bald:62, ; frie:63, ; kiha:63, ; wein:64, ; aono:65, ; goul:67, ).

Until today, no one has derived a practical kinetic equation beyond the previous limits. The generalization of the Boltzmann equation to higher densities remains an unsolved challenge (much progress on the kinetic theory of dense gases and liquids has been made, but for systems in thermal equilibrium). It was shown that the systematic inclusion of many-body collisions through a density expansion similar to the virial expansion for computing the equilibrium properties of dense gases is plagued by unphysical divergences.dorf:67 This is because $N$-body collisions and correlations cannot be treated separately from those of all higher orders ($N\!+\!1$-body, $N\!+\!2$-body, etc.). This difficulty is symptomatic of many perturbative approaches in physics and could in principle be dealt with using so-called renormalization techniques. While such renormalization schemes have been proposed,dali:11 their exceeding complexity has thus far restricted their usefulness to formal, but not practical, solutions.

In this paper, we present a closure scheme for the BBGKY hierarchy that, like the standard equations, results from the identification of a small expansion parameter but, unlike the concentration and strength parameters, does not explicitly depend on the characteristics of the interaction potential like $l$ and $\phi_{0}$. The present expansion parameter [Eq. (11) below] is a measure of the perturbation of the distribution function about thermal equilibrium. The closure retains the correct equilibrium limit at all orders, and ensures that screening is captured in the near-equilibrium limit, while at the same time allowing inclusion of short-range interactions. It results in a convergent and tractable kinetic equation that provides a rigorous extension of plasma kinetic theory into the regime of strong Coulomb coupling. Here, strong Coulomb coupling refers to plasmas in which the average Coulomb potential energy of interacting particles $(Ze)^{2}/4\pi\epsilon_{0}a$ ($Ze$ is the electric charge, $a=(3/4\pi n)^{1/3}$ is the average interparticle spacing) exceeds their average kinetic energy $k_{\textrm{B}}T$ ($T$ is the temperature), i.e. $\Gamma\gtrsim 1$, where

[TABLE]

The resulting collision operator is the same as what is obtained from the Boltzmann equation if the interparticle force in binary collisions is taken to be the potential of mean force, hill:60 rather than the Coulomb potential. In the weakly coupled limit, the potential of mean force is the Debye-Hückel potential deby:23 and the theory gives the same results as standard plasma kinetic equations. baal:14 The transport coefficients resulting from the convergent kinetic equation are those of the recent “effective potential theory” (EPT), which was shown to successfully extend the conventional plasma transport theory into the strongly coupled regime.baal:13 ; baal:14 ; baal:15 Previous work has shown that EPT can extend plasma kinetic theory well into the regime of strong Coulomb coupling. dali:14 ; bezn:14 ; haxh:14 ; stri:16 ; dali:16 ; shaf:17

In addition, in the resulting hydrodynamic equations the pressure and internal energy are composed of two parts; a kinetic (or ideal) part and a potential (or excess) part. The ideal part, which is obtained with all kinetic equations, is the standard ideal gas component that represents the transfer of momentum or energy due to the flow of particles. The excess part, which is typically absent from standard kinetic theories, represents the transfer of momentum or energy between particles by the particle interactions.

The paper is organized as follows. In Section II, we introduce the expansion parameter at the basis of the closure of the BBGKY hierarchy. In Section III, the convergent kinetic equation resulting from this closure is derived and its properties are discussed. To this end, we closely follow the method used by Grad in his derivation of the Boltzmann equation. grad:58 In Section IV, the convergent kinetic equation is discussed in comparison to the standard plasma kinetic equations.

II Basic expansion parameter

Kinetic theories can be derived from the BBGKY hierarchy ferz:72

[TABLE]

where

[TABLE]

defines the $n^{\textrm{th}}$-order reduced distribution functions in terms of the $N$-particle distribution function $f^{[N]}$. Here, $\mathbf{r}^{N}=(\mathbf{r}_{1},\mathbf{r}_{2},\ldots\mathbf{r}_{N})$, where $\mathbf{r}_{i}$ denotes the spatial location of particle $i$, $\mathbf{v}^{N}=(\mathbf{v}_{1},\mathbf{v}_{2},\ldots\mathbf{v}_{N})$ where $\mathbf{v}_{i}$ denotes the velocity of particle $i$, $d\mathbf{\Gamma}_{n}\equiv d\mathbf{r}_{n}d\mathbf{v}_{n}$ is a shorthand notation for the 6-dimensional phase-space, and $d\mathbf{\Gamma}^{(N-n)}\equiv d\mathbf{\Gamma}_{n+1}\ldots d\mathbf{\Gamma}_{N}$. The BBGKY hierarchy follows from integrating the Liouville equation

[TABLE]

in order to obtain an evolution equation for each reduced distribution function. In equation (2),

[TABLE]

where the second term is associated with the Lorentz force due to an external magnetic field, and

[TABLE]

where $\mathbf{F}_{i}^{\textrm{ext}}=-\nabla_{i}\phi_{\textrm{ext}}(\mathbf{r}_{i})$ is an external force, associated with an external potential. The external potential is added here for later convenience to introduce spatial inhomogeneities in the thermal equilibrium state. The term

[TABLE]

is associated with the electrostatic Coulomb interactions between particles

[TABLE]

Any kinetic theory aims to describe the evolution of the first-order reduced distribution function $f^{(1)}(\mathbf{r},\mathbf{v},t)$ ($n=1$):

[TABLE]

which depends on $f^{(2)}(1,2)$ [$(1,2)$ is shorthand notation for $(\mathbf{r}_{1},\mathbf{v}_{1},\mathbf{r}_{2},\mathbf{v}_{2},t)$]. The task of a collision operator is to provide an approximate expression for $f^{(2)}(1,2,t)$. For example, this can be related to $f^{(3)}(1,2,3,t)$ via the $n=2$ equation

[TABLE]

but approximations must be introduced to close the hierarchy.

The most famous closure is to neglect the third order distribution function $f^{(3)}\rightarrow 0$ in Eq. (10), which is an excellent approximation for dilute gasses that interact via a short-range force (in the notation of the introduction, $nl^{3}\ll 1$). One then solves the homogeneous equation for $f^{(2)}$ subject to a constraint on the length scale over which binary collisions occur, as well as a lack of initial correlations in a binary scattering event. This method was shown to lead to the Boltzmann equation. The calculation is not trivial and several variations exist; below we shall rely on the method proposed by Grad. grad:58 However, as recalled in the introduction, this closure is a poor approximation for plasmas. The reason is that the physics of screening is contained in $f^{(3)}$. Neglecting screening not only misses an important physical process, but the resulting kinetic equation diverges because the two-body Coulomb force has an infinite range.

The closure scheme proposed here enforces that the correct equilibrium (i.e., thermodynamic) limit is maintained at all orders. This ensures that screening is captured, while at the same time allowing inclusion of short-range interactions. This can be accomplished by taking the basic expansion parameter to be

[TABLE]

Equation (11) measures the difference of the $n+1$ and $n$ probability distributions, referenced to their equilibrium values. It is a measure of the perturbation of non-equilibrium correlations about equilibrium; $\Delta f^{(n+1)}\rightarrow 0$ represents that correlations approach their value at equilibrium. In Eq. (11),

[TABLE]

is the equilibrium reduced distribution function,

[TABLE]

is the Maxwellian velocity distribution function,

[TABLE]

is the $n-$particle density distribution function, $V_{N}(\mathbf{r}^{N})=\sum_{i=1}^{N}\sum_{j>i}^{N}\phi(r_{ij})$ is the electrostatic potential energy, $V_{\textrm{ext}}(\mathbf{r}^{N})=\sum_{i=1}^{N}{\phi_{\textrm{ext}}(\mathbf{r}_{i})}$ is the interaction energy with the external potential $\phi_{\textrm{ext}}$, and $\mathcal{Z}_{N}=\int\exp\left[-(V_{\textrm{ext}}+V_{N})/k_{\textrm{B}}T\right]d\mathbf{r}^{N}$ is the configurational integral.

The BBGKY hierarchy from Eq. (2) can be rearranged in a way that expresses $\Delta f^{(n+1)}$ as an expansion parameter

[TABLE]

where

[TABLE]

is an operator associated with force

[TABLE]

The force (II) has a simple physical interpretation. It is the mean force acting on particle $i$ obtained when keeping both it and a set of other particles ($j=1,\ldots,n$, excluding $j=i$) at fixed positions and averaging over all equilibrium configurations of the other $N-n$ particles. As shown in Appendix A, this statistical force can be expressed as the gradient of a potential,

[TABLE]

where the potential of mean force is

[TABLE]

where $\rho=\rm lim_{V\to\infty}{\int_{V}{d\mathbf{r}\rho^{(1)}({\bf r})}/V}$ is the average particle density. In the absence of the external potential, $\phi_{\textrm{ext}}\equiv 0$, $w^{(n)}(\mathbf{r}^{n})=-k_{\textrm{B}}T\ln g^{(n)}(\mathbf{r}^{n})$, where $g^{(n)}(\mathbf{r}^{n})=\rho^{(n)}(\mathbf{r}^{n})/\rho^{n}$ is the $n$-particle distribution function (when $\phi_{\textrm{ext}}\neq 0$, $g^{(n)}(\mathbf{r}^{n})=\rho^{(n)}(\mathbf{r}^{n})/\Pi_{i=1}^{n}\rho^{(1)}(\mathbf{r}_{i})$).

Although Eq. (15) is equivalent to Eq. (2), writing it in this way affords certain pedagogical clarities. The right-hand side, at any order $n$, can now be interpreted as a “collision operator” in the sense that it vanishes at equilibrium and is small compared to the left side of the equation for slight perturbations from equilibrium. This contrasts with the right side of Eq. (2), which does not vanish at equilibrium. Of course, solving Eq. (15) still requires a closure. The scheme suggested here is that the dynamical evolution of $f^{(n)}$ be closed at order $n$ by taking $\Delta f^{(n+1)}\rightarrow 0$. However, a closure for the equilibrium distribution $\rho^{(n)}(\mathbf{r}^{n})$ (or $g^{(n)}(\mathbf{r}^{n})$ in the homogeneous case) is still required in order to determine the potential of mean force arising on the left side of the equation. Determining $\rho^{(n)}$ is a more tractable problem because one can rely on methods of equilibrium statistical mechanics. At equilibrium, the BBGKY hierarchy reduces to the Yvon-Born-Green (YBG) limit hans:06

[TABLE]

Although this is also a hierarchical equation, accurate approximations have been developed for most potentials of interest and for quite general conditions of density and temperature. hans:06 We also note that Eq. (19) follows directly from Eq. (20). To proceed with the kinetic theory derivation, $\rho^{(n)}$ will be considered a known quantity.

In particular, in the following section a kinetic equation is derived from the second order ($n=2$) term of Eq. (15) assuming that $\Delta f^{(3)}=0$. The derivation closely follows Grad’s method. grad:58 Many methods have been used to derive the Boltzmann equation from Eq. (2), ferz:72 ; cohe:61 and any of these could also be used for our purposes. We choose Grad’s method because it follows directly from the closure $f^{(3)}\rightarrow 0$ (or $\Delta f^{(3)}\rightarrow 0$ in our modification) without the need to introduce additional complications associated with a cluster expansion. It will also give a clear description of the effective collision volume in a Coulomb system. The closure for $\rho^{(2)}(\mathbf{r}_{1},\mathbf{r}_{2})$ will be assumed known. We will focus on the near-homogenous limit that is relevant to fluid theory. Here, the two-particle density consists of a component that varies on large scales indicative of the fluid-scale gradients, and another short spatial scale indicative of the scale of interactions: $\rho^{(2)}\approx\rho(\mathbf{r}_{1})g^{(2)}(\mathbf{r}_{1},r)$, where $r=|\mathbf{r}_{1}-\mathbf{r}_{2}|$. The radial distribution function $g^{(2)}(r)$ can be provided by the hypernetted chain approximation, hans:06 which is known to be accurate for plasmas at conditions spanning weak to strong coupling. However, the theory does not depend on the method used to obtain $g^{(2)}$. At weak coupling, $g^{(2)}$ simply asymptotes to the Debye-Hückel limit and conventional plasma kinetic theory will result from this unified framework.

III Plasma Kinetic Equation

III.1 Modified Grad Method

The plasma kinetic equation (9) in terms of the expansion parameter of Eq. (11) is

[TABLE]

Following Grad’s method, grad:58 space is divided into a volume within which collisions occur $V_{\sigma}$, and a volume outside of this in which the observable (truncated) distribution functions are defined

[TABLE]

Here, $\sim V_{\sigma}$ denotes that the spatial integral excludes the small volume $V_{\sigma}$. By integrating Liouville’s equation (4), the first-order BBGKY hierarchy equation for the “truncated reduced distribution function” is

[TABLE]

where $\mathbf{u}\equiv\mathbf{v}_{1}-\mathbf{v}_{2}$. Here, an extra term arises in comparison to Eq. (21) that depends on the two-body distribution function evaluated at the surface of the collision volume, $\mathbf{S}_{2}$. The method relies on the assumption that the collision volume, $V_{\sigma}$, is sufficiently small that there is an appropriate limit in which $f_{\sigma}^{(n)}\rightarrow f^{(n)}$ (i.e., that $f_{\sigma}^{(n)}$ is the observable of interest).

A key question is, what determines the scale of $V_{\sigma}$? For a dilute gas, this is typically associated with the range of the short-range intermolecular forces. In these theories, $f^{(2)}$ replaces $\Delta f^{(2)}$ in Eq. (III.1) and the last term is considered negligible because the intermolecular force (inside the $\mathcal{L}_{12}$ operator) is small in the region of space outside 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 to determine $f^{(2)}$ on its surface. This leads to the Boltzmann equation. grad:58

Clearly, the same argument does not apply to plasmas because the Coulomb force is long-range and effectively extends over all of space. Indeed, this is why application of the Boltzmann equation to plasma results in a divergence. However, writing Eq. (III.1) in terms of $\Delta f^{(2)}$ introduces the notion of screening because the effective range of interaction is associated with the spatial correlation scale of $\Delta f^{(2)}$ itself. Of course $\Delta f^{(2)}$ vanishes at equilibrium at all spatial scales, so the relevant correlation arises only away from equilibrium. Considering a small perturbation from equilibrium, $\Delta f^{(2)}$ will always become small on large scales at which two particles become decorrelated. Only at sufficiently small spatial scales (i.e., the collision scale $|\mathbf{r}_{1}-\mathbf{r}_{2}|<V_{\sigma}$) will $\Delta f^{(2)}$ be appreciable in magnitude. For instance, consider a homogenous plasma. Starting at large scales where spatial correlations are small $f^{(2)}(1,2)\approx\rho^{2}f^{(1)}(\mathbf{v}_{1})f^{(2)}(\mathbf{v}_{2})$, the spatial component of $\Delta f^{(2)}$ will have the possibility of contributing at the scale associated with spatial correlation of $f_{o}^{(2)}$; i.e., $\Delta f^{(2)}\sim\rho^{2}-\rho^{(2)}\approx\rho^{2}[1-g(r)]$. This scale is directly determined by the potential of mean force via Eq. (19), and thus the screening length. For example, in a weakly coupled plasma the potential of mean force is the Debye-Hückel potential deby:23

[TABLE]

where $\lambda_{D}$ is the Debye length. This corresponds directly to the conventional notion that $V_{\sigma}$ is a “Debye interaction sphere” at weak coupling. However, the potential of mean force is a more general concept that extends to strongly coupled plasmas as well.

Since $V_{\sigma}$ is defined as the excluded region of the volume integral, the last term in Eq. (III.1) is negligible and the kinetic equation is

[TABLE]

Solving for the collision operator in Eq. (25) requires solving for $f^{(2)}$ on the surface of the collision volume. As discussed in the previous section, this will be determined from the second order equation [$n=2$ in Eq. (15)] using $\Delta f^{(3)}=0$ as a closure

[TABLE]

As does the derivation of the Boltzmann equation, Eq. (26) describes the dynamical evolution of two interacting particles. However, retention of the higher order terms arising in the operators $\bar{\mathcal{L}}_{i}^{(n)}$ shows that the particles effectively interact via the mean force, rather than the bare Coulomb force $\mathbf{F}_{ij}^{C}$. In this way, surrounding particles “mediate” the binary interactions; a concept Rostoker referred to as “dielectric dressing”. rost:64 Thus, although the dynamics are two-body, the interaction is more general than that used in the Boltzmann equation. In this sense, it relaxes the binary collision assumption to a certain degree by including the statistical influence of all other $N-2$ particles.

III.2 Collision Operator

The vast majority of plasmas are magnetized in the sense that the Lorentz force may influence macroscopic dynamics, but it does not influence microscopic dynamics within the collision volume. baal:17 In a weakly coupled plasma, this condition is satisfied if the gyroradius of most particles is signficantly larger than the Debye length: $r_{c}\gg\lambda_{D}$, where $r_{c}=v_{T}/\omega_{c}$ is the thermal gyroradius and $\omega_{c}=|Ze|B/m$ is the gyrofrequency. In a strongly coupled plasma, it requires $r_{c}\gg a$. baal:17 Most of plasma kinetic theory is based on this scale separation, and this section concentrates on this situation. In this case, the Lorentz force term in $\mathcal{L}_{1}$ contributes to the plasma kinetic equation (21), but it is negligible in the operators $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ in Eq. (26), which describes $f^{(2)}$ inside the collision volume.

Following the method of Ref. grad:58, , we choose a coordinate system such that a diametral plane intersecting the sphere is aligned perpendicular to the relative velocity vector $\mathbf{u}$. Using polar coordinates in this plane, an area element is denoted $d\omega=rdrd\epsilon$. This transformation maps the projection of $\mathbf{S}_{2}$ onto the disk $0\leq r\leq\sigma$ and $0\leq\epsilon\leq 2\pi$. The disk $\omega$ is covered twice: once by the hemisphere $S_{2}^{+}:\,\mathbf{u}\cdot d\mathbf{S}_{2}>0$ and once by $S_{2}^{-}:\,\mathbf{u}\cdot d\mathbf{S}_{2}<0$. Points that map onto $S_{2}^{+}$ represent particles moving away from one another (post-collision particles), while those that map onto $S_{2}^{-}$ represent particles moving toward one another (pre-collision particles). Denoting the points that map onto $S_{2}^{+}$ as $\mathbf{r}_{2}^{+}$ and those that map onto $S_{2}^{-}$ as $\mathbf{r}_{2}^{-}$, we find $\mathbf{u}\cdot d\mathbf{S}_{2}=-ud\omega$ on $S_{2}^{+}$ and $\mathbf{u}\cdot d\mathbf{S}_{2}=ud\omega$ on $S_{2}^{-}$. In terms of this transformation, Eq. (25) becomes

[TABLE]

To obtain an explicit collision operator, we invoke four assumptions that are similar to those underlying the Boltzmann equation: (1) truncation (finite size of $V_{\sigma}$), (2) binary collisions, (3) modified molecular chaos, and (4) slow variation of $f^{(1)}$. Assumption (1) was already applied above, which justified the neglect of the last term in Eq. (III.1), and that $f_{\sigma}^{(1)}$ and $f_{\sigma}^{(2)}$ accurately approximate $f^{(1)}$ and $f^{(2)}$. We point out again that assumption (1) has a different origin in the context of the theory presented here, in comparison to the Boltzmann equation. Here, the small size of $V_{\sigma}$ is justified only by the limited range of the effective interaction implied by $\Delta f^{(2)}$. In contrast, in terms of $f^{(2)}$ the neglected term in Eq. (III.1) would diverge for a plasma. Neglecting this is justified in derivations of the Boltzmann equation because they focus on dilute gases which have short-range potentials.

Since the collision volume is small, the number of particle inside of it is very small compared to the total number of particles. Furthermore, we assume that within the collision volume particle interactions are predominately binary [assumption (2)], but that they occur via the mean force. This justifies $\Delta f^{(3)}=0$ as a closure, leading to Eq. (26) within $V_{\sigma}$. The solution of this equation is $f^{(2)}=\textrm{const}$ on a 2-particle trajectory. Thus, the points that map into the hemisphere $S_{2}^{+}$ can be equated with post-collision coordinates: $f^{(2)}(\mathbf{r}_{1},\mathbf{v}_{1},\mathbf{r}_{2}^{+},\mathbf{v}_{2},t)\rightarrow f^{(2)}(\hat{\mathbf{r}}_{1},\hat{\mathbf{v}}_{1},\hat{\mathbf{r}}_{2},\hat{\mathbf{v}}_{2},\hat{t})$.

To solve for the particle trajectories inside the sphere, the Boltzmann equation typically makes the “molecular chaos” approximation (3) whereby the two-particle distribution is assumed to be uncorrelated in both the initial and final states: $f^{(2)}(1,2)=f^{(1)}(1)f^{(1)}(2)$. This is a rather more subtle point associated with the introduction of irreversibility into the theory, which is discussed in depth in Ref. grad:58, . The procedure presented here attempts to develop a model that is consistent with the equilibrium state. As such, we modify this approximation slightly to account for the statistical spatial correlation of particles at the surface of the collision volume: $f^{(2)}(1,2)|_{r=\sigma}\approx\chi f^{(1)}(1)f^{(1)}(2)$. This concept, first introduced by Enskog ensk:22 in the context of a hard-sphere gas, accounts for the “excluded volume” associated with the fact that hard spheres cannot overlap. For hard spheres, $\chi=g(r=\sigma)$. For a soft potential, $\sigma$ is not a single value for all interactions, and as such the concept is understood as a statistical analog. hanl:72 This is discussed in Ref. baal:15, for single-component plasmas. In this model, $\sigma$ is associated with the location where $g(r=\sigma)=0.87$ and $\chi$ is determined by using this effective diameter in the virial expansion of the Enskog equation of state; see also Ref. shaf:17, for an extension to binary mixtures. The value $0.87$ comes from a fit to molecular dynamics simulations of the OCP, but is expected to be a universal value for repulsive potentials. baal:15 ; dali:16 From this model $\chi\rightarrow 1$ for weakly coupled plasmas $\Gamma\lesssim 1$, and $\chi\approx 1.3-1.4$ for $1\lesssim\Gamma\lesssim 30$.

Finally, we invoke approximation (4), which states that the collision time is much shorter than the timescale associated with the evolution of $f^{(1)}$. In this limit, we expect the spatial coordinates $\hat{\mathbf{r}}_{1}$, $\mathbf{r}_{2}$ and $\hat{\mathbf{r}}_{2}$ to differ from $\mathbf{r}_{1}$ by a small amount on the order of $\sigma$. Likewise, the short collision time implies that $\hat{t}$ differs from $t$ by a negligible amount. Putting the results of these four approximations into Eq. (27) leads to the plasma kinetic equation

[TABLE]

where it is understood that each occurrence of $f^{(1)}$ is evaluated at the same local position $\mathbf{r}_{1}$ and time $t$. Equation (28) was derived for a one-component plasma. The generalization to a multicomponent system

[TABLE]

follows directly by tracking the derivation above and distinguishing species $s$ from $s^{\prime}$.

Since the mean force is central and conservative, the differential area of the disk can be related to a scattering cross section, $d\omega=\sigma_{ss^{\prime}}\,d\Omega$ where $d\Omega=d\phi d\theta\sin\theta$, or an impact parameter, $d\omega=bdbd\phi$, in the usual way. grad:58 The results of the standard two-body scattering problem are determined by the scattering angle $\theta=\pi-2\Theta$, where

[TABLE]

Here, $m_{ss^{\prime}}=m_{s}m_{s^{\prime}}/(m_{s}+m_{s^{\prime}})$ is the reduced mass, and $r_{o}$ is the distance of closest approach, which is determined by the largest root of the term in square brackets. Note that this is the same as the scattering angle used in the Boltzmann equation, except that the potential of mean force $w^{(2)}$ replaces the bare interaction potential $\phi(r)$.

III.3 Comment on Strongly Magnetized Plasmas

If the external magnetic field is sufficiently strong that the gyromotion of particles fits within the collision volume, the situation is much more complicated. O’Neil considered very strongly magnetized one-component plasmas (such that $r_{c}\ll r_{o}$), and derived a kinetic equation by aligning the surface $d\mathbf{S}_{2}$ along the magnetic field (rather than $\mathbf{u}$). onei:83 A collision operator was then derived based on an assumption that the adiabatic invariant $|u_{+}|^{2}/B$ is preserved during the interaction. Here $u_{+}\equiv u_{x}+iu_{y}$ are Cartesian components of $\mathbf{u}$ and the magnetic field is in the $\hat{z}$ direction. This adiabatic assumption is restrictive, and a more general theory is desirable.

The method outlined above may provide a fruitful starting point for developing such a theory. The main advantage is that it self-consistently accounts for screening. A matter of continuing investigation in strongly magnetized plasmas is the relative importance of short and long range collisions. dubi:14 Early theories modified Landau-Spitzer based approaches by changing the Debye length to the gyroradius in the Coulomb logarithm. sili:63 ; mont:74 Others modified the Lenard-Balescu equation to account for magnetization in the plasma dielectric function. rost:60b ; hass:77 However, recent molecular dynamics results have shown trends inconsistent with either prediction at strong magnetization. baal:17 By accounting for screening via the potential of mean force, and for gyromotion via the 2-body interaction inside the collision sphere, one may address in a self-contained way the combined influence of short and long-range interactions. The Bohr-van Leeuwen theorem ensures that the magnetic field does not influence any statistical property at equilibrium, so the potential of mean force is expected to remain unchanged from that presented above.

III.4 Equation of State

In addition to providing a convergent collision operator, Eq. (25) introduces a term in the convective derivative ($\bar{\mathcal{L}}_{1}^{(1)}$) that is absent from the Boltzmann equation. This term is associated with non-ideal components of the equation of state, i.e. the excess (or potential) component of the pressure and internal energy. Although this is negligible in a weakly coupled plasma, it is the dominant component in moderately and strongly coupled plasmas. hans:06 The fact that these terms arise naturally should be expected because the closure is designed to ensure that the exact equilibrium state consistent with the YBG hierarchy, Eq. (20), is maintained.

To make the connection with pressure, it is first useful to notice that the mean force acting on one particle can be written as

[TABLE]

where

[TABLE]

is a second-rank tensor, $\mathbf{r}\equiv\mathbf{r}_{2}-\mathbf{r}_{1}$ and $\phi^{\prime}(r)\equiv d\phi/dr$. A derivation of Eq. (31) from (II) is provided in Appendix B. In the common limit that the plasma is sufficiently homogeneous that $\rho^{(1)}(\mathbf{r}_{1})$ is constant on the spatial scale of the two-body correlations, then $\rho^{(2)}(\mathbf{r}_{1},\mathbf{r}_{2})\approx[\rho^{(1)}(\mathbf{r}_{1})]^{2}g(\mathbf{r}_{1},r)$; see Pohl, et al pohl:04 for a related discussion. In this case, the divergence of the pressure tensor in Eq. (31) reduces to the gradient of a scalar, $\nabla_{1}p_{\Phi}$, where

[TABLE]

is the well-known equilibrium expression for the potential component of the scalar pressure. hans:06

We note that Eq. (33) makes use of the assumption of “local thermodynamic equilibrium” whereby there is a scale separation between large (fluid-scale) gradients represented by $\mathbf{r}_{1}$ and the interaction scale represented by $r=|\mathbf{r}_{1}-\mathbf{r}_{2}|$. Thus, the force represented by $\nabla_{1}p_{\Phi}$ arises due to a slight inhomogeneity of the background density across the collision volume ($L\gg\sigma$). This origin of the excess pressure can be compared with Enskog’s theory of hard sphere gases. ensk:22 ; ferz:72 Enskog’s theory expands $f^{(1)}(\mathbf{r}_{1}\pm\sigma_{\textrm{s}}\mathbf{k})$ about the local position $\mathbf{r}_{1}$, where $\sigma_{s}$ is the sphere diameter and $\mathbf{k}$ is vector connecting the center of two spheres at the instant of contact. The excess pressure (and excess internal energy) arising from the lowest-order term in this expansion comes about due to any inhomogeneity of $f^{(1)}(\mathbf{r}_{1})$ across the collision volume, which is simply the volume of the sphere in a hard-sphere gas. The higher order terms proportional to $\nabla_{1}f^{(1)}$, etc., are associated with the non-local nature of the collision itself (a hard sphere gas is a singular case in that $\sigma_{\textrm{s}}$ is both the maximum and minimum scale of an interaction). Equation (33) represents a generalization of the first of these contributions (local collisions taking place in a finite sized collision volume) to the case of arbitrary potential. Its existence in the theory is enforced by the requirement that the exact equilibrium, or local equilibrium, limit be maintained.

A closed fluid description follows from a Chapman-Enskog solution of the kinetic equation, ferz:72 but the basic features of non-ideality can be illustrated directly from the conservation equations. In particular, the density moment of Eq. (29) leads to the usual continuity equation

[TABLE]

where $\mathbf{V}\equiv\int d\mathbf{v}\,\mathbf{v}f^{(1)}/\rho$ is the hydrodynamic velocity and $\rho=\rho^{(1)}(\mathbf{r}_{1})$. This is the same as what results from the Boltzmann equation. However, the momentum moment leads to

[TABLE]

in which the total pressure tensor now includes both kinetic (ideal) and potential (excess) components

[TABLE]

Here, $d/dt=\partial/\partial t+\mathbf{V}\cdot\nabla$, $\mathcal{P}_{K}=\sum_{s}m_{s}\int d\mathbf{v}\,\mathbf{v}_{r}\mathbf{v}_{r}f_{s}^{(1)}/\rho$ is the kinetic part of the pressure tensor, $\mathbf{v}_{r}=\mathbf{v}-\mathbf{V}$, $\rho_{m}=\sum_{s}m_{s}\rho_{s}$ is the mass density, $\rho_{q}=\sum_{s}q_{s}\rho_{s}$ the charge density, $\mathbf{J}=\sum_{s}q_{s}\rho_{s}(\mathbf{V}_{s}-\mathbf{V})$ the current density, and $\mathbf{E}^{\prime}=\mathbf{E}+\mathbf{V}\times\mathbf{B}$ the electric field in the reference frame of the fluid.

The mean force will also contribute to the energy balance. Taking the $\frac{1}{2}m_{s}\mathbf{v}_{r}^{2}$ moment of the kinetic equation, and summing over species, gives

[TABLE]

where $u_{K}=\sum_{s}\frac{1}{2}m_{s}\int d\mathbf{v}\,v_{r}^{2}f_{s}^{(1)}$ is the internal kinetic energy, and $\mathbf{q}=\sum_{s}\frac{1}{2}m_{s}\int d\mathbf{v}\,\mathbf{v}_{r}v_{r}^{2}f_{s}^{(1)}$ is the kinetic component of the heat flux. Noting that $\mathbf{V}\cdot(\nabla\cdot\mathcal{P}_{\Phi})=\nabla\cdot(\mathbf{V}\cdot\mathcal{P}_{\Phi})-\mathcal{P}_{\Phi}\colon\nabla\mathbf{V}$, Eq. (37) can alternatively be written

[TABLE]

where $\mathbf{q}=\mathbf{q}_{K}+\mathbf{q}_{\Phi}$, and

[TABLE]

is a contribution to the heat flux from the mean force that is absent in the Boltzmann equation.

An evolution equation for the total internal energy must also include the contribution from the average interparticle potential energy

[TABLE]

The evolution of $u_{\Phi}$ is described by the $\frac{1}{2}\phi(r)$ moment of the second order equation of the BBGKY hierarchy

[TABLE]

Explicit evaluation of this requires a time-dependent solution for $f^{(2)}$. However, if we apply a Kirkwood-like superposition approximation kirk:46 at second order, $f^{(2)}(1,2)\approx\rho^{(2)}(\mathbf{r}_{1},\mathbf{r}_{2})f^{(1)}(\mathbf{v}_{1})f^{(2)}(\mathbf{v}_{2})$, as in Sec. III.2 above, Eq (40) reduces to $\rho u_{\Phi}=\frac{1}{2}\int\phi(r)\rho^{(2)}(\mathbf{r}_{1},\mathbf{r}_{2})d\mathbf{r}_{2}|_{\mathbf{r}_{1}=\mathbf{r}}$, and Eq. (III.4) to $du_{\Phi}/dt=0$. In this limit, the total energy evolution equation then has the form

[TABLE]

where $u=u_{K}+u_{\Phi}$ is the total internal energy. Note that because of the non-ideal terms, thermodynamic relations must be invoked to relate this to temperature. ferz:72 In the limit that the plasma is locally homogeneous, the familiar expression for the potential energy density hans:06

[TABLE]

is obtained.

IV Comparison with previous models

Common plasma theories can also be derived from closures of the BBGKY hierarchy. Here, these closures are compared with the one presented above. Specific attention is paid to the approximations that lead to divergences by considering the equilibrium limit of the second order distribution function resulting from each closure.

A common first step in the derivation of these equations is to define a cluster expansion of the form

[TABLE]

The cluster expansion allows aspects of the triplet distribution function to be retained via products of $f^{(1)}$ and $g_{2}$. Applying this, Eq. (10) is

[TABLE]

Each of the theories discussed below are based on a binary collision approximation in which the triplet correlation term Eq. (45d) is neglected, but where each keeps a different subset of terms (45a)–(45c) selected, as illustrated below for the Landau and Boltzmann equations, by identifying the leading order in the expansion parameters discussed in the introduction.Balescu_book

Divergences arise in each, not because of the neglect of triplet correlations, but rather because of neglecting one or more of the terms (45b)–(45c). Although the collision operator vanishes at equilibrium due to the velocity dependence of the Maxwellian distribution, the divergences are associated with the spatial component of the distribution function. Thus, they can be illustrated by considering just the spatial dependence of the collision operator at equilibrium

[TABLE]

where

[TABLE]

Landau Equation: As recalled in the introduction, the Landau equation is obtained in the weakly coupled limit characterized by $\lambda=\phi_{0}/k_{\textrm{B}}T\ll 1$. In terms of this perturbative parameter, the Liouville operators satisfy $\mathcal{L}_{i}=\mathcal{O}(\lambda^{0})$ and $\mathcal{L}_{ij}=\mathcal{O}(\lambda^{1})$, while the distribution functions $f=\mathcal{O}(\lambda^{0})$ (since it carries the complete normalization regardless of $\lambda$) and it is assumed $g_{n}=\mathcal{O}(\lambda^{n-1})$ (e.g., at least one direct interaction is required to create a two-body correlation from an uncorrelated state). Using these orderings in Eq. (45) and keeping the lowest order in $\lambda$, Landau’s seminal kinetic equation land:36 can be derived from Eq. (45) by keeping only term (45a) to model $g_{2}$. Solving Eq. (45a) at equilibrium implies

[TABLE]

Using this result in Eq. (46) shows that

[TABLE]

which diverges logarithmically in both the large and small $r$ limits. This indicates that the Landau closure neglects both short range physics [term (45b)] and screening [term (45c)]. Analogously, the same logarithmic divergences would be observed in the terms (45b) and (45c), but these are dropped from the analysis. The same result is obtained for other plasma kinetic equations that are equivalent to the Landau equation, such as Rosenbluth’s Fokker-Planck equation. rose:57 Landau argued that the limits of integration should range from the thermal distance of closest approach in a binary collision, $r_{\textrm{L}}=e^{2}/k_{B}T$, to the Debye length $\lambda_{D}$.

Boltzmann Equation: The Boltzmann equation assumes $\lambda=nl^{3}\ll 1$. In terms of this perturbative parameter, the Liouville operators satisfy $\mathcal{L}_{i}=\mathcal{O}(\lambda^{0})$ and $\mathcal{L}_{ij}=\mathcal{O}(\lambda^{0})$, and, from their dependence on the particle density, the distribution functions $f=\mathcal{O}(\lambda^{1})$ and $g_{n}=\mathcal{O}(\lambda^{n})$. Using these orderings in Eq. (45) and keeping the lowest order in $\lambda$, the Boltzmann equation can be derived from Eq. (45) by keeping terms (45a) and (45b), while dropping terms (45c) and (45d) (this is equivalent to taking a closure $f^{(3)}=0$ and applying the method of Sec. III). ferz:72 Substituting Eq. (47) into the equation resulting from this approximation shows that it implies

[TABLE]

at equilibrium. Using this result in Eq. (46) shows that

[TABLE]

This expression converges in the close interaction limit ($r\rightarrow 0$), but diverges logarithmically in the far interaction limit. Analogously, the same logarithmic divergence would be observed in the term Eq. (45c) that is dropped from the analysis. Physically, this is because screening is neglected along with Eq. (45c). When the Boltzmann equation is used in plasma theory, interactions are limited to within a Debye length $\lambda_{D}$ ad hoc.

Lenard-Balescu Equation: The Lenard-Balescu equation lena:60 ; bale:60 ; guer:60 can be derived from Eq. (45) by keeping terms (45a) and (45c) [while dropping terms (45b) and (45d)]. nich:83 Substituting Eq. (47) into the resultant approximation shows that this implies

[TABLE]

Using this result in Eq. (46) shows that

[TABLE]

This expression converges in the far interaction limit ($r\rightarrow\infty$), but since the integrand is proportional to $1/r$ for small $r$, it diverges logarithmically in the close interaction limit. The neglected term (45b) contains the physics of the close interaction. This limit is typically resolved by truncating the spatial integral at the thermal distance of closest approach in a binary collision, $r_{\textrm{L}}$.

Frieman-Book Equation: Self-consistently convergent kinetic equations have previously been developed. Frieman and Book frie:63 derived such an equation by matching solutions of Eq. (45) in asymptotic limits of close interaction ($r\lesssim r_{\textrm{L}}$), far interaction ($r\gtrsim\lambda_{D}$) and intermediate scale ($r\sim n^{-1/3}$). The solution in each limit is obtained by neglecting the terms described above for the Boltzmann equation, Lenard-Balescu equation and Landau equation, respectively. The resultant expression for $g_{2}^{\textrm{eq}}$ consists of the sum of the Boltzmann equation plus the Lenard-Balescu equation minus the Landau equation. The radial distribution function then has the form

[TABLE]

(see also Refs. shur:64, ; dewi:65, for further discussion). Using this result in Eqs. (46) and (47) shows that

[TABLE]

This expression converges in both the close and far interaction limits. This shows that the divergences in the other theories are associated with neglecting one of the terms (45b)-(45c), rather than the triplet correlation (45d).

Relation with the theory of Sec. III: Since this procedure preserves the exact equilibrium limit, the relationship implied by the second order equation is simply the exact statistical statement hill:60

[TABLE]

However, a closure of Eq. (20) is still required to determine $g(r)$. Since this is an equilibrium quantity, the tools of equilibrium statistical mechanics are available for this task, and accurate approximations have been developed. In principle, any such approximation for $g(r)$, or even experimental data, can serve as the input to this theory.

One example of an accurate closure for Coulomb systems is the combination of the hypernetted chain and Ornstein-Zernike equations hans:06

[TABLE]

where $h(r)\equiv g(r)-1$, and $\hat{h}(k)$ denotes the Fourier transform of $h(r)$. This closed set of equations is known to accurately describe $g(r)$ in plasmas spanning from asymptotically weak coupling well into the strong coupling regime. Further extensions have been provided, via models for the bridge function that extend this even further to near solid-state plasma conditions. tana:86 ; ichi:92 Thus, very accurate methods are available to determine $g(r)$ at equilibrium at essentially any conditions.

At weak coupling, Eq. (57) reduces to the Debye-Hückel limit

[TABLE]

O’Neil and Rostoker showed that this is valid to order $\Gamma^{3}\ln\Gamma$ for $\Gamma\ll 1$. onei:65 None of the approximate theories described above capture Eq. (58). However, the convergent Frieman-Book result is consistent with it to order $\Gamma^{2}$; see Refs. shur:64, ; dewi:65, for further discussion.

Other previous kinetic theories, including Liboff libo:59 and Paquette, paqu:86 have modified the Boltzmann equation in an ad-hoc manner by modeling binary collisions as occurring via the Debye-Hückel potential (rather than the Coulomb potential). This leads to a convergent kinetic equation that accurately models weakly coupled plasmas. Indeed, if $\phi(r)$ is replaced by the Debye-Hückel potential in Eq. (50) the correct result for $g_{2}$ at equilibrium [Eq. (58)] is obtained.

In a sense, these theories, as well as Lenard-Balescu theory, can be understood as variants of the “dressed test particle” concept. rost:64 Differences lie in the details of how interactions are modeled. Lenard-Balescu type approaches model interactions via the correlation of linear fluctuations. This loses the “particle” concept, and thus does not include close collisions, but it does account for dynamics in the dielectric “dressing”, which can be important for super-thermal particles. For instance, the dynamic aspect of the dielectric response results in order-unity contributions to fast particle stopping in a plasma. zwic:99 In contrast, the modified Boltzmann based approaches retain the particle concept, and physics of close collisions. However, the dielectric response in this case is input in an ad hoc manner, via the screened Coulomb potential, which is accurate for sub-thermal particles nich:83 and thus expected to apply to near equilibrium transport processes.

The approach of Sec. III essentially formalizes the latter line of reasoning by enforcing the viewpoint that the exact equilibrium limit should be maintained in the closure to the BBGKY hierarchy. Since this is the Debye-Hückel result at weak coupling, the theory reduces to the result of Liboff libo:59 in that limit. It also extends this line of reasoning by revealing that the appropriate interaction potential is the potential of mean force, which differs substantially from the Debye-Hückel potential at strong Coulomb coupling. Recent work has shown that modeling the potential of mean force via Eq. (57) enables an extension of traditional plasma transport theory well into the strongly coupled regime (typically for $\Gamma\lesssim 20$). baal:13 ; baal:14 ; baal:15 ; dali:14 Similar results have also been obtained in non-equilibrium two-temperature plasmas. seuf:89 ; shaf:17b ; shaf:19

V Conclusions and outlook

The closure presented in this paper makes two advances. First, it provides a formal expansion parameter [Eq. (11)] for the BBGKY hierarchy that enables a self-contained derivation of a plasma kinetic equation that includes both static screening and close interactions. Second, it provides a means to extend traditional plasma theory to strong Coulomb coupling, not only in the collision operator, but also in the equation of state. It does so by ensuring that the equilibrium limit is preserved in the spatial correlation at second order of the hierarchy.

This formalizes and extends the EPT collision operator obtained previously from physical arguments. baal:13 ; baal:14 ; baal:15 It also provides a conceptual basis that may prove useful for extending the theory to other regimes of plasma physics. For example, kinetic theories have been developed to include three-body interactions. onei:65 ; choh:58 A similar analysis could be applied while using $\Delta f^{(4)}=0$ to close the hierarchy, thereby including the static screening response in effective three-body dynamics. It may provide a basis for addressing strongly magnetized plasmas in which the Lorentz force acts at the collision scale. Finally, the approach was motivated by plasmas, which emphasize the need to account for many-body screening self-consistently. However, the theory may also be applied to other types of systems, particularly in the kinetic theory of dense gases. ferz:72

Appendix A Derivation of Equation (18)

Here, it is shown that the mean force given by Eqs. (18) and (19) is equivalent to (II). First, note that for $i\leq n$,

[TABLE]

so that the gradient of the $n-$particle density $\rho^{(n)}(\mathbf{r}^{n})$ defined by Eq. (14) can be written

[TABLE]

Using,

[TABLE]

we find

[TABLE]

which completes the proof.

Appendix B Derivation of Equation (31)

Equation (31) follows from the $n=1$ component of Eq. (II)

[TABLE]

This derivation follows an analogous method presented in Ref. ferz:72, . First, note that

[TABLE]

where the second line follows from the property that $\rho^{(2)}$ is constant under the interchange $1\leftrightarrow 2$. Combining these two equivalent expressions, and applying the substitution $\mathbf{r}\rightarrow-\mathbf{r}$ as the integration variable in Eq. (60b), gives

[TABLE]

Next, observing that

[TABLE]

shows that Eq. (61) can be written as the divergence of a tensor

[TABLE]

Putting Eq. (64) into (59) completes the derivation of Eq. (31).

Acknowledgements.

The authors thank Nathaniel Shaffer and Louis Jose for helpful discussions. This work was supported by the U. S. Air Force Office of Scientific Research under Award No. FA9550-16-1-0221; and by the U. S. Department of Energy, Office of Fusion Energy Sciences, under Award No. DE-SC0016159. The work of J. D. was supported by the US Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U. S. Department of Energy (Contract No. 89233218CNA000001).