The Barkas Effect in Plasma Transport

TL;DR

This study uncovers a charge-sign dependent asymmetry in plasma collision rates at strong coupling, driven by screening effects, which significantly impacts plasma transport models and experimental interpretations.

Contribution

It introduces the first evidence of the Barkas effect in plasma transport, showing how screening causes charge-sign asymmetry in collision rates at moderate to strong Coulomb coupling.

Findings

Collision rates depend on charge sign at strong coupling.

The Barkas effect can increase electron-ion collision rates by up to 70%.

Traditional models do not predict this asymmetry.

Abstract

Molecular dynamics simulations reveal that a fundamental symmetry of plasma kinetic theory is broken at moderate to strong Coulomb coupling: the collision rate depends on the signs of the colliding charges. This symmetry breaking is analogous to the Barkas effect observed in charged-particle stopping experiments and gives rise to significantly enhanced electron-ion collision rates. It is expected to affect any neutral plasma with moderate to strong Coulomb coupling such as ultracold neutral plasmas (UNP) and the dense plasmas of ICF and laser-matter interaction experiments. The physical mechanism responsible is screening of binary collisions by the correlated plasma medium, which causes an asymmetry in the dynamics of large-angle scattering. Because the effect pertains only to close interactions, it is not predicted by traditional transport models based on cut-off Coulomb collisions or linear dielectric response. A model for the effective screened interaction potential is presented that is suitable for the coupling strengths achieved in UNP experiments. Transport calculations based on this potential and the effective potential kinetic theory agree with simulated relaxation rates and predict that the Barkas effect can cause up to a 70% increase in the electron-ion collision rate at the conditions of present UNP experiments. The influence of the Barkas effect in other transport processes is also considered.

Equations

$${\partial }_t{\mathbf{V}}_e=-{\nu }_{ei}{\mathbf{V}}_e,$$

$$T_s(t)=\frac{m_s}{3k_B{N}_s}\sum _{i=1}^{N_s}|{\mathbf{v}}_s^{(i)}(t)-{\mathbf{V}}_s(t){|}^2,$$

$${\mathbf{V}}_s(t)=\frac{1}{N_s}\sum _{i=1}^{N_s}{\mathbf{v}}_s^{(i)}(t)$$

$${\varphi }_{s{s}^{\prime }}(r)=\frac{q_s{q}_{s^{\prime }}}{r}\times \{\begin{array}{ll}1& q_s{q}_{s^{\prime }}>0\\ 1-e^{-r^2/{\alpha }^2}& q_s{q}_{s^{\prime }}<0\end{array}.$$

$$\tau \approx \frac{2\pi }{{\omega }_{pe}}\sqrt{\frac{3}{Z+1}}(\alpha /a{)}^{3/2},$$

$${\theta }_{s{s}^{\prime }}(b,u)=\pi -2b{\int }_{r_{\min }}^{\infty }\frac{dr}{r^2\sqrt{1-{\Phi }_{s{s}^{\prime }}(r,b,u)}}$$

$${\Phi }_{s{s}^{\prime }}(r,b,u)=\frac{b^2}{r^2}+\frac{2{\varphi }_{s{s}^{\prime }}(r)}{m_{s{s}^{\prime }}{u}^2}$$

$${\Xi }_{s{s}^{\prime }}^{(l,k)}=\frac{1}{2}{\int }_0^{\infty }\frac{{\sigma }_{s{s}^{\prime }}^{(l)}(\xi )}{{\sigma }_{0,s{s}^{\prime }}}{\xi }^{2k+3}{e}^{-{\xi }^2}d\xi ,$$

$${\sigma }_{s{s}^{\prime }}^{(l)}(v)={\int }_0^{\infty }[1-{\cos }^l{\theta }_{s{s}^{\prime }}(v,b)]2\pi bdb$$

$${\hat{h}}_{s{s}^{\prime }}(k)={\hat{c}}_{s{s}^{\prime }}(k)+\sum _{\sigma }{n}_{\sigma }{\hat{c}}_{s\sigma }(k){\hat{h}}_{\sigma s^{\prime }}(k),$$

$$g_{s{s}^{\prime }}(r)=\exp \left[-\frac{{\varphi }_{s{s}^{\prime }}(r)}{k_{B}T}+h_{s{s}^{\prime }}(r)-c_{s{s}^{\prime }}(r)\right],$$

$$g_{s{s}^{\prime }}(r)=\exp \left[-\frac{{\varphi }_{s{s}^{\prime }}(r)}{k_B{T}_{s{s}^{\prime }}}+h_{s{s}^{\prime }}(r)-c_{s{s}^{\prime }}(r)\right]$$

$${\hat{h}}_{s{s}^{\prime }}={\hat{c}}_{s{s}^{\prime }}+\frac{T_{s{s}^{\prime }}}{m_{s{s}^{\prime }}}\sum _{\sigma }{n}_{\sigma }\left(\frac{T_{s\sigma }}{m_s}{\hat{c}}_{s\sigma }{\hat{h}}_{\sigma s^{\prime }}+\frac{T_{\sigma s^{\prime }}}{m_{s^{\prime }}}{\hat{h}}_{s\sigma }{\hat{c}}_{\sigma s^{\prime }}\right)$$

$$T_{s{s}^{\prime }}=\frac{m_{s^{\prime }}{T}_s+m_s{T}_{s^{\prime }}}{m_s+m_{s^{\prime }}}$$

$${\varphi }_{ei}(r)\to -\frac{Z{e}^2}{r}(1-e^{-r/\alpha }),$$

$$S_{ii}=\frac{1}{1-n_i{\hat{c}}_{\mathrm{scr}}}$$

$$S_{ei}=\frac{\sqrt{n_e{n}_i}{\hat{c}}_{ei}}{1-n_e{\hat{c}}_{ee}}{S}_{ii}$$

$$S_{ee}=\frac{1}{1-n_e{\hat{c}}_{ee}}+\frac{n_e{n}_i{\hat{c}}_{ei}^2}{(1-n_e{\hat{c}}_{ee}{)}^2}{S}_{ii}$$

$${\hat{c}}_{\mathrm{scr}}={\hat{c}}_{ii}+\frac{T_e}{T_i}\frac{n_e{\hat{c}}_{ei}^2}{1-n_i{\hat{c}}_{ii}}.$$

$${\hat{c}}_{ei}(k)\approx -\frac{{\hat{\varphi }}_{ei}(k)}{k_B{T}_{ei}}=\pm \frac{1}{\sqrt{n_e{n}_i}}\frac{\sqrt{Z}}{{\lambda }_{De}^2{k}^2}$$

$${\hat{c}}_{ee}(k)\approx -\frac{{\hat{\varphi }}_{ee}(k)}{k_B{T}_e}=-\frac{1}{n_e}\frac{1}{{\lambda }_{De}^2{k}^2},$$

$$S_{ei}(k)\approx \pm \frac{\sqrt{Z}{S}_{ii}(k)}{1+{\lambda }_{De}^2{k}^2},$$

$$w_{ei}(r)$$

$$\approx \mp \frac{Z{e}^2}{r}\cdot \frac{2}{\pi }{\lambda }_{De}^2{\int }_0^{\infty }\frac{k\sin (kr)}{1+{\lambda }_{De}^2{k}^2}{S}_{ii}(k)dk,$$

$$\begin{array}{rl}{\nu }_{ei}=const.\times & \iint {\left|\frac{{\hat{\varphi }}_{ei}(k)}{D(k,\omega )}\right|}^2[1-G_{ei}(k)]\\ & \times \mathrm{Im}{\chi }_e^0(k,\omega )\mathrm{Im}{\chi }_i^0(k,\omega )d\omega d\mathbf{k}\end{array}$$

$$D(k,\omega )=\det \left\{{\delta }_{s{s}^{\prime }}-{\hat{\varphi }}_{s{s}^{\prime }}(k)[1-G_{s{s}^{\prime }}(k)]{\chi }_s^0(k,\omega )\right\}$$

$${\hat{c}}_{s{s}^{\prime }}(k)=-\frac{{\hat{\varphi }}_{s{s}^{\prime }}(k)}{k_B{T}_{s{s}^{\prime }}}[1-G_{s{s}^{\prime }}(k)].$$

$${\mathbf{R}}_{ei}=-m_{ei}\iint u\mathbf{u}{\sigma }_{ei}^{(1)}(u)f_e({\mathbf{v}}_e)f_i({\mathbf{v}}_i)d{\mathbf{v}}_{i}d{\mathbf{v}}_e,$$

$${\mathbf{R}}_{ei}=-m_{ei}{n}_e{\bar{\nu }}_{ei}({\mathbf{V}}_e-{\mathbf{V}}_i),$$

$${\bar{\nu }}_{ei}(\eta )=\frac{16\sqrt{\pi }{Z}^2{e}^4{n}_i}{3m_e{m}_{ei}{\bar{v}}_{ei}^3}{\bar{\Xi }}_{ei}(\eta ),$$

$${\bar{\Xi }}_{ei}(\eta )=\frac{3}{32{\eta }^3}{\int }_0^{\infty }\frac{{\sigma }_{ei}^{(1)}(\xi )}{{\sigma }_{0,ei}}{\xi }^{2}H(\xi ,\eta )d\xi ,$$

$${\bar{\Xi }}_{ei}(\eta )={\Xi }_{ei}^{(1,1)}\left[1+c{\eta }^2+O({\eta }^4)\right],$$

$${\bar{\nu }}_{ei}(0)={\nu }_{ei}\approx \sqrt{\frac{2}{3\pi }}Z{\Gamma }_e^{3/2}{\Xi }_{ei}^{(1,1)}{\omega }_{pe}.$$

$$\begin{array}{rl}{Q}_{ei}=-m_{ei}& \iint u\mathbf{u}\cdot \left({\mathbf{v}}_i+\frac{m_{ei}}{m_i}\mathbf{u}\right){\sigma }_{ei}^{(1)}(u)\\ & \times f_e({\mathbf{v}}_e)f_i({\mathbf{v}}_i)d{\mathbf{v}}_{i}d{\mathbf{v}}_e-{\mathbf{V}}_e\cdot {\mathbf{R}}_{ei}.\end{array}$$

Full text

I Introduction

Traditional plasma theory obeys a fundamental symmetry: Coulomb collision rates and associated transport coefficients (such as diffusivity, conductivity, and thermal relaxation) are independent of the signs of the interacting charges.Spitzer, Jr. (1956) That is, the values of the transport coefficients do not change if the electrons are instead modeled as positrons. The origin of this symmetry is the assumption that interactions occur either through binary Coulomb collisionsMontgomery and Tidman (1964) or via linear response to weak electrostatic fluctuations.Ichimaru (1992) Using molecular dynamics (MD) simulations, we demonstrate that this symmetry is broken in moderately and strongly coupled plasmas. In particular, when $\Gamma_{e}\gtrsim 0.1$ the electron-ion relaxation rate is significantly larger than the positron-ion relaxation rate at the same conditions. Here, $\Gamma_{e}=e^{2}/a_{e}k_{B}T_{e}$ is the ratio of the Coulomb potential energy of two electrons at the mean inter-electron separation $a_{e}=(3/4\pi n_{e})^{1/3}$ to the average kinetic energy per electron. The effect is explained using an extension of the Effective Potential Theory (EPT)Baalrud and Daligault (2013) to charge-neutral plasmas, which quantitatively predicts the simulated relaxation rate. It is expected to influence any moderately or strongly coupled neutral plasma, including ultracold neutral plasmas (UNP)Killian et al. (1999); Chen, Witte, and Roberts (2017) and dense plasmas.Chen et al. (2013); Cho et al. (2016)

This breaking of charge-sign symmetry is found to be a many-body effect associated with how screening influences low-energy, large-angle collisions. Such collisions are rare (and thus negligible) at weak coupling, but they become the dominant type of collision in strongly coupled plasmas. Despite this, the current leading approaches to collisional transport based on linear response theory do not treat such interactions.Daligault and Dimonte (2009); Vorberger et al. (2010); Benedict et al. (2012, 2017) Consequently, such models do not predict the large sign-asymmetry reported here. In contrast, EPT argues that the essential aspects of many-body screening can be accounted for by treating collisions via the potential of mean force. This has been shown to provide highly accurate transport rates in one-component plasmas, extending from the weakly coupled limit into the strongly coupled regime (for $\Gamma\lesssim 20$).Baalrud and Daligault (2013) However, the theory has not yet been applied to a charge-neutral plasma because a model for the potential of mean force in such a plasma was, until now, not available. The agreement with MD simulations for charge-neutral plasmas validates key concepts of EPT that are neglected in competing transport models. Its accuracy also provides compelling evidence that the Barkas effect in collision rates is a large-angle scattering phenomenon.

The general idea that collision rates can depend on the sign of the charges involved is often referred to as “the Barkas effect”. Barkas, Dyer, and Heckman measured that the stopping power of a medium depends on the sign of the projectile’s charge,Barkas, Dyer, and Heckman (1963); Andersen et al. (1989) an observation at odds with standard models based on Coulomb scattering or linear dielectric response.Bohr (1913); Bethe (1930); Pines and Bohm (1952); de Ferrariis and Arista (1984) The Barkas effect is conventionally understood not as a plasma physics phenomenon but as an atomic physics one, where the projectile either excites or polarizes the bound electrons in solids.Jackson and McCarthy (1972); Lindhard (1976); Ashley, Ritchie, and Brandt (1972) The stopping power of a weakly coupled plasma also exhibits a charge-sign asymmetry when nonlinear polarization of the free electrons is considered, but the effect is small unless the projectile is slow or very highly charged.Peter and Meyer-ter-Vehn (1991) This work presents a new kind of Barkas effect in the macroscopic transport of classical strongly coupled plasmas.

The proposed effect may be observed in UNP experiments, which create a plasma that is simultaneously neutral, classical, and able to achieve moderate to strong coupling in both species.Killian et al. (1999) This makes UNPs the ideal context in which to test the presence of the Barkas effect in transport processes such as ambipolar diffusion into a vacuum,Kulin et al. (2000) electron-ion temperature relaxation,McQuillen et al. (2015) and frictional damping of electron oscillations.Chen, Witte, and Roberts (2017) In particular, the recent measurement of anomalously large electron center-of-mass oscillation damping by Chen et al. at $\Gamma_{e}\sim 0.1-0.35$ (compared with collision models based on repulsive collisions) suggests that the Barkas effect may be responsible.Chen, Witte, and Roberts (2017) Its influence is also expected to extend to electron-ion transport in warm dense plasmas with partially Fermi-degenerate electrons. The existence of a smooth transition from classical to quantum kinetic behavior with increasing electron degeneracy implies that the effect should be present if the electrons are mildly degenerate.Daligault (2017)

This paper is organized as follows. Section II presents the setup, methodology, and results of MD simulations of drift velocity relaxation that show the onset of a Barkas effect. Section III reviews the EPT concept for collisional transport and proposes a semi-analytic model for the effective potential appropriate for UNP conditions. The results of this model are shown to accurately reproduce the relaxation rates in the MD simulations, allowing the Barkas effect to be understood in simple physical terms. It is also explained why the effect is absent from competing models. Section IV describes how the Barkas effect is expected to influence other transport processes. Section V offers some concluding remarks and outlook.

II Molecular Dynamics Simulations

To investigate the Barkas effect at UNP conditions, classical MD simulations of drift velocity relaxation were performed using LAMMPS.lam ; Plimpton (1995) Two separate cases were considered: the relaxation of electrons or positrons, each on ions.

The essential challenge to simulating electron-ion collisional processes in UNPs is that they spend the entirety of their lifetimes far from thermodynamic equilibrium, with $T_{e}>T_{i}$. As a consequence, transport coefficients that depend on electron-ion collisions cannot be studied using the techinques of equilibrium MD, i.e., Green-Kubo relations. Instead, a non-equilibrium methodology must be adopted, but this introduces the difficulty of deciding which transport process to investigate. Temperature relaxation seems like a natural choice, but it is too slow. On the long timescale of temperature relaxation, UNPs are subject to other heating and cooling mechanisms (especially three-body recombination heating) that obscure the electron-ion collision physics.McQuillen et al. (2015) Velocity relaxation, however, is faster than temperature relaxation by a factor of about $m_{i}/m_{e}$. This makes it a suitable probe of electron-ion collisional transport (and thus the Barkas effect), since a drift velocity can relax on a timescale that is short compared to both temperature relaxation and three-body recombination.

All simulations described here followed the same basic methodology: a uniform plasma with unequal temperatures was imparted with a small, uniform drift velocity in the electrons or positrons, which subsequently relaxed due to electron-ion or positron-ion collisions. Section II.1 decribes the simulation setup: the initial and boundary conditions, the interaction potentials, and the choice of time step. Special attention is given to the question of accurately treating the attractive electron-ion interaction.

In order to connect with the theoretical predictions in Sec. III, the simulations were conducted such that the relaxation rate could be identified by modeling the drift velocity with a memoryless linear friction law,

[TABLE]

assuming a constant relaxation rate $\nu_{ei}$ obtained from an exponential fit to the component of $\boldsymbol{V}_{e}(t)$ parallel to the perturbation. This model is justified as long as three criteria are met. First, the drift kinetic energy had to be negligible compared to the thermal kinetic energy, $\frac{1}{2}m_{e}V_{e}^{2}\ll k_{B}T_{e}$. Second, the transient nonlinear response of the electrons or positrons to the velocity perturbation had to be excluded from the analysis. Third, the electrons or positrons had to be in a quasi-steady state, meaning that their temperature and radial distribution functions varied little during the relaxation process. Section II.2 details how the simulations were choreographed to meet these conditions and allow for the identification of a constant relaxation rate suitable for comparing with theory.

The main results of the simulations are presented in Fig. 1, which shows that there exists a Barkas effect in the velocity relaxation rate. The fitting procedure for obtaining the MD relaxation rate as well as the meaning of the error bars are described in Appendix A. The theoretical predictions (“EPT”) are those of the model described in Sec. III. The influence of the effect is evident in two consistent trends in the simulation results: the velocity relaxation rate for an electron-ion plasma is faster than that for a positron-ion plasma, and the difference between the two increases the more strongly coupled the electrons/positrons are. The onset of the asymmetry is only weakly dependent on the ion coupling strength, but it depends strongly on the electron coupling strength. Indeed, the simulation results show that the electrons need only be mildly non-ideal ($\Gamma_{e}\gtrsim 0.2$) in order for the Barkas effect to influence macroscopic transport. This degree of electron coupling has already been surpassed in the UNP experiments on electron center-of-mass oscillation experiments by Chen et al.Chen, Witte, and Roberts (2016, 2017) Thus, it is reasonable to expect that electron-ion transport in present-day UNP experiments may be influenced by the Barkas effect.

The remainder of this section discusses the simulation methodology in detail. An extensive discussion of the motivation, background, and techinical aspects of this work may also be found in Ref. Shaffer, 2018.

II.1 Simulation Setup

All simulations used 5,000 particles of each species. The ion mass and charge were set to $m_{i}=1836m_{e}$ and $Z=1$, respectively, i.e., those of a proton. The particles were initially distributed uniformly throughout a cubic domain with number densities $n_{e}=n_{i}=10^{16}\mathrm{m^{-3}}$. Periodic boundary conditions were employed to emulate an infinite, homogeneous plasma. Initial velocities were drawn from distinct Maxwell-Boltzmann distributions for each species so that initially $T_{e}\geq T_{i}$. Since the number density was constant across simulations and throughout their duration, the coupling strength of each species, $s$, was set by the instantaneous value of the kinetic temperature

[TABLE]

where $\boldsymbol{v}_{s}^{(i)}$ is the velocity of an individual particle of species $s$ and

[TABLE]

is the species drift velocity.

The force on each particle was computed by the particle-particle/particle-mesh method using a cutoff radius of $10a_{e}$ and a $50\times 50\times 50$ $\boldsymbol{k}$-space mesh.Hockney and Eastwood (1981); Plimpton, Pollock, and Stevens (1997) The pairwise interactions for like charges used the Coulomb potential, while the electron-ion pair interactions used the Coulomb potential supplemented by a Gaussian-shaped repulsive core:

[TABLE]

The repulsive core width, $\alpha$, used in the electron-ion simulations is a purely numerical device to prevent rare encounters between electrons and ions that pass too close to resolve with a fixed time step, leading to poor energy conservation.Kuzmin and O’Neil (2002); Hockney and Eastwood (1981) A previous study on the thermodynamics of a plasma interacting with the above potential showed that as $\alpha$ decreased, the non-Coulombic features of the electron-ion interaction could be made sufficiently short-ranged as to not affect the macroscopic properties of the plasma.Tiwari, Shaffer, and Baalrud (2017) However, since the basic physics of the present velocity relaxation problem are different (non-equilibrium versus equilibrium), it is still necessary to establish convergence of the relevant macroscopic observables with respect to $\alpha$.

Figure 2 shows convergence tests with initially moderately coupled electrons, $\Gamma_{e}^{0}=1$. Pictured are the instantaneous electron drift velocity and coupling strength as a function of time after the electron drift was induced. By successively halving $\alpha$, it was found that it could be made sufficiently small not to impact the electron temperature evolution or drift velocity relaxation. It was found that simulations with initially strongly coupled electrons ($\Gamma_{e}^{0}\geq 1$) were most sensitive to $\alpha$, but that a repulsive core width of $\alpha=a/80$ was sufficiently small to reach convergence even in these cases. All simulations of velocity relaxation used $\alpha=a/80$.

The potential well formed by the repulsive core allowed some electrons to fall into Rydberg-atom-like orbits around ions, but this recombination process occurred slowly enough and the lifetime of bound states was typically short enough that only a small fraction of electrons were bound at any given time. Even so, the orbital motion of these few bound electrons set the shortest dynamical timescale in the electron-ion simulations. A rough estimate for this timescale is the orbital period of an electron in circular motion about an ion,

[TABLE]

where the centripetal acceleration is based on the Coulomb force at an orbital radius $\alpha$. For the repulsive core width $\alpha=a/80$ used in the simulations, the characteristic orbital timescale evaluates to $\tau\approx 0.011\omega_{pe}^{-1}$. Accordingly, all electron-ion simulations were performed with a time step $\delta t=0.00025\omega_{pe}^{-1}$, corresponding to $\delta t\approx\tau/43$. Figure 3 shows the variation in total energy of typical electron-ion simulations using this time step, one with initially weakly coupled electrons and one with initially moderately coupled electrons. There was no appreciable drift in the total energy (which would indicate an overall failure to resolve orbital motion), but sharp fluctuations did occur as a result of occasional close electron-ion interactions. The frequency and amplitude of these fluctuations increased with the electron coupling, since less energetic electrons were more likely to form deep and/or long-lived bound states. Bound states did not form in the positron-ion simulations, so a larger time step could be used: $\delta t=0.0025\omega_{pe}^{-1}$ for $\Gamma_{e}^{0}\geq 0.1$ and $\delta t=0.00125\omega_{pe}^{-1}$ for $\Gamma_{e}^{0}<0.1$.

II.2 Non-Equilibrium Methodology

In order to create conditions conducive to identifying a constant relaxation rate, it was necessary for the simulated plasma to reach a quasi-steady state before inducing a drift velocity. Physically, this is the condition that the two species are separately in thermal equilibrium with themselves (but not with one another) and that there is minimal exchange of energy between them. A direct way to test that this is the case is to inspect the evolution of the instantaneous temperature of each species as well as the radial distribution functions (RDFs). Checking the temperatures diagnoses if there is kinetic energy exchange between species, while the RDFs diagnose if there is potential energy exchange.Hansen and MacDonald (1976)

Figure 4 shows the evolution of the drift velocity and instantaneous coupling strength in two sets of simulations, one of electrons and ions and one of positrons and ions. In the electron-ion relaxation, it is seen that the ions slowly heated throughout. This was because the run time of the electron-ion simulations had to be kept shorter than the ion disorder-induced heating (DIH) timescale of $\omega_{pi}^{-1}$ in order to avoid the influence of three-body recombination. However, the $15\omega_{pe}^{-1}$ run time was long enough for electron DIH to saturate to nearly constant $\Gamma_{e}$. In the positron-ion relaxation, the absence of three-body recombination meant that the simulation could have a much longer run time. Both species were given ample time ($60\omega_{pe}^{-1}$) to finish DIH and reach nearly constant temperatures before the positrons were kicked. With these long run times, one can see the slight influence of positron-ion temperature relaxation in that $\Gamma_{i}$ slightly decreased and $\Gamma_{e}$ slightly increased in the late stages of the simulation. However, the positron drift velocity decayed away well before there was any appreciable temperature relaxation.

The need to kick the positrons at a much later time than the electrons was informed by two observations. The first was oscillations in the positron coupling strength at early time shown in Fig. 4, a phenomenon which is analagous to the ion kinetic energy oscillations measured in UNPs.McQuillen et al. (2015) The positron kinetic energy oscillations typically required $10$ to $20\omega_{pe}^{-1}$ to fully decay, depending on their initial coupling strength. The second and more significant reason for delaying the positron kick was the observation that positrons take a much longer time to form screening clouds around ions than do electrons. This was observed in the evolution of the RDFs, pictured in Figure 5, which contrasts the evolution of the RDFs in typical electron-ion and positron-ion simulations. In both cases, $g_{ee}(r)$ quickly reached an essentially static value, signifying that the electrons reached equilibrium with themselves. Also seen in both cases is that $g_{ii}(r)$ relaxed slowly, since the ions did not yet finish disorder-induced heating to their potential energy minimum on the time scale of the figure. Where positrons and electrons differed was in $g_{ei}(r)$, which represents the formation of screening clouds around the ions; $g_{ei}(r)$ becomes static when there is no longer appreciable potential energy exchange between the species. The electron-ion RDF rapidly saturated after DIH, on the order of a few $\omega_{pe}^{-1}$, meaning that quasi-static screening clouds around ions were established. However, the positron-ion RDF evolved much more slowly, requiring a time on the order of the ion plasma period ($\omega_{pi}\approx 43\omega_{pe}^{-1}$) to stabilize. This was in spite of the fact that positron DIH and kinetic energy oscillations subsided much earlier (recall Fig. 4, which is for the same data set shown here).

The fact that positrons formed screening clouds much more slowly than electrons has a simple qualitative explanation. At the start of a simulation, all particles’ positions are distributed uniformly randomly. Electrons rush toward their nearest ion, while positrons rush away from theirs. Thus, after a short time ($\sim\omega_{pe}^{-1}$), the electrons are already screening the charge of their original nearest ions (which have remained essentially stationary), but the positrons have been redirected to some other ion’s vicinity. The positrons continue to be jostled around by ions until they have settled into the interstitial spaces between ions. This requires a time on the order of $\omega_{pi}^{-1}$ because where there are clusters of a few ions (as will be frequently be the case for an initially random configuration), it is unlikely that a positron will be able to occupy the space between them until enough time has passed for the ions to repel one another. Only once they have moved apart and their potential barrier lowered can a positron enter the space.

After the initial waiting period ($5\omega_{pe}^{-1}$ for electrons, $60\omega_{pe}^{-1}$ for positrons) the drift velocity was introduced. An impulse was delivered to the electrons or positrons over a single time step in order to induce a drift velocity small in magnitude compared to the instantaneous thermal speed, $V_{e}(t_{\mathrm{kick}})=0.1v_{Te}(t_{\mathrm{kick}})$, where $v_{Ts}(t)=\sqrt{2k_{B}T_{s}(t)/m_{s}}$. The plasma was then evolved free from external perturbation to allow the drift velocity to relax. For the electrons, simulations were carried out an additional $10\omega_{pe}^{-1}$, which was long enough for their drift velocity to appreciably relax but not so long that a significant population of bound electron-ion pairs could form. Otherwise, inelastic electron-“Rydberg” collisions would contribute an unknown but potentially significant amount to the total damping rate, which would impede comparison with theory. The positrons were allowed to relax for $100\omega_{pe}^{-1}$ after the kick, since recombination was not a concern.

For later comparison with theoretical models, it was important to ensure the velocity distribution functions, $f_{s}(\boldsymbol{v},t)$, did not significantly deviate from a flow-shifted Maxwellian distribution throughout the relaxation process. A simple test of this assumption was done by inspecting the time-evolution of the normalized skewness and excess kurtosis of the reduced distributions, $\bar{f}(v_{z},t)=\iint f_{s}(\boldsymbol{v},t)dv_{x}dv_{y}$. The only non-Maxwellian feature observed was a slight positive excess kurtosis that would develop during the DIH of either species. As a stringent cross-check, the kinetic pressure tensor $\boldsymbol{\mathsf{p}}_{s}(t)=\left<m_{s}[\boldsymbol{v}-\boldsymbol{V}_{s}(t)][\boldsymbol{v}-\boldsymbol{V}_{s}(t)]\right>$ and heat flux $\boldsymbol{q}_{s}(t)=\left<\frac{1}{2}m_{s}[\boldsymbol{v}-\boldsymbol{V}_{s}(t)]|\boldsymbol{v}-\boldsymbol{V}_{s}(t)|^{2}\right>$ were computed from the full distribution functions of an electron-ion simulation with $\Gamma_{i}^{0}=1$ and $\Gamma_{e}^{0}=0.1$. From these, it was confirmed that throughout the simulation, equipartition of kinetic energy was satisfied ($|p_{ii,s}/n_{s}k_{B}T_{s}|<1.05$), viscous stresses were negligible ($|p_{i\neq j,s}/n_{s}k_{B}T_{s}|<0.02$), and heat flux was negligible ($|q_{i,s}/n_{s}k_{B}T_{s}v_{Ts}|<0.05$). Thus, to a good approximation, the MD relaxation process could be modeled assuming flow-shifted Maxwellian distributions throughout.

III Theory

III.1 Effective Potential Theory

EPT is a plasma kinetic theory based on a Boltzmann-like collision operator in which transport occurs through binary interactions between particles, but where many-body effects are included by modeling those interactions as occurring via the potential of mean force.Baalrud and Daligault (2013) In a collision between particles of species $s$ and $s^{\prime}$ with impact parameter $b$ and incident relative speed $u=|\boldsymbol{v}_{s}-\boldsymbol{v}_{s^{\prime}}|$, the angle of deflection is

[TABLE]

where

[TABLE]

is the dimensionless scattering potential containing both centrifugal repulsion and the inter-particle potential. Above, $r_{\min}$ is the distance of closest approach (largest solution to $\Phi_{ss^{\prime}}=1$ at fixed $b$ and $u$) and $m_{ss^{\prime}}=m_{s}m_{s^{\prime}}/(m_{s}+m_{s^{\prime}})$ is the reduced mass.Landau and Lifshitz (1976) Transport coefficients may be expressed in terms of the functions

[TABLE]

which are generalizations of the “Coulomb logarithm” appearing in weakly coupled theory.Baalrud and Daligault (2013) In Eq. (8),

[TABLE]

is the order-$l$ momentum-transfer cross-section,Ferziger and Kaper (1972) $\xi=u/\bar{v}_{ss^{\prime}}$ is a dimensionless relative speed, $\bar{v}_{ss^{\prime}}=\sqrt{v_{Ts}^{2}+v_{Ts^{\prime}}^{2}}$ is the thermal speed associated with the distribution of relative velocities, and $\sigma_{0,ss^{\prime}}=\pi(q_{s}q_{s^{\prime}}/m_{ss^{\prime}}\bar{v}_{ss^{\prime}}^{2})^{2}$ is a reference cross-section.

The premise of EPT is that the appropriate choice of the interaction potential is not the Coulomb interaction, but rather the potential of mean force, $w_{ss^{\prime}}(r)$.Baalrud and Daligault (2013) The potential of mean force is defined as the potential one obtains by holding two particles at fixed positions and ensemble-averaging over the positions of the other particles. Usually, it is defined only for statistical ensembles with a single temperature, in which case it is related to the RDF by $w_{ss^{\prime}}(r)=-k_{B}T\ln g_{ss^{\prime}}(r)$.Hill (1986) For two-temperature plasmas with weak electron-ion coupling, the Boercker-More ensemble may be used to define the potentials of mean force, e.g., $w_{ei}(r)=-k_{B}T_{e}\ln g_{ei}(r)$.Boercker and Pollock (1987)

In the limit of weak coupling, the potential of mean force is the Debye-Hückel potential, $w_{ei}^{\textsc{dh}}(r)=\pm\frac{Ze^{2}}{r}e^{-r/\lambda_{D}}$.Hansen and MacDonald (1976) Binary scattering through a DH potential does result in a momentum-transfer cross section that depends on whether the collision is attractive ($-$) or repulsive ($+$).Khrapak et al. (2003); Baalrud (2012) However, this asymmetry is present only in low-velocity, large-angle collisions such that $\frac{1}{2}m_{ei}u^{2}\ll Ze^{2}/\lambda_{D}$, which are infrequent in a weakly coupled plasma. The resulting order-unity corrections to the Coulomb logarithm are negligible at weak coupling.Liboff (1959); Baalrud (2012) One may expect a larger effect at increased coupling, but the DH potential is no longer valid in this regime, necessitating a more accurate model for the potential of mean force.

III.2 Model for the Potential of Mean Force

For plasmas in thermal equilibrium (i.e., those characterized by a single temperature $T$), an accurate potential of mean force may be obtained from the RDFs that solve the Ornstein-Zernike (OZ) relations

[TABLE]

paired with the hypernetted chain (HNC) closure

[TABLE]

where $h_{ss^{\prime}}(r)=g_{ss^{\prime}}(r)-1$ is the total correlation function, $c_{ss^{\prime}}(r)$ is the direct correlation function, hats denote Fourier transforms, and the sum in Eq. (10) runs over all species labels. There is no exact extension of the theory to two-temperature plasmas, but the proposal by Seuferling, Vogel, and Toepffer (SVT)

[TABLE]

was recently shown to yield RDFs and static structure factors that compare favorably with MD simulations of two-temperature positron-ion plasma.Seuferling, Vogel, and Toepffer (1989); Shaffer, Tiwari, and Baalrud (2017)

The derivation of the SVT equations relies on an assumption that the $N$-body phase-space distribution may be taken in its static limit, that its velocity and position dependence may be decoupled, and that the velocity dependence may be further decoupled into a product of Maxwellians for each species with appropriate temperatures for each. It is also assumed that any relative drift between the species is negligible. These assumptions are consistent with the quasi-steady state constructed in the MD simulations.

It is difficult to compute the RDFs of an electron-ion plasma from either the SVT or OZ equations. The reason is uncontrolled numerical errors that arise from using the attractive Coulomb interaction in the HNC closure. Since the $r\to 0$ behavior of the RDFs is controlled by the Coulomb interaction, $g_{ei}(r)\sim\exp{(1/r)}$ for attractive interactions. By rearranging the HNC closure as $c_{ss^{\prime}}=-\phi_{ss^{\prime}}/k_{B}T_{ss^{\prime}}+g_{ss^{\prime}}-\ln g_{ss^{\prime}}-1$, it is evident that $c_{ei}(r)$ also diverges exponentially as $r\to 0$. On the one hand this behavior seems physically reasonable, since the potential of mean force should recover the bare interaction at close separation. On the other hand, since $g_{ei}(r)$ and $c_{ei}(r)$ each have a non-integrable singularity at $r=0$, they do not have Fourier transforms. Consequently, the SVT or OZ equations cannot be solved in the usual way.Springer, Pokrant, and Stevens (1973); Ng (1974)

An attempt was made to circumvent these numerical difficulties by artificially softening the electron-ion interaction

[TABLE]

which takes the value $\phi_{ei}(0)=-Ze^{2}/\alpha$ at zero separation. This introduces a new parametric dependence on the softening length scale $\alpha$, and it is the $\alpha\to 0$ limit which is relevant to UNPs. To reach the lowest numerically feasible value of $\alpha$, a sequence of HNC calculations at fixed $\Gamma_{i}$ and $\Gamma_{e}$ were performed with successively decreasing $\alpha$. It was found that the iterative solution to the HNC-SVT equations could only be made to converge when $\alpha/a_{e}\gtrsim\Gamma_{e}/10$. At the moderate electron coupling strengths where the Barkas effect was observed in MD, this was not small enough to extrapolate to an accurate $\alpha\to 0$ limit. Details of this attempt and its failure to describe the MD simulation results are given in Appendix B.

A more fruitful approach is to exploit the fact that the electron coupling strength is not very high in UNPs, which enables a connection between the SVT equations and the screened one-component plasma model.Shaffer, Tiwari, and Baalrud (2017) The SVT equations in the limit $m_{e}/m_{i}\to 0$ can be rearranged in terms of the partial static structure factors $S_{ss^{\prime}}(k)=\delta_{ss^{\prime}}+\sqrt{n_{s}n_{s^{\prime}}}\hat{h}_{ss^{\prime}}(k)$ as

[TABLE]

Taking $n_{e}=Zn_{i}$ and treating the electrons (or positrons) in the Debye-Hückel approximation amounts to setting

[TABLE]

where plus or minus signs refer to electrons and positrons, respectively. One then finds that the electron-ion structure factor simplifies to

[TABLE]

and the ion-ion structure factor becomes that of a one-component plasma with coupling strength $\Gamma=\Gamma_{i}$ and screening parameter $\kappa=a_{i}/\lambda_{De}=Z^{1/3}\sqrt{3\Gamma_{e}}$.Shaffer, Tiwari, and Baalrud (2017) Note that in this treatment, only the electrons (or positrons) are treated in the Debye-Hückel approximation; no assumption on the ion coupling strength is made.

In the present approximation, the structure factors are in fact the same as those derived by Boercker and More from their two-temperature statistical ensemble.Boercker and Pollock (1987) The HNC closure, Eq. (12a), may then used to compute an approximate electron-ion potential of mean force, $w_{ei}(r)=\phi_{ei}(r)-k_{B}T_{e}[h_{ei}(r)-c_{ei}(r)]$ which is consistent with its statistical definition in terms of the Boercker-More ensemble. In this approximation,

[TABLE]

meaning the electron-ion potential of mean force may be computed with only the knowledge of the structure factor of an appropriately parameterized screened OCP.

The approximations made in Eq. (15) are appropriate when the electrons are not strongly coupled. This is the case in all UNP experiments done to date, where $\Gamma_{e}\leq 0.35$,Chen, Witte, and Roberts (2017) as well as in the majority of the MD simulations described in Sec. II. The accuracy of Eq. (18) for the potential of mean force has been verified both by comparison with fully two-component SVT solutions for positron-ion plasmas,Shaffer, Tiwari, and Baalrud (2017) as well as by comparison with MD simulations, shown in Fig. 6.

An advantageous side-effect of using this approximation is that Eq. (18) differs between electrons and positrons only in overall sign of the potential. In other words, even though positrons and electrons have different charge, in the limit that they are not too strongly coupled, they screen the ions identically. This means that whatever Barkas effect is exhibited in the EPT results that follow from this approximation arise solely due to this leading sign of the effective potential. This is an important point for identifying the physical mechanism responsible for the observed Barkas effect.

III.3 Origin of the Barkas Effect

The solid curves in Fig. 1 show the result of evaluating EPT for the velocity relaxation rate at the conditions of the MD simulations, $\nu_{ei}=\sqrt{2/3\pi}\Gamma_{e}^{3/2}\Xi_{ei}^{(1,1)}\omega_{pe}$, derived from the momentum moment of the collision operator (see Sec. IV.1 for more details). Since the final value of the ion coupling strength in the MD simulations depends slightly on $\Gamma_{e}$ (and therefore not all points in the same frame have the same $\Gamma_{i}$), the EPT relaxation rates are shown as a band of values spanning the highest and lowest values that occurred in the simulations. Across all coupling strengths investigated, the EPT predictions qualitatively reproduce the Barkas effect seen in the MD simulations. Furthermore, for $\Gamma_{e}\lesssim 0.5$, EPT accurately predicts the numerical value of the relaxation rate. At higher electron coupling strength, the approximation of Sec. III.2 used in the EPT calculations is no longer valid, so it is unsurprising that the model for the potential of mean force used here does not reproduce the MD relaxation rates in that regime.

The accuracy of EPT in predicting the Barkas effect sheds light on its physical origin, namely that it can be understood in terms of the kinematics of attractive versus repulsive screened collisions. The central dynamical quantity for understanding the Barkas effect is the scattering potential, $\Phi_{ei}(r,b,u)=b^{2}/r^{2}+2w_{ei}(r)/m_{ei}u^{2}$. Figure 7 shows how the sign of $w_{ei}(r)$ alters the structure of the scattering potential and how this carries through to the deflection angle and cross-section. For high-velocity collisions ($\xi\geq 1$), the collision dynamics are controlled mainly by the centrifugal repulsion $b^{2}/r^{2}$, which is sign-independent. The Barkas effect instead is prominent at lower relative velocities ($\xi<1$). For positron-ion collisions, decreasing velocity causes $\Phi_{ei}(r)$ to smoothly transition from being centrifugally dominated to being dominated by $w_{ei}(r)$. The cross-section likewise transitions from being dominated by small-angle scattering to large-angle scattering. Electron-ion collisions, however, allow the possibility of spiral scattering when $\xi^{2}\lesssim Z\Gamma_{e}^{3/2}$ whereby the centrifugal and attractive terms of $\Phi_{ei}(r)$ compete at intermediate range to form a potential well.Eisberg and Porter (1961); Khrapak et al. (2003) This behavior is illustrated in the solid $\xi=0.2$ curve of Fig. 7a. As an electron traverses this well, it experiences a large angular deflection, sometimes involving multiple full rotations of the electron around the ion before eventually scattering away, as seen in the $\xi=0.2$ curve of Fig. 7b. The net effect is that electrons scatter through large angles over a wider range of impact parameters than do positrons, leading to a bigger electron-ion cross-section for momentum transfer compared to positron-ion collisions. In turn, this leads to enhanced transport rates.

III.4 Absence from Other Models

Besides EPT, there are two main theoretical approaches to collisional transport in strongly coupled plasma: those that extend the Landau-Spitzer theoryLee and More (1984); Li and Petrasso (1993); Dharma-wardana and Perrot (1998); Gericke, Murillo, and Schlanges (2002); White et al. (2017) and those that extend the Lenard-Balescu theory.Daligault and Dimonte (2009); Vorberger et al. (2010); Benedict et al. (2012, 2017) Neither of these approaches predict the effect shown here. The reason is that neither accurately models how screening influences close interactions.

Extensions of the Landau-Spitzer theory retain the idea that the binary Coulomb collision is the basic unit of transport. That is, the deflection angle is given by the Rutherford formula $\theta_{ss^{\prime}}=2\arctan(q_{s}q_{s^{\prime}}/m_{ss^{\prime}}u^{2}b)$, so that the (first) momentum-transfer cross-section takes the form $\sigma^{(1)}_{ss^{\prime}}(u)=2\pi(q_{s}q_{s^{\prime}}/m_{ss^{\prime}}u^{2})^{2}\ln[1+(m_{ss^{\prime}}u^{2}b_{\max}/q_{s}q_{s^{\prime}})^{2}]$ after truncating the divergent integral in Eq. (9) at a characteristic interaction range $b_{\max}$. It is in the choice of $b_{\max}$ that these approaches attempt to insert many-body screening physics, typically by interpolating between $\lambda_{D}$ for weak coupling and $a_{i}$ for strong coupling.Lee and More (1984); Gericke, Murillo, and Schlanges (2002) However, the resulting cross-section is independent of whether the collision is attractive or repulsive. This is a peculiarity of inverse-square forces; the assumption of Coulomb collisions prevents Landau-Spitzer-based models from exhibiting a Barkas effect.

The lack of a Barkas effect in Landau-Spizter-based models is a symptom of not accounting for screening in the kinematics of close collisions. The effect cannot be adequately recovered by means of altering $b_{\max}$. Recall that the EPT model of Sec. III.2 was able to reproduce the Barkas effect observed in MD by changing only the sign — not the range — of the interaction potential. While one could imagine prescribing a $b_{\max}$ that reproduces the Barkas effect, the underlying collision physics would be dubious. This shortcoming of the Landau-Spitzer approach is evocative of the Salpeter enhancement of nuclear reaction rates, where screening between ions decreases the energy barrier for close collisions, a phenomenon also missed in a treatment based on impact parameter cutoffs.Salpeter (1954); Anderegg et al. (2009)

The Barkas effect is also not accurately treated in linear response theory, such as extensions of the Lenard-Balescu collision operator. In these models, strong-coupling effects are included through static local field corrections (LFCs), $G_{ss^{\prime}}(k)$, to the polarization potential around charges. For instance, the electron-ion relaxation rate derived by Daligault and Dimonte takes the form (for classical plasmas)

[TABLE]

where $\chi_{s}^{0}(k,\omega)$ is the free-particle density-density response function, and the function $\hat{\phi}_{ei}(k)/|D(k,\omega)|$ may be interpreted as a screened electron-ion interaction.Daligault and Dimonte (2009) (See Eqs. (23-26) in Ref. Daligault and Dimonte, 2009 for comparison with related models.) Since the electron-ion interaction only appears squared, a charge-sign asymmetry can only arise via a difference between electron-ion versus positron-ion LFCs. That is, generalized Lenard-Balescu models predict a sign asymmetry due to the shape of particles’ dielectric dressing, but not the overall sign of the interaction. Physically, this is because approaches based on linear dielectric response model interactions via the correlation of linear fluctuations. This excludes large-angle close collisions, which are a nonlinear phenomenon in this respect.

This basic limitation of linear response-based collision models becomes especially clear in comparison to the EPT model derived in Sec. III.2. The point of connection is the relationship between the direct correlation functions and the static LFCs,Ichimaru (1992); Hansen and MacDonald (1976)

[TABLE]

The assumption of weak electron coupling in the present model (Sec. III.2) translates directly to setting $G_{ei}(k)=G_{ee}(k)=0$. Thus, to the same level of approximation used in the EPT computations shown in Fig. 1, the linear polarization potential approach does not predict a Barkas effect. Such an approximation is necessary for making quantitative predictions because it is not known how to compute the LFCs of a classical electron-ion plasma for reasons described in Ref. Daligault and Dimonte, 2009 and in Appendix B. Even if one had a good model for $G_{ei}(k)$, it seems unlikely that using it in generalized Lenard-Balescu theory will accurately capture the Barkas effect, since it arises primarily due to the sign of the screened potential, not the shape.

IV Impact on Transport Processes

IV.1 Velocity Relaxation in UNPs

In a uniform plasma, the electron drift velocity evolves according to $m_{e}n_{e}\partial_{t}\boldsymbol{V}_{e}=\boldsymbol{R}_{ei}$, where $\boldsymbol{R}_{ei}$ is the rate of change of the electron momentum density due to electron-ion collisions, i.e., the momentum moment of the collision operator. For the Boltzmann collision operator, this may be evaluated asFerziger and Kaper (1972)

[TABLE]

where $\boldsymbol{u}=\boldsymbol{v}_{e}-\boldsymbol{v}_{i}$. Eq. (21) is general for any form of the distribution functions, but when they are flow-shifted Maxwellians, $\boldsymbol{R}_{ei}$ takes the simple form

[TABLE]

where the velocity-dependent relaxation rate is

[TABLE]

in which $\eta=|\boldsymbol{V}_{e}-\boldsymbol{V}_{i}|/\bar{v}_{ei}$ is the dimensionless relative drift speed, $\bar{\Xi}_{ei}$ is a velocity-dependent generalization of $\Xi_{ei}^{(1,1)}$,

[TABLE]

and $H(\xi,\eta)=(2\xi\eta+1)e^{-(\xi+\eta)^{2}}+(2\xi\eta-1)e^{-(\xi-\eta)^{2}}$.Baalrud and Daligault (2016) When the drift kinetic energy is small compared to the thermal energy ($\eta^{2}\ll 1$),

[TABLE]

where $c=[6\Xi_{ei}^{(1,2)}-15\Xi_{ei}^{(1,1)}]/5\Xi_{ei}^{(1,1)}$. Since ratios of generalized Coulomb logarithms are order-unity quantitiesBaalrud and Daligault (2014) and the simulations in Sec. II were performed at conditions such that $\eta\leq 0.1$, the drift-velocity dependent corrections to the generalized Coulomb logarithm may be neglected here. With the further simplifications that $m_{e}\ll m_{i}$, $T_{e}\geq T_{i}$, and $Zn_{i}=n_{e}$, a theoretical prediction for the velocity relaxation rate pertinent to the present MD simulations may be computed from

[TABLE]

The EPT data shown in Fig. 1 is the result of evaluating Eq. (26) using the potential of mean force model described in Sec. III.2.

Velocity relaxation was the basic process measured in the experiments by Chen et al. on electron center-of-mass oscillation damping in UNPs.Chen, Witte, and Roberts (2016, 2017) Comparison of the measured damping rates with those obtained by hybrid MD-Monte Carlo simulations implementing various models for the electron-ion collision operator showed that the theories systematically underestimated the damping rate. However, all the theories considered make either a weak-coupling approximation that precluded any Barkas effect or were based on models for strongly coupled plasmas with repulsive interactions. Based on Fig. 8, though, at the strongly coupled conditions reported in those experiments, $\Gamma_{e}=0.35$, EPT predicts that the Barkas effect could contribute as much as a $50-70\%$ increase to the predicted collision rate compared to theories based on repulsive interactions, depending on the ion coupling strength. If binary collisions were the dominant mode of damping in the experiments, then correcting for the sign of the electron-ion interaction could close the gap between collision models and measurements reported in Ref. Chen, Witte, and Roberts, 2017. However, at the high $\Gamma_{e}$ they obtained, inelastic collisions between electrons and Rydberg atoms are also expected to influence the damping of the electron motion,Roberts ; Gupta et al. (2007) and it is not clear at present which is the dominant damping mechanism as $\Gamma_{e}$ increases.

IV.2 Temperature Relaxation in UNPs

The temperature evolution in a uniform plasma is given by $\frac{3}{2}m_{e}n_{e}k_{B}\partial_{t}T_{e}=Q_{ei}$, where $Q_{ei}$ is the rate of change of the electron internal energy density due to collisions, i.e., the $\frac{1}{2}m_{e}|\boldsymbol{v}_{e}-\boldsymbol{V}_{e}|^{2}$ moment of the collison operator. For the Boltzmann collision operator,Ferziger and Kaper (1972)

[TABLE]

Again, taking the distribution functions to be flow-shifted Maxwellians, one findsBaalrud and Daligault (2016)

[TABLE]

The first term describes temperature relaxation between the electrons and ions at the velocity-dependent rate

[TABLE]

which involves a different velocity-dependent Coulomb logarithm

[TABLE]

with $I(\xi,\eta)=e^{-(\xi-\eta)^{2}}-e^{-(\xi+\eta)^{2}}$.Baalrud and Daligault (2016) For small relative drift speeds,

[TABLE]

where $c=[2\Xi_{ei}^{(1,2)}-3\Xi_{ei}^{(1,1)}]/3\Xi_{ei}^{(1,1)}$ is again of order unity. The second term in Eq. (28) represents frictional heating; since it is $O(\eta^{2})$ compared to the temperature relaxation, it is disregarded here. Comparing Eq. (29) with Eq. (23) for small drift velocity, it is evident that both the temperature and velocity relaxation depend on the same quantity $\bar{\nu}_{ei}(0)=\tilde{\nu}_{ei}(0)=\nu_{ei}$, given by Eq. (26). Thus the Barkas effect influences the temperature relaxation rate in exactly the same way as shown in Fig. 8 for the velocity relaxation rate.

IV.3 Electrical Conductivity in Dense, Thermal Plasmas

The Barkas effect is also expected to extend to electron-ion transport in dense plasmas with partially Fermi-degenerate electrons. Since EPT has not yet been formulated for quantum-mechanical transport, it would be incorrect to directly apply the results obtained so far for classical plasmas to dense plasmas. Nevertheless, it is expected that the transport coefficients of dense plasmas are smooth functions of the electron degeneracy.Daligault (2017, 2018) This implies that a Barkas effect should still be present in these systems at partial electron degeneracy.

In dense plasmas, theoretical models for electron transport are difficult to validate due to a lack of a first principles quantum simulation method. An alternative approach is to treat the dynamics of the plasma classically but to model quantum-mechanical effects by altering the interaction potential between particles. The main reason to do so is that it allows the use of classical MD to test theory. A recent example is the semi-classical model of electrical conductivity presented by Whitley et al.,Whitley et al. (2015) who performed MD simulations of hydrogen plasmas using the Dunn-Broyles potentialDunn and Broyles (1967)

[TABLE]

with $\lambda_{ss^{\prime}}=\hbar/\sqrt{2m_{ss^{\prime}}k_{B}T}$. They compared their MD results with theoretical predictions based on the classical Lenard-Balescu collision operator. Due to the use of the Lenard-Balescu collision operator, there is no charge-sign dependence in their model for the electrical conductivity.

For comparison, EPT predictions for the electrical conductivity were computed from the second-order Chapman-Enskog formula,

[TABLE]

taken in the approximation $m_{e}\ll m_{i}$.Ferziger and Kaper (1972) The factor $(1-\Delta)^{-1}$ is the “Spitzer” correction accounting for skewness in the electron distribution function due the applied electric field.Spitzer and Härm (1953) The electron-electron and electron-ion Coulomb logarithms were computed using the potentials of mean force obtained from equilibrium HNC using the Dunn-Broyles potential above in Eqs. (10) and (11). That is, no weak-coupling approximations were made. This was possible because the Dunn-Broyles potential is very soft as $r\to 0$ for the densities and temperatures considered here, so these HNC calculations did not suffer from the numerical difficulties described in Sec. III.2 or Appendix B.

Figure 9 compares EPT with the Lenard-Balsecu model and MD simulations of Ref. Whitley et al., 2015 for a fully ionized hydrogen plasma with mass density $\rho=40\mathrm{g/cm^{3}}$ and $k_{B}T=500,700,900\mathrm{eV}$. Two sets of EPT results are shown; they differ only in the overall sign of the electron-ion potential of mean force. At these conditions, the plasma is moderately coupled ($\Gamma=e^{2}/ak_{B}T\sim 0.1$) and the electron-ion and electron-electron interactions are quite soft ($\lambda_{ei}\sim 2a_{e}$). Both the Lenard-Balescu model and attractive EPT offer reasonable predictions for the conductivity in light of the large variability in the MD results. However, in comparing the attractive and repulsive EPT calculations, it is seen that the repulsive interaction leads to a conductivity that is systematically about $10\%$ smaller. That is, softening the Coulomb interaction results in attractive interactions having a smaller collision rate compared to repulsive ones. This means that within this simplistic model, not only does a charge-sign asymmetry effect persist, but in fact it inverts compared to the classical case. This is in part because in the Dunn-Broyles interaction the large-angle and spiral scattering events described in Sec. III.3 onset at lower velocities compared to the attractive Coulomb interaction, so they do not contribute significantly to transport at the conditions considered here. It remains to be seen if a more accurate treatment of quantum-mechanical effective-potential scattering would exhibit a similar inversion of the Barkas effect.

V Conclusions

It has been demonstrated that the transport rates of plasmas exhibit a dependence on the sign of the electron charge analogous to the Barkas effect in charged particle stopping. This is an emergent consequence of strong coupling that arises when one makes a detailed account of how screening in the plasma alters the trajectories of binary encounters between electrons and ions; the effect cannot be recovered by treating collisions as Coulomb interactions with truncated range, nor is it captured by the leading theories of classical transport in strongly coupled plasma based on linear response. The effect is small at weak electron coupling where these approaches still approximately hold, but as the electron coupling strength approaches unity, the electron-ion collision rate is significantly enhanced compared to that for positron-ion collisions. UNP experiments may be able to measure this fundamental physical effect in plasma transport. It is also expected to influence essentially all neutral plasmas at strong coupling, including dense plasmas (though electron degeneracy must also be accounted for to accurately treat these systems).

Acknowledgements

We wish to thank Jérôme Daligault and Jacob Roberts for stimulating conversations, as well as Lorin Benedict for sharing the LB and MD electrical conductivity data shown in Fig. 9. This material is based upon work supported by the National Science Foundation under Grant No. PHY-1453736. It used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF Grant No. ACI-1053575 under Project Award No. PHY-150018.

Appendix A Fitting and Error Estimation of MD Results

Because the drift velocity was necessarily small compared to the electron/positron thermal speed, each individual run contained significant fluctuations in $\boldsymbol{V}_{e}(t)$. Thus, fitting individual drift-velocity time series was inappropriate. Instead, the sample mean and standard deviation of the drift velocity were used to obtain the relaxation. The “sample” here means that the mean and standard deviation were taken over 100 independent replicas of a single combination of $\Gamma_{i}^{0}$ and $\Gamma_{e}^{0}$. Specifically, the relaxation rate was determined by minimizing Press et al. (2007)

[TABLE]

with respect to the fit parameters $V_{0}$ and $\nu_{ei}$, where $N_{t}$ is the number of time steps in the fitting interval,

[TABLE]

is the sample mean drift velocity across $N_{r}$ runs, and

[TABLE]

is the sample variance. Confidence intervals for the fit parameters were taken from the diagonals entries of the fit parameter covariance matrix; the off-diagonal entries were always negligible in comparison. The uncertainties in the fit parameters were typically very small (parts per hundred), as a result of the fit being highly over-constrained (thousands of time steps to determine two fit parameters).

The start of the time interval used to fit the drift velocity was chosen to be late enough after the kick that the non-exponential transient response was excluded. Ideally, the fit should be started as late as possible, just not so late that the drift velocity falls below the thermal noise floor. However, in the electron-ion systems, the competing requirement that the total run time be short constrains how late the fitting window can start. It was found that starting $5\omega_{pe}^{-1}$ after the kick was suitable. This was long enough after the kick that non-exponential decay was not apparent, yet still early enough that a substantial time was included ($10\omega_{pe}^{-1}$). For the longer positron-ion simulations, there was much more flexibility in choosing the fitting interval. The start was chosen to be $10\omega_{pe}^{-1}$ after the kick, and the end was chosen to be either the end of the simulation or the time at which the drift velocity had decayed by three $e$-foldings, whichever occurred earlier. The coupling strengths associated with the fitted relaxation rate were the average $\Gamma_{e}$ and $\Gamma_{i}$ over the fitting interval. Conservative estimates for the “error” in the coupling strengths were computed from the minimum and maximum values of $\Gamma_{e}$ and $\Gamma_{i}$ over the whole fit interval and across all included runs.

Appendix B Numerical Details of Softened Electron-Ion HNC Calculations

The HNC calculations presented in this paper were done using Fozzie, a free and open-source program to numerically solve the OZ or SVT equations.Foz The solution method is an iterative scheme in the same vein as Refs. Springer, Pokrant, and Stevens, 1973; Ng, 1974. The problem of solving the OZ or SVT equations subject to the HNC closure is cast as a fixed-point problem for the direct correlation functions, $\mathcal{A}[\boldsymbol{\mathsf{c}}(r)]=\boldsymbol{\mathsf{c}}(r)$, where $\boldsymbol{\mathsf{c}}$ is the two-by-two matrix whose elements are $c_{ss^{\prime}}$, and the operator $\mathcal{A}$ represents a sequence of steps: (1) Fourier transforming each $c_{ss^{\prime}}(r)\to\hat{c}_{ss^{\prime}}(k)$, (2) solving the OZ or SVT relations for the indirect correlation functions, $\hat{\gamma}_{ss^{\prime}}(k)\equiv\hat{h}_{ss^{\prime}}(k)-\hat{c}_{ss^{\prime}}(k)$, (3) inverse Fourier transforming each $\hat{\gamma}_{ss^{\prime}}(k)\to\gamma_{ss^{\prime}}(r)$, and (4) evaluating $\exp(-\phi_{ss^{\prime}}/k_{B}T_{ss^{\prime}}+\gamma_{ss^{\prime}})-\gamma_{ss^{\prime}}-1$. Given a trial solution $\boldsymbol{\mathsf{c}}^{(i)}$, the application of $\mathcal{A}$ produces a new trial solution $\boldsymbol{\mathsf{c}}^{(i+1)}=\mathcal{A}[\boldsymbol{\mathsf{c}}_{i}]$. This trial solution was then linearly mixed with the previous one: $\boldsymbol{\mathsf{c}}^{(i+1)}\leftarrow\zeta\boldsymbol{\mathsf{c}}^{(i+1)}+(1-\zeta)\boldsymbol{\mathsf{c}}^{(i)}$, where the mixing parameter $\zeta$ was set a priori. Self-consistency was judged by computing $\delta^{(i)}=\sqrt{\int_{0}^{\infty}\frac{1}{4}\sum_{s,s^{\prime}=1}^{2}|c_{ss^{\prime}}^{(i+1)}-c_{ss^{\prime}}^{(i)}|^{2}\,dr}$. If $\delta^{(i)}<10^{-8}$, then the calculation was ended and $\boldsymbol{\mathsf{c}}^{(i+1)}$ was taken to be the solution. Otherwise, a new trial solution was generated from $\boldsymbol{\mathsf{c}}^{(i+1)}$.

The rate of convergence depended both on the initial trial solution $\boldsymbol{\mathsf{c}}^{(0)}$ and the mixing parameter $\zeta$. When the system to solve had strong attractive interactions (as in the case $\alpha\to 0$), the iteration did not converge at all unless $\boldsymbol{\mathsf{c}}^{(0)}$ was already close to a solution and the iterative refinement was forced to progress slowly by setting $\zeta\ll 1$. To obtain HNC solutions with the smallest feasible values of $\alpha$, a sequence of solutions were generated, each with a smaller value of $\alpha$ than its predecessor. To aid convergence, each run was initialized with the converged result of the preceding one. If a run failed to converge, it was retried with a value of $\alpha$ closer to the previous run’s and a more conservative choice of the mixing parameter $\zeta$. This continued either until 50 attempts had been made or until the values of $\alpha$ between two successive attempts differed by less than $0.05\%$. It was observed that this routine tended to terminate at $\alpha/a_{e}\approx\Gamma_{e}/10$. The resulting potentials of mean force with those from MD showed unacceptable non-Coulombic behavior when $\Gamma_{e}\sim 0.1-1$, even when $r\sim a(\gg\alpha)$.