Temperature Anisotropy Relaxation of the One-Component Plasma

TL;DR

This study evaluates the relaxation rate of temperature anisotropy in a one-component plasma using molecular dynamics simulations, comparing results with various kinetic theories across different coupling regimes.

Contribution

It provides a comprehensive comparison of kinetic theories with simulations for anisotropy relaxation, including strong coupling effects often neglected by traditional models.

Findings

Dynamic screening has negligible effect on relaxation rate.

Oscillations and delayed relaxation occur at strong coupling, not captured by existing theories.

Simulation results align with weakly coupled theories but reveal limitations at strong coupling.

Abstract

The relaxation rate of a Maxwellian velocity distribution function that has an initially anisotropic temperature $(T_\parallel \neq T_\perp)$ is an important physical process in space and laboratory plasmas. It is also a canonical example of an energy transport process that can be used to test theory. Here, this rate is evaluated using molecular dynamics simulations of the one-component plasma. Results are compared with the predictions of four kinetic theories; two treating the weakly coupled regime (1) the Landau equation, and (2) the Lenard-Balescu equation, and two that attempt to extend the theory into the strongly coupled regime (3) the effective potential theory and (4) the generalized Lenard-Balescu theory. The role of dynamic screening is studied, and is found to have a negligible influence on this transport rate. Oscillations and a delayed relaxation onset in the temperature profiles are observed at strong coupling, which are not described by the kinetic theories.

Equations

$$f=\frac{n}{{\pi }^{3/2}{v}_{T\parallel }{v}_{T\perp }^2}{e}^{-v_{\parallel }^2/v_{T\parallel }^2}{e}^{-v_{\perp }^2/v_{T\perp }^2}$$

$$\frac{d{T}_{\perp }}{dt}=-\frac{1}{2}\frac{d{T}_{\parallel }}{dt}=-\nu (T_{\perp }-T_{\parallel })$$

$$\nu =\frac{m}{n}\int d^{3}v{v}_{\perp }^2\frac{C(f)}{T_{\perp }-T_{\parallel }}=-\frac{1}{2}\frac{m}{n}\int d^{3}v{v}_{\parallel }^2\frac{C(f)}{T_{\perp }-T_{\parallel }}$$

$$A\equiv \frac{T_{\perp }}{T_{\parallel }}-1$$

$$\frac{dA}{dt}=-\nu (3A+2A^2).$$

$$\Gamma \equiv \frac{e^2/a}{k_{B}T},$$

$$C_{\textrm{L}}=-\frac{2\pi e^{4}\ln\Lambda}{m^{2}}\frac{\partial}{\partial\mathbf{v}}\cdot\int d^{3}v^{\prime}\frac{u^{2}\mathcal{I}-\mathbf{u}\mathbf{u}}{u^{3}}\cdot\biggl{[}f(\mathbf{v})\frac{\partial f(\mathbf{v}^{\prime})}{\partial\mathbf{v}^{\prime}}-f(\mathbf{v}^{\prime})\frac{\partial f(\mathbf{v})}{\partial\mathbf{v}}\biggr{]}$$

$$\frac{\nu_{\textrm{L}}}{\bar{\nu}}=\frac{(1+\frac{2}{3}A)^{3/2}}{A^{2}}\biggl{[}-3+(A+3)\frac{\arctan(\sqrt{A})}{\sqrt{A}}\biggr{]}\ln\Lambda$$

$$\bar{\nu }=2\sqrt{\frac{\pi }{m}}\frac{n{e}^4}{(k_{B}T{)}^{3/2}}={\omega }_p\sqrt{\frac{3}{4\pi }}{\Gamma }^{3/2}$$

$$A_{1}(t)=A_{o}\exp\biggl{(}-\frac{4}{5}\tilde{t}\biggr{)}$$

$$A_2(t)=\frac{21A_o}{\exp (\frac{4}{5}\tilde{t})(17A_o+21)-17A_o}.$$

$$C_{\textrm{LB}}=-\frac{2e^{4}}{m^{2}}\frac{\partial}{\partial\mathbf{v}}\cdot\int d^{3}v^{\prime}\int d^{3}k\frac{\delta(\mathbf{k}\cdot\mathbf{u})}{|\hat{\varepsilon}(\mathbf{k},\mathbf{k}\cdot\mathbf{v})|^{2}}\frac{\mathbf{k}\mathbf{k}}{k^{4}}\cdot\biggl{[}f(\mathbf{v})\frac{\partial f(\mathbf{v}^{\prime})}{\partial\mathbf{v}^{\prime}}-f(\mathbf{v}^{\prime})\frac{\partial f(\mathbf{v})}{\partial\mathbf{v}}\biggr{]}.$$

$$\nu_{\textrm{LB}}=\frac{ne^{4}}{\pi}\int d^{3}k\int_{-\infty}^{\infty}d\omega\frac{1}{k^{4}}\frac{1}{|\hat{\varepsilon}(\mathbf{k},\omega)|^{2}}\frac{k_{\parallel}^{2}k_{\perp}^{2}}{(k_{\parallel}^{2}k_{B}T_{\parallel}+k_{\perp}^{2}k_{B}T_{\perp})^{2}}\exp\biggl{(}-\frac{m\omega^{2}}{k_{\parallel}^{2}k_{B}T_{\parallel}+k_{\perp}^{2}k_{B}T_{\perp}}\biggr{)}$$

$$\hat{\varepsilon}(\mathbf{k},\omega)=1+\frac{\omega_{p}^{2}}{k^{2}}\int d^{3}v\frac{\mathbf{k}\cdot\partial f/\partial\mathbf{v}}{\omega-\mathbf{k}\cdot\mathbf{v}}=1-\frac{\omega_{p}^{2}}{k_{\parallel}^{2}v_{T\parallel}^{2}+k_{\perp}^{2}v_{T\perp}^{2}}Z^{\prime}\biggl{(}\frac{\omega}{\sqrt{k_{\parallel}^{2}v_{T\parallel}^{2}+k_{\perp}^{2}v_{T\perp}^{2}}}\biggr{)}.$$

$$\frac{\nu_{\textrm{LB}}}{\bar{\nu}}=\biggl{(}\frac{1+\frac{2}{3}A}{1+A}\biggr{)}^{3/2}\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}dxe^{-x^{2}}\int_{-1}^{1}dy\frac{y^{2}(1-y^{2})}{(y^{2}\tilde{A}+1)^{3/2}}\int_{0}^{\bar{k}_{\max}}\frac{d\bar{k}}{\bar{k}}\biggl{|}1+\frac{\Phi(x)}{\bar{k}^{2}(y^{2}\tilde{A}+1)}\biggr{|}^{-2}.$$

$$\frac{{\nu }_{\text{LB}}}{\bar{\nu }}$$

$$\displaystyle-\frac{1}{2}\frac{\Phi_{R}}{\Phi_{I}}\biggl{[}\arctan\biggl{(}\frac{\bar{k}_{\max}^{2}(y^{2}\tilde{A}+1)+\Phi_{R}}{\Phi_{I}}\biggr{)}-\arctan\biggl{(}\frac{\Phi_{R}}{\Phi_{I}}\biggr{)}\biggr{]}\biggr{\}},$$

$${\bar{k}}_{\max }=\frac{{\lambda }_{D\perp }}{b_{\min }}=\frac{\sqrt{k_B{T}_{\perp }/(4\pi e^{2}n)}}{e^2/(2k_{B}T)}=\sqrt{\frac{1+A}{1+\frac{2}{3}A}}\frac{2}{\sqrt{3}{\Gamma }^{3/2}},$$

$$2\int_{0}^{1}dy\frac{y^{2}(1-y^{2})}{(y^{2}\tilde{A}+1)^{3/2}}=\frac{1}{r^{3/2}}\frac{1}{A^{2}}\biggl{[}-3+(A+3)\frac{\arctan(\sqrt{A})}{\sqrt{A}}\biggr{]}.$$

$$\frac{\nu_{\textrm{LB,I}}}{\bar{\nu}}=\frac{(1+\frac{2}{3}A)^{3/2}}{A^{2}}\biggl{[}-3+(A+3)\frac{\arctan(\sqrt{A})}{\sqrt{A}}\biggr{]}\Xi_{\textrm{LB,I}}$$

$$\Xi_{\textrm{LB,I}}=\frac{1}{2}\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}dxe^{-x^{2}}\biggl{\{}\ln\biggl{(}\frac{|\bar{k}_{\max}^{2}+\Phi|}{|\Phi|}\biggr{)}-\frac{\Phi_{R}}{\Phi_{I}}\biggl{[}\arctan\biggl{(}\frac{\bar{k}_{\max}^{2}+\Phi_{R}}{\Phi_{I}}\biggr{)}-\arctan\biggl{(}\frac{\Phi_{R}}{\Phi_{I}}\biggr{)}\biggr{]}\biggr{\}}.$$

$$\frac{{\nu }_{\text{LB,s}}}{\bar{\nu }}$$

$$\displaystyle-\frac{1}{2}\frac{\bar{k}_{\max}^{2}(y^{2}(r-1)+1)}{1+\bar{k}_{\max}^{2}(y^{2}(r-1)+1)}\biggr{]}.$$

$$\frac{\nu_{\textrm{DH}}}{\bar{\nu}}=\frac{(1+\frac{2}{3}A)^{3/2}}{A^{2}}\biggl{[}-3+(A+3)\frac{\arctan(\sqrt{A})}{\sqrt{A}}\biggr{]}\Xi_{\textrm{DH}}$$

$$\Xi_{\textrm{DH}}=\ln\bigl{(}\sqrt{\Lambda^{2}+1}\bigr{)}-\frac{1}{2}\frac{\Lambda^{2}}{\Lambda^{2}+1}.$$

$$C_{\text{B}}=\int d^3{v}^{\prime }\int d\Omega \sigma u[f(\hat{\mathbf{v}})f({\hat{\mathbf{v}}}^{\prime })-f(\mathbf{v})f({\mathbf{v}}^{\prime })],$$

$$\frac{d{T}_{\parallel }}{dt}=\frac{m}{n}\int d^{3}v{v}_{\parallel }^2{C}_B(f)=\frac{m}{n}\int d^{3}v\int d^3{v}^{\prime }\{\Delta (v_{\parallel }^2)\}f(\mathbf{v})f({\mathbf{v}}^{\prime })$$

$$\frac{dT_{\parallel}}{dt}=\frac{n}{\pi^{3/2}v_{T\perp}^{2}v_{T\parallel}}\frac{m}{2^{9/2}}\int d^{3}u\,u(u^{2}-3u_{\parallel}^{2})\bar{\sigma}^{(2)}(u)\exp\biggl{(}-\frac{1}{2}\frac{u_{\perp}^{2}}{v_{T\perp}^{2}}\biggr{)}\exp\biggl{(}-\frac{1}{2}\frac{u_{\parallel}^{2}}{v_{T\parallel}^{2}}\biggr{)}$$

$${\bar{\sigma }}^{(l)}=2\pi {\int }_0^{\infty }dbb[1-{\cos }^l(\pi -2\Theta )]$$

$$\Theta =b{\int }_{r_o}^{\infty }dr{r}^{-2}[1-b^2/r^2-(e\varphi (r)/k_{B}T)/{\xi }^2{]}^{-1/2}.$$

$$\frac{\nu_{\textrm{EPT}}}{\bar{\nu}}=\chi\frac{3\sqrt{\pi}}{16}\frac{(1+\frac{2}{3}A)^{3/2}}{\sqrt{\alpha}A^{5/2}}\int_{0}^{\infty}d\xi\,\xi^{2}e^{-\alpha\xi^{2}}\frac{\bar{\sigma}^{(2)}}{\sigma_{o}}\biggl{[}\frac{2}{3}\xi^{2}\alpha A\textrm{erf}(\xi\sqrt{\alpha A})-\psi(\xi^{2}\alpha A)\biggr{]}.$$

Full text

1 Introduction

The collisional relaxation of a temperature anisotropy is a canonical example of energy transport in plasmas [1]. It serves as a test for transport theories, as well as playing an essential role in many experiments [2, 3, 4, 5] and natural occurring plasmas [6]. For instance, plasmas that are preferentially heated, or cooled, in one direction will form a temperature anisotropy. Magnetized plasmas often have different energy confinement times either along or against the magnetic field, and can form a temperature anisotropy as a result [7]. It is often essential to know the rate at which this temperature anisotropy relaxes in order to accurately model these plasmas. Here, we use molecular dynamics (MD) simulations to test theories of the collisional temperature anisotropy relaxation rate. The work concentrates on the unmagnetized one-component plasma (OCP) [8] and considers a broad range of coupling parameters. Although this concentrates on unmagnetized plasmas, the results are also expected to apply to magnetized plasmas as long as $\omega_{c}/\omega_{p}\ll 1$, where $\omega_{c}=eB/m$ is the gyrofrequency and $\omega_{p}=\sqrt{4\pi e^{2}n/m}$ is the plasma frequency [1, 9]. Examples of systems in which moderate to strong coupling arises and magnetization may lead to temperature anisotropies include magnetized ultracold neutral plasmas [10] and in magnetized inertial confinement fusion [11].

The initial distribution function is assumed to have the form

[TABLE]

where $n$ is the number density, $v_{T\parallel}^{2}=2k_{B}T_{\parallel}/m$ and $v_{T\perp}^{2}=2k_{B}T_{\perp}/m$. Here, $T_{\parallel}$ and $T_{\perp}$ are the parallel and perpendicular temperatures, which are related to the total temperature by $T=(T_{\parallel}+2T_{\perp})/3$. In weak and moderately coupled plasmas, the rate of relaxation of the temperature anisotropy can be modeled from the $v_{\parallel}^{2}$ or $v_{\perp}^{2}$ moment of a kinetic equation $df/dt=C(f)$. This provides the anisotropy evolution equation

[TABLE]

where

[TABLE]

is the collisional equipartition rate and $C(f)$ is the collision operator.

A consequence of the underlying assumptions of kinetic theories of weak or moderately coupled plasmas is that kinetic energy is a conserved quantity throughout the relaxation $dT/dt=0$. However, this assumption becomes suspect at strong coupling where the potential energy of interactions in the system is high enough that the exchange between kinetic and potential energy may significantly influence the relaxation processes. The MD simulation results shown here reveal that, indeed, kinetic energy is well conserved at weak and moderate coupling, but that small dynamic oscillations in the kinetic energy, due to exchange with potential energy, are observed at strong coupling. Recent MD simulations have also shown that external energy perturbations can lead to the spontaneous development of temperature anisotropies [7].

Each of the theories that will be discussed assume kinetic energy conservation, and thus satisfy Eq. (2). This is a set of two coupled ordinary differential equations that can be solved for $T_{\parallel}(t)$ and $T_{\perp}(t)$. For comparing different collision theories, it is convenient to note that the equations can be decoupled to a single equation by defining

[TABLE]

as a measure of the temperature anisotropy. Note that $T>0$ implies $A>-1$. In terms of $A$, Eq. (2) is

[TABLE]

Thus, $A(t)$ along with $dT/dt=0$ provides the solution for both the parallel $T_{\parallel}(t)=3T/(3+2A(t))$ and perpendicular $T_{\perp}(t)=3T(1+A(t))/(3+2A(t))$ temperatures.

Here, we compute $\nu$ from four different collision operators and compare the results with MD simulations. Simulations were carried out spanning a wide range of coupling strengths. The coupling strength of the OCP is described by the Coulomb coupling parameter [8]

[TABLE]

where $e$ is the elementary electron charge, $T$ is the total temperature and $a=(3/4\pi n)^{1/3}$ is the average interparticle spacing. Two of the collision operators we compare with apply to weakly coupled plasmas ($\Gamma\ll 1$): the Landau equation [12] and the Lenard-Balescu equation [13]. The distinguishing feature of these is that the Lenard-Balescu equation treats dynamic (velocity-dependent) screening, whereas the Landau equation is based on a static screening model. We find that both give very similar predictions for $\Gamma\lesssim 0.1$ and that both models rapidly break down for larger $\Gamma$ values. Thus, there is no significant contribution associated with dynamic screening for this process over the range of coupling strength in which these theories are applicable. Physically, this is because temperature anisotropy relaxation is a thermal processes associated with particles in the thermal part of the distribution, which do not have significant velocity-dependent components to their screening clouds. We also compare with two theories that extend these collision operators to strong coupling: effective potential theory (EPT) [14, 15, 16] and generalized Lenard-Balescu theory (GLB) [17, 18, 19]. Again, GLB includes dynamic effects in the long distance component of the linear response function, whereas EPT is based on a static description. We find that each of these give similar results that are in agreement with MD from weak to moderate coupling, but that EPT is able to extend to significantly larger $\Gamma$ values than the GLB theory. Again, the conclusion is that dynamic screening does not significantly influence the anisotropy relaxation rate.

The main approximation in this analysis is that the distribution function maintains the anisotropic Maxwellian form from Eq. (1) throughout the relaxation, with only the magnitudes of $T_{\parallel}$ and $T_{\perp}$ evolving. It is known from other transport processes, such as self-diffusion of the OCP, that deviations from Maxwellian distributions can contribute corrections on the order of $20\%$ in the weakly coupled limit, but are typically insignificant in the strongly coupled limit [20, 15]. To the best of our knowledge, a perturbation expansion of the distribution function, such as the Chapman-Enskog [20, 21] or Grad [22] expansions for hydrodynamic transport coefficients, has not yet been developed for temperature anisotropy relaxation. This aspect of the analysis will be discussed further in Sec. 5.

In addition to comparing with MD simulations, explicit expressions are derived for $\nu$ from each collision operator. Various limits of each are also considered, including early time, weak anisotropy and limits connecting the theories based on static or dynamic screening. The results quantify the accuracy of each approach, and the various reduced models in asymptotic limits. This provides useful information for modeling a variety of different plasmas with temperature anisotropies.

Finally, time-dependent temperature profiles from the MD simulations are used to highlight interesting features of correlation effects that arise at strong coupling, which are not described by any of the theories considered. These include delays in the temperature relaxation at early times, substantial oscillations in the temperature profiles in time, and small violations of the kinetic energy conservation assumption ($dT/dt=0$), indicating correlation between kinetic and internal energies during the relaxation process. These aspects are discussed in Sec. 5.

2 Weakly Coupled Plasma Theories

2.1 Landau theory

Plasma transport theory is typically derived from the Landau kinetic equation [12], which for the OCP is

[TABLE]

where $\mathbf{u}\equiv\mathbf{v}-\mathbf{v}^{\prime}$ and $\ln\Lambda$ is the Coulomb logarithm. In Landau’s theory $\Lambda=b_{\textrm{max}}/b_{\textrm{min}}=\lambda_{D}/r_{L}$ is the ratio of the maximum and minimum length scales over which particles interact. The maximum is determined by the Debye screening cloud $\lambda_{D}=\sqrt{k_{B}T/(4\pi e^{2}n)}$ and the minimum by the distance of closest approach $r_{L}=e^{2}/(2k_{B}T)$. In terms of the Coulomb coupling parameter $\Lambda=2/(\sqrt{3}\Gamma^{3/2})$. The standard plasma physics result for the temperature anisotropy relaxation rate follows from substituting the distribution from Eq. (1) into Eq. (3) with Eq. (7) as the collision operator. This leads to Eq. (2) with the equipartition rate [1, 23]

[TABLE]

where

[TABLE]

is a reference collision rate. Equation (8) is valid for any sign of $A$, but for $A<0$ is often written in the equivalent form by replacing $\arctan(\sqrt{A})/\sqrt{A}=\textrm{arctanh}(\sqrt{-A})/\sqrt{-A}$ [23].

Equation (8) was derived from Eq. (7) without explicit reference to the magnitude of the anisotropy. In the weak anisotropy limit $A\ll 1$, Eqs. (5) and (8) reduce to $dA/d\tilde{t}=-\frac{4}{5}A-\frac{68}{105}A^{2}+\mathcal{O}(A^{3})$, where $\tilde{t}\equiv t\bar{\nu}\ln\Lambda$. Thus, to first order in the anisotropy expansion the relaxation is exponential

[TABLE]

and at second order has the form

[TABLE]

Figure (1) shows results of the temperature evolution for initial anisotropies of $A_{o}=0.2$ and $A_{o}=1$ predicted from Eq. (8). Also shown are the small anisotropy expansions from Eqs. (10) and (11). This shows that for a weak initial anisotropy, the simple expression obtained from the lowest order term in the expansion, Eq. (10), provides an accurate approximation of the full result. As $A$ increases to 1, deviations from the exponential form become important, but these are accurately modeled by Eq. (11).

2.2 Lenard-Balescu theory

An alternative to Landau’s kinetic theory for weakly coupled plasmas is the Lenard-Balescu equation [13], which for the OCP is

[TABLE]

The primary distinction between this and Landau’s collision operator is that the long-range screening is self-consistently accounted for via the linear dielectric response function $\hat{\varepsilon}(\mathbf{k},\mathbf{k}\cdot\mathbf{v})$. In contrast to the static Debye-Hückel screening assumed in the Landau equation, this is a generalized theory that can, in principle, account for velocity dependence (aka dynamics) of the screening. Following the same procedure as the previous section by using the distribution from Eq. (1) in Eq. (3) with Eq. (12) as the collision operator provides

[TABLE]

where

[TABLE]

The expression following the second equality in Eq. (14) follows from inserting the anisotropic Maxwellian distribution from Eq. (1), $Z$ is the plasma dispersion function [24] and the “prime” denotes the first derivative with respect to the argument.

Lenard-Balescu theory contains an ultraviolet divergence because it does not self-consistently account for the close interaction of colliding particles. This physics is imposed by limiting the minimum interaction length scale to be the classical distance of closest approach, or equivalently, the maximum wavenumber to be $k_{\textrm{max}}=1/b_{\textrm{min}}=2k_{B}T/e^{2}$. Representing $k$ in a spherical coordinate system $\mathbf{k}=k\cos\theta\sin\phi\hat{x}+k\sin\theta\sin\phi\hat{y}+k\cos\phi\hat{z}$, Eq. (13) can be expressed as

[TABLE]

Here, the following variables have been defined: $y=\cos\phi$, $\bar{k}=k\lambda_{D\perp}$, $x=(\omega/\omega_{p})/\sqrt{\bar{k}_{\parallel}^{2}r+\bar{k}_{\perp}^{2}}$, $r=T_{\parallel}/T_{\perp}=1/(1+A)$, $\tilde{A}\equiv-A/(1+A)$, and $\Phi(x)\equiv-\frac{1}{2}Z^{\prime}(x/\sqrt{2})$, where $\lambda_{D\perp}^{2}\equiv k_{B}T_{\perp}/(4\pi e^{2}n)$ is a Debye length associated with the perpendicular temperature. Equation (15) can be further reduced by carrying out the $\bar{k}$ integral, leading to

[TABLE]

where

[TABLE]

$\Phi_{R}=\textrm{Re}\{\Phi\}$ and $\Phi_{I}=\textrm{Im}\{\Phi\}$. Equation (16) follows directly from the collision operator, Eq. (12), assuming only that the distribution function remains an anisotropic Maxwellian throughout the evolution. This will be compared directly with the Landau equation result from Eq. (8), but first the following subsections discuss how the integrals in Eq. (16) can be simplified via accurate approximations.

2.2.1 Isotropic Screening Approximation

Note that the wavenumber ($\bar{k}$) and azimuthal ($y$) integrals in Eq. (15) would be decoupled if not for the anisotropy of the screening (i.e., take $\tilde{A}=0$ in the last term). This feature arises in the dielectric function of Eq. (14) because the temperature anisotropy leads to a different screening length in the parallel and perpendicular directions. However, the distance of closest approach, which is also temperature dependent, must be added in an ad-hoc way in Lenard-Balescu theory, so it is not clear what the overall anisotropy dependence of this term should be. In any case, since the dependence of transport rates on the screening length is logarithmic in the weakly coupled limit, it has a negligible influence if the coupling parameter $\Gamma$ is sufficiently small. If the screening is assumed isotropic, one can connect the results from the Landau equation based theory, Eq. (8), with the results from the Lenard-Balescu equation based theory via the Coulomb logarithm. Taking $\tilde{A}=0$ in the screening term of Eq. (15), or equivalently in the terms in curly braces in Eq. (16), we note that

[TABLE]

In this isotropic screening limit, the Lenard-Balescu result from Eq. (16) can be written in the same form as the Landau result from Eq. (8)

[TABLE]

but where the Coulomb logarithm is replaced by

[TABLE]

This is an effective Coulomb logarithm that includes dynamic screening (via the $\Phi$ terms), as well as an $\mathcal{O}(1)$ term in addition to the logarithmic term. Section 2.3 will show that the latter contribution is significant only when $\Gamma\gtrsim 0.1$, which also happens to be where the theory breaks down due to the onset of strong correlations.

2.2.2 Static Screening Limit

The primary difference between the Landau and Lenard-Balescu based theories is the role of dynamic screening. To quantify the role of dynamic screening within the Lenard-Balescu theory, Eq. (16), we compare with the static screening limit. In this limit, $\Phi(x)=\Phi(0)$, where $\Phi_{R}(0)=1$, $\Phi_{I}(0)=0$ and the $x$ integral in Eq. (16) can be evaluated explicitly, leading to

[TABLE]

This is compared with the full expression from Eq. (16) in Fig. 2 (discussion in Sec. 2.3).

2.2.3 Debye-Hückel Screening

Finally, to connect the Lenard-Balescu and Landau results, the Deybe-Hückel screening limit is considered. The linear dielectric corresponding to Debye-Hückel screening is $\hat{\varepsilon}=1+1/(k\lambda_{D})^{2}=1+3/(\bar{k}^{2}(2+r))$. Using this in Eq. (13) leads to the same form for the relaxation rate as Eqs. (8) and (19),

[TABLE]

but where

[TABLE]

Note that for $\Lambda\gg 1$, $\Xi_{\textrm{DH}}\rightarrow\ln\Lambda+\mathcal{O}(1)$, returning the conventional result from Eq. (8).

2.3 Comparison of weakly coupled plasma theories

Figure 2 shows a comparison of the results of the weakly coupled theories described in this section. Panel (a) shows that both the Landau result from Eq. (8) and the Lenard-Balescu result from Eq. (16) give essentially indistinguishable predictions for $\Gamma\lesssim 0.1$ at $A=-0.2$ (a weak temperature anisotropy $T_{\perp}/T_{\parallel}=0.8$). This is because dynamic screening leads to a $\mathcal{O}(1)$ correction to the Coulomb logarithm; only super-thermal particles exhibit significant dynamic screening, but this transport process is controlled by thermal particles. Since $\ln\Lambda$ becomes large in the weakly coupled regime, any contribution from dynamic screening is inconsequential. For $\Gamma\gtrsim 0.1$ the two results differ, and the predictions of the Landau theory becomes negative near $\Gamma\simeq 1$. However, we will see in the next section that neither of these theories are accurate for $\Gamma\gtrsim 0.1$ because correlation effects onset, which require a strongly coupled plasma theory to describe. Thus, there seems to be no advantage to using the comparatively complicated Lenard-Balescu theory to describe this process. Figure 2b shows the ratio of each of the approximations of the Lenard-Balescu result with the full numerical result. The results of the static screening model, and the Debye-Hückel screening model are within a few percent of the full results for $\Gamma\lesssim 0.1$; further emphasizing the conclusion that dynamic screening does not significantly influence this process. The isotropic screening assumption is excellent for $\Gamma\lesssim 1$, beyond which the Lenard-Balescu theory is not expected to be accurate. The strongly coupled regime requires a theory that can account for correlation effects, which is the topic of the next section.

Figure 2 shows that the influence of dynamic screening is negligible when the anisotropy is small. However, it may be expected to become more important as the anisotropy increases. Figure 3 shows how the relaxation rate predicted by each theory depends on the anisotropy parameter $A$. The relaxation rate is plotted in units of $\bar{\nu}\ln\Lambda$, so the Landau result from Eq. (8) provides a single curve. The Taylor expansion of this curve up to linear order, $\nu/(\bar{\nu}\ln\Lambda)=4/15+(4/105)A$, shows that, as expected, the expansion holds in a narrow neighborhood near $A=0$. The accuracy of the expansion extends much further on the $A>0$ side than for $A<0$. The dashed lines show the predictions of the Lenard-Balescu result from Eq. (16) at $\Gamma=10^{-1}$, $10^{-2}$, $10^{-4}$ and $10^{-6}$. Comparison of these curves shows that, indeed, dynamic screening has a larger influence as the anisotropy increases or the coupling strength increases. For example, at $A=10$ and $\Gamma=10^{-2}$ there is an approximately 10% difference between the Landau and Lenard-Balescu predictions. Two notes should be kept in mind when interpreting this result.

First, the Lenard-Balescu theory treats the long-distance behavior of Coulomb interactions, which corrects the infrared divergence of Landau theory and incorporates temperature anisotropy in the screening length. However, it still contains an ultraviolet divergence, which is resolved by limiting the distance of closest approach to be the thermal Landau length $r_{L}=e^{2}/(2k_{B}T)$. The temperature arising in this length scale is assumed to be the total temperature, but it is not clear if it should also account for temperature anisotropy.

Second, the potential role of dynamic screening should be analyzed keeping in mind that both theories are making an assumption that the distribution function remains an anisotropic Maxwellian throughout the evolution. Expectations from other transport processes, such as diffusion, suggest that the error in this approximation is on the order of 20% at weak coupling [20, 15]. Some evidence for this from comparison with MD simulations will also be shown in Sec. 5. Thus, the error associated with approximations related to the distribution function are often more significant than changes in the transport rates associated with anisotropy in the screening.

3 Strongly Coupled Plasma Theories

Both the Landau and Lenard-Balescu based theories from the previous section address the weakly coupled limit. In this section, two theories are evaluated that attempt to extend these methods into the strongly coupled regime. Each of these is based on phenomenological assumptions, rather than rigorous asymptotic expansions, so the subsequent comparison with MD serves as a test of these assumptions.

3.1 Effective Potential Theory

Like the Landau approach, the effective potential theory (EPT) [14, 15, 16] is based on the Boltzmann collision operator, which for the OCP is

[TABLE]

where $\sigma$ is the differential scattering cross section, the hat denotes the velocity vector after a binary collision: $\hat{\mathbf{v}}=\mathbf{v}+\Delta\mathbf{v}$ and $\mathbf{u}\equiv\mathbf{v}-\mathbf{v}^{\prime}$. Unlike the other theories discussed in this paper, the anisotropy relaxation rate has not previously been calculated from this theory, so we provide details of the calculation here. Starting with the $v_{\parallel}^{2}$ moment of Eq. (24) provides

[TABLE]

where the curly brackets are defined as $\{\chi\}=\int d\Omega\sigma u\Delta\chi$ [25], for any function $\chi$, $\Delta(v_{\parallel}^{2})=\hat{v}_{\parallel}^{2}-v_{\parallel}^{2}$, and momentum conservation implies $\Delta\mathbf{v}=\Delta\mathbf{u}/2$, so $\Delta(v_{\parallel}^{2})=v_{\parallel}\Delta u_{\parallel}+\frac{1}{4}(\Delta u_{\parallel})^{2}$.

Assuming that the distribution function has the anisotropic Maxwellian form described by Eq. (1), the three velocity integrals for the variable $\mathbf{v}$ in Eq. (25) can be computed analytically leading to

[TABLE]

in which

[TABLE]

is the $l^{\textrm{th}}$ momentum scattering cross section and $\theta=\pi-2\Theta$ is the scattering angle, where

[TABLE]

Obtaining Eq. (26) from Eq. (25) makes use of the relations $\{\Delta\mathbf{u}\}=-u\mathbf{u}\bar{\sigma}^{(1)}(u)$ and $\{\Delta\mathbf{u}\Delta\mathbf{u}\}=u[2\mathbf{u}\mathbf{u}\bar{\sigma}^{(1)}+\frac{1}{2}(u^{2}\mathcal{I}-3\mathbf{u}\mathbf{u})\bar{\sigma}^{(2)}(u)]$, which imply $u_{\parallel}\{\Delta u_{\parallel}\}+\frac{1}{2}\{\Delta(u_{\parallel}^{2})\}=\frac{1}{4}u(u^{2}-3u_{\parallel}^{2})\bar{\sigma}^{(2)}$. Applying a spherical coordinate system $\mathbf{u}=u[\cos\theta^{\prime}\,\sin\phi^{\prime}\hat{x}+\sin\theta^{\prime}\sin\phi^{\prime}\hat{y}+\cos\phi^{\prime}\hat{z}]$ both the $\theta^{\prime}$ and $\phi^{\prime}$ integrals can be carried out. This leads to Eq. (2) where the anisotropy relaxation rate is

[TABLE]

Here, $\sigma_{o}\equiv\pi e^{4}/(2k_{B}T)^{2}$, $\xi^{2}\equiv\frac{1}{2}u^{2}/v_{T}^{2}=u^{2}/(4k_{B}T/m)$, $\alpha\equiv T/T_{\perp}=\frac{1}{3}(3+A)/(1+A)$ and

[TABLE]

is the Maxwell integral. The right side of Eq. (29) has been multiplied by a factor $\chi$, which is an aspect of EPT that will be explained below.

The expression in Eq. (29) follows directly from inserting the anisotropic distribution function in the Boltzmann collision operator. The EPT extends this Boltzmann-based result by relaxing (but not eliminating) two of its underlying assumptions: binary collisions and molecular chaos. First, the main premise of the theory is that binary collisions between particles do not occur in isolation, but rather in a “sea” of background particles. This “relaxes” Boltzmann’s assumption by accounting for aspects of many-body physics via the interaction potential. The appropriate potential is one which takes into account the fields associated with nearby charges. In particular, the potential of mean force is the potential obtained when taking two particles at fixed positions and averaging over the positions of all other particles [26]

[TABLE]

Here $\mathcal{Z}=\int\exp(-U/k_{B}T)d\mathbf{r}_{N}$ is the configurational integral and $U\equiv\sum_{i,j}v(|\mathbf{r}_{i}-\mathbf{r}_{j}|)$. This depends only on the bare interaction potential of the particles $v_{ij}$, which is the Coulomb potential. Thus, once $g(r)$ is determined, the interaction potential is modeled as $\phi(r)=-k_{B}T\ln[g(r)]$.

Second, the factor $\chi$ in Eq. (29) accounts for an increased collision rate due to the excluded volume surrounding particles interacting via a repulsive (Coulomb) potential [16]. This excluded volume is sometimes referred to as the Coulomb hole. It decreases the volume of space that particles can occupy, resulting in an increased scattering rate. Accounting for this relaxes Boltzmann’s molecular chaos assumption somewhat because it introduces an aspect of correlation in the initial positions of scattering particles. We model this factor using a modified version of Enskog’s theory of hard sphere gases that has been adapted to plasmas [16]. In this theory, the collision frequency enhancement factor is

[TABLE]

where $b\rho\simeq\pi n\bar{\sigma}^{3}/3$ and $\bar{\sigma}$ is the effective particle diameter. This is determined from the location where $g(r=\bar{\sigma})=0.87$. For the OCP at $\Gamma\lesssim 1$, $\chi\simeq 1$, while for $\Gamma\simeq 1-100$, $\chi\simeq 1.2-1.4$. This 20-40% increase in the collision frequency can be seen in figure 7 of [16].

Equations (27)–(29) along with $\phi=-k_{B}T\ln[g(r)]$ and (32) provide a closed set of equations based on one input: $g(r)$. To compute $g(r)$ we use the hypernetted chain (HNC) approximation [27]

[TABLE]

where $v(\mathbf{r})/k_{B}T=\Gamma a/r$ is the bare Coulomb potential, $h(\mathbf{r})\equiv g(\mathbf{r})-1$ and “hats” denote Fourier transforms in the spatial coordinate. In Eq. (33), $b(\mathbf{r})$ is a bridge function that that increases the accuracy of the HNC approximation at high $\Gamma$. For this term, we apply the model of Iyetomi et al[28]. The HNC approximation has long been benchmarked for the Coulomb OCP, and it has been shown to provide accurate input for the EPT over the range of coupling strengths considered here [14, 15].

Finally, we note that by modeling the interaction potential using an equilibrium theory, we have implicitly assumed that the effective interaction potential is spherically symmetric. For a large enough temperature anisotropy this approximation should be expected to break down. Next, we consider the weak anisotropy limit and the weakly coupled limit of the EPT theory.

3.1.1 Weak anisotropy expansion

If the temperature anisotropy is small, $|A|\ll 1$, expanding Eq. (29) gives

[TABLE]

in which the generalized Coulomb logarithms

[TABLE]

have been identified in the analogous way as Refs. [25, 14]. Equation (34) returns the same result as the small anisotropy expansions of the Landau theory from Eq. (8) if the appropriate Coulomb logarithm is substituted for $\ln\Lambda$. This can be obtained by applying the weakly coupled limit of the generalized Coulomb logarithms from [15], $\Xi^{(l,k)}=l\Gamma(k)\ln\Lambda$, where $\Gamma(k)$ is the Gamma function.

3.1.2 Weakly coupled limit

To connect Eq. (29) with the standard result from weakly coupled plasma theory, Eq. (8), consider taking the weak coupling limit of Eq. (29). In this limit, $\bar{\sigma}^{(l)}=\lambda_{D}^{2}4\pi l\ln(\Lambda\xi^{2})/(\Lambda\xi^{2})^{2}$ [15], and the $\xi$ dependence of the Coulomb logarithm is a small correction because $\Lambda$ is large, so $\bar{\sigma}^{(l)}/\sigma_{o}\simeq 4l\ln\Lambda/\xi^{4}$. With this, Eq. (29) is

[TABLE]

The first of these integrals can be written in terms of the $\arctan$ function $\int_{0}^{\infty}dx\,e^{-x^{2}}\textrm{erf}(x\sqrt{A})=\arctan(\sqrt{A})/\sqrt{\pi}$. The second integral can be evaluated by first integrating by parts, then identifying a similar integral that can be evaluated in terms of $\arctan$. The result returns Eq. (8). In the weak anisotropy limit $|A|\ll 1$, $\nu/\bar{\nu}\rightarrow[4/15+(4/105)A+\mathcal{O}(A^{3/2})]\ln\Lambda$. Alternatively, this can be obtained from Eq. (34) by applying the weakly coupled limit of the generalized Coulomb logarithms.

3.2 Generalized Lenard-Balescu Theory

The second theory that we consider is Ichimaru’s generalization of the Lenard-Balescu equation [18]

[TABLE]

Here, $G(k)$ is the local field correction (LFC) which accounts for static correlations. For the OCP [18]

[TABLE]

where

[TABLE]

is the static structure factor. This will be determined by the HNC approximation described in Eq. (33). Inserting the anisotropic distribution function from Eq. (1) leads to the evolution equation Eq. (3) with the relaxation rate given by

[TABLE]

where

[TABLE]

is the dielectric function.

Unlike the Lenard-Balescu equation, Eq. (37) does not suffer an ultraviolet divergence because the local field correction accounts for the close interaction physics. For computational simplicity, one of the $k$ integrals can be computed in an analogous way to what was done in Eq. (15) providing a similar result that includes the LFC

[TABLE]

Limiting values of Eq. (42) that are similar to what is listed for the Lenard-Balescu equation in Secs. 2.2.1–2.2.3 can also be taken, providing similar results but where the LFC enters the expressions for the generalized Coulomb logarithm in place of the cutoff $b_{\min}$. The full expression is compared with the EPT predictions in Sec. 5.

4 MD Simulation Description

The MD simulations were performed as follows. $N$ charged particles interacting through the pure Coulomb interaction in a uniform, neutralizing background were placed in a cubic box of volume $V$. Periodic conditions were imposed on all boundaries. Particle trajectories were determined by solving Newton’s equations of motion with the velocity Verlet integrator [29]. The force on an ion that results from its interaction with the ions in the simulation box and with those in the periodically replicated cells were calculated using the Ewald summation technique. For numerical efficiency, the Ewald sum was calculated with a parallel implementation of the particle-particle-particle-mesh ($\rm P^{3}M$) method that simultaneously provides high resolution for individual encounters combined with rapid, mesh-based, long range force calculations [30]. Time was normalized to the plasma frequency $\omega_{p}$. The integration time step $\delta t$ and the number of particles $N$ given below were chosen in order to conserve energy to $<10^{-5}$, and to ensure high enough collision ability in the simulation cell (the calculations are more demanding at small coupling due to long collision mean-free path). Specifically, for $\Gamma\leq 0.1$, we used $N=10^{5}$ and $\delta t=10^{-3}/\omega_{p}$; for $0.1<\Gamma<1$, $N=50000$ and $\delta t=0.01/\omega_{p}$; for $\Gamma\geq 1$, $N=5000$ and $\delta=0.01/\omega_{p}$.

Each simulation consisted of an equilibration phase of length $t_{\textrm{eq}}=N_{\textrm{eq}}\delta t$ followed by a relaxation phase of length $t_{\textrm{run}}=N_{\textrm{run}}\delta t$, with $t_{\textrm{eq}}=1000/\omega_{p}$ and $t_{\textrm{run}}=300/\omega_{p}$. The initial particle positions at time $t=-t_{\textrm{eq}}$ were assigned randomly in the simulation box, with a small region surrounding each particle excluded to avoid initial explosion. The initial particle velocities ${\bf v}_{i}$ were sampled from a Maxwell-Boltzmann distribution at the desired temperature $T$, i.e. at the desired value of $\Gamma$. During the equilibration phase, velocity scaling [29] was used to maintain the desired $\Gamma$ value. At $t=0$, velocity scaling was turned off and the particle velocities were rescaled to the desired initial parallel and perpendicular temperatures as follows

[TABLE]

The system was then left to evolve freely, i.e. in the microcanonical ensemble, for a duration $t_{\textrm{run}}$. During this period the instantaneous parallel and perpendicular temperatures, defined as follows in terms of the particles kinetic energies

[TABLE]

were monitored. In order to sample the initial statistical distribution function, Eq. (1), and to compare the MD simulations to the predictions of kinetic theories, $N_{\textrm{sim}}=25$ independent runs were performed with different initial conditions and averaged. As illustrated in Fig. 4, the averaging smooths out the fluctuations in the kinetic energies inherent to each microcanonical dynamics.

To obtain the relaxation rate $\nu$, we assume that the MD temperatures evolve according to the rate equations (2). The latter implies $\frac{d}{dt}\Delta T=-3\nu\Delta T$ with $\Delta T=T_{\parallel}-T_{\perp}$. In general, as predicted by kinetic theory, $\nu$ depends on the time $t$ through its dependence upon $T_{\parallel}$ and $T_{\perp}$. If one assumes that the dependence is weak enough that $\nu$ is a constant on a short enough time scale $\Delta t$ beyond an initial time $t_{0}$ , then $\Delta T(t)=\Delta T(t_{0})e^{-3\nu(t-t_{0})}$ on this time scale (more sophisticated treatments that take into account the temperature dependence of $\nu$ are possible (e.g. [31]) but we find that they do not affect the results in any significant manner). The MD simulations confirm that $T_{\parallel}-T_{\perp}$ indeed decays exponentially for coupling parameters $\Gamma\!<\!1$ but only following a short transient period of duration $t_{0}$ of order $\omega_{p}^{-1}$; see Fig. 4b. For $\Gamma\geq 1$, as discussed in the next section, the temperature evolution is never exponential and the notion of relaxation rate defined based on the rate equations (2) does not apply. The initial transient period describes the dependence on initial correlations, which are discarded in the kinetic theories discussed before. In practice, the relaxation rates $\nu$ are obtained by fitting the MD data to the analytical solution $\Delta T(t_{0})e^{-3\nu(t-t_{0})}$ over the time interval $[t_{0},t_{0}+\Delta t]$ with the time $t_{0}$ chosen right after the early transient behavior (the data shown below were obtained with $t_{0}=1/\omega_{p}$) and with $t_{0}+\Delta t$ chosen where the exponential behavior switches slope as a consequence of the dependence of $\nu$ on the time $t$. Figure 4b illustrates this procedure.

5 Results

Figure 5 shows parallel and perpendicular temperature profiles at four values of $\Gamma$ from MD simulations that were initiated with temperatures of $T_{\parallel}/T=1.15$ and $T_{\perp}/T=0.923$ ($A=-0.2$). Also shown are the profiles predicted by the EPT theory of Eq. (29). The comparison shows good agreement with the predicted monotonically decreasing profiles at the lower $\Gamma$ values. As the coupling strength increases the rate of the relaxation remains well modeled, but features emerge in the MD data that are not predicted by the model. One feature is that the profiles from MD become non-monotonic, exhibiting oscillations. Such oscillations are common in time correlation functions at strong coupling. For instance, they can be observed in the velocity autocorrelation function, which when integrated provides the macroscopic diffusion coefficient [32]. Dubin [4] has previously discussed how the temperature anisotropy relaxation rate can also be considered as the time integral of an energy correlation function.

Another interesting feature of the MD data at strong coupling is a short initial delay, compared to the theoretical prediction, before relaxation onsets. This is also common feature of other correlation functions at strong coupling and is associated with initial conditions. For instance, the velocity autocorrelation function at short times behaves as $Z(t)/Z(0)=1-\Omega_{E}^{2}t^{2}$, where $\Omega_{E}$ is the Einstein frequency [27]. For an OCP, $\Omega_{E}=\omega_{p}/\sqrt{3}$. At weak coupling the relaxation time is much longer than $\omega_{p}^{-1}$, so this initial stage is so short that it is inconsequential in the overall relaxation profile. For instance, Fig. 5a shows that the relaxation time is approximately 200$\omega_{p}^{-1}$ at $\Gamma=0.05$. However, at strong coupling the relaxation rate is only a few $\omega_{p}^{-1}$, so the short time rolloff on the $\Omega_{E}^{-1}$ timescale is a more pronounced feature of the overall relaxation profile. The data shown in Fig. 5 suggests that similar arguments can also be applied to the temperature anisotropy relaxation process. All of the theories discussed above are Markovian, and as a consequence do not include information about these initial conditions. They all predict exponential profiles at short times, which leads to a linear-in-time term at the relaxation onset rather than a quadratic one.

Figure 6 shows the results of the relaxation rate obtained from the MD data using the fitting procedure described in Sec. 4. The initial temperature anisotropy for this data set was the same $A=-0.2$ as from the data shown in Fig. 5. Predictions of the four theories described in Sections 2 and 3 are shown for comparison. As expected, the weakly coupled theories compare well with the MD data at sufficiently small coupling strength ($\Gamma\lesssim 0.2$). As described in Sec. 2, there is apparently no significant advantage to accounting for dynamic screening at these conditions. Both strongly coupled theories extend the standard plasma theories into the strongly coupled regime, but the generalized Lenard-Balescu theory from Sec. 3.2 predicts a negative collision frequency for $\Gamma\gtrsim 20$ (negative values are not shown on this log-log scale). This feature of the theory has also been noted in computing the viscosity coefficient in [33]. Although MD simulations were carried out at higher $\Gamma$ values than shown in this figure, the oscillatory nature of the temperature profiles made the fitting procedure unreasonable at sufficiently high coupling strength.

Figures 5a and 6 show that for $\Gamma\lesssim 0.1$ the relaxation rate predicted by all of the theories is slightly larger than what is observed in the MD simulations (for example, by approximately 20% at $\Gamma=0.05$). This is likely associated with the assumption made in the theoretical analysis that the distributions maintain the anisotropic Maxwellian form. It is well known in other transport processes that distortions of the distribution function contributes an order unity correction to the transport rates at weak coupling. This is accounted for in hydrodynamic theories, such as Chapman-Enskog or Grad, through the higher order terms of the Sonine polynomial expansions of the distribution function. However, the magnitude of these terms is known to rapidly diminish as $\Gamma\gtrsim 1$. For instance, it was shown in [15] that for self-diffusion of the OCP the approximately 20% contribution of the second order term of the expansion vanishes for $\Gamma>1$ (see figure 9 of that reference). Figure 6 shows evidence for the same effect in the temperature anisotropy relaxation rate. To the best of our knowledge, such an expansion procedure has not yet been applied to temperature anisotropy relaxation.

Finally, Fig. 7 shows temperature relaxation profiles from MD simulations at $\Gamma=0.1$, 1, 10 and 100 conducted with a large initial temperature anisotropy of $T_{\perp}/T_{\parallel}=10$ ($A=9$). The comparison with the EPT predictions is similar to what was observed at smaller initial anisotropy. The accuracy is, perhaps, slightly less than was observed from the small anisotropy cases in Fig. 5 near $\Gamma\simeq 1$, but the generally good agreement over this entire range of coupling strength provides strong evidence that the theory is robust even at a very large initial anisotropy. The predicted rate compares well even at $\Gamma=100$. Dynamic screening apparently is not significant even at $T_{\perp}/T_{\parallel}=10$. Panels e-h show the total temperature as a function of time for the same simulations. These show that at weak coupling the total temperature, and thus total kinetic energy, does not fluctuate more than approximately 0.1% throughout the evolution of the system. The higher $\Gamma$ values show that the temperature fluctuates in time by as much as 0.5%. Recall that each curve is obtained from 25 independent simulations. The total energy of the system was confirmed to be conserved to better than 1 part in $10^{6}$ throughout each simulation. These oscillations are not numerical artifacts, but rather represent oscillations in the exchange between kinetic and potential energy associated with correlations as the system relaxes. None of the theories discussed address potential energy of the system, and thus do not model this effect.

6 Conclusions

This comparison between two common approaches to plasma kinetic theory and MD simulations has suggested a few general conclusions with regard to the temperature anisotropy relaxation, and has also revealed a few remaining gaps in current understanding. The data comparison suggests that dynamic screening, which is modeled in the Lenard-Balescu and generalized Lenard-Balescu theories, does not significantly influence the relaxation rate. The Boltzmann-based approaches (Landau or EPT) can be accurately applied in most situations of practical interest, rather than Lenard-Balescu based theories which are comparatively difficult to evaluate. Comparison with MD revealed that the EPT provides an accurate approach to modeling the anisotropy relaxation rate over a similar range of coupling strength as has been encountered for other processes, such as diffusion [16]. However, none of the theories capture kinetic energy oscillations (associated with kinetic-potential energy exchange) or early time delay (associated with initial conditions) in the time-dependent temperature profiles that were observed in the MD simulations. These effects are associated with strong correlations. Finally, one additional limitation of current models is the assumption that the distribution maintains the specified anisotropic Maxwellian form throughout the evolution. Comparison with MD reveals that improvements to predicted relaxation rates may be expected from accounting for deviations from this assumed distribution.

{acknowledgement}

This material is based upon work supported by LDRD project 20150520ER at Los Alamos National Laboratory. The work of SDB was also supported by the Air Force Office of Scientific Research under award number FA9550-16-1-0221 and by the National Science Foundation under Grant No. PHY-1453736.