Using molecular dynamics simulations, we perform the first direct tests of three proposed models for the pair correlation functions of strongly coupled plasmas with species of unequal temperature. The models are all extensions of the Ornstein–Zernike/hypernetted-chain theory used to good success for equilibrium plasmas. Each theory is evaluated at several coupling strengths, temperature ratios, and mass ratios for a model plasma in which the electrons are positively charged. We show that the model proposed by Seuferling et al. [Phys. Rev. A 40, 323 (1989)] agrees well with molecular dynamics over a wide range of mass and temperature ratios, as well as over a range of coupling strength similar to that of the equilibrium hypernetted-chain (HNC) theory. The SVT model also correctly predicts the strength of interspecies correlations and exhibits physically reasonable long-wavelength limits of the static structure factors. Comparisons of the SVT model with the Yukawa one-component plasma (YOCP) model are used to show that ion-ion pair correlations are well described by the YOCP model up to Γe1, beyond which it rapidly breaks down.

Strongly coupled plasmas produced in experiments are often far from thermal equilibrium. Ultracold neutral plasmas and the dense plasmas in sonoluminescent bubbles, for instance, can have electron and ion temperatures that differ by an order of magnitude or more.1,2 In inertial confinement fusion plasmas and ultracold plasma mixtures, there can be significant differences in temperatures not just between ions and electrons, but also between the different species of ions.3,4 The combination of strong coupling and multiple temperatures makes these plasmas especially challenging to model, since one cannot freely call upon results from equilibrium statistical mechanics to make predictions about the plasma's thermodynamic and transport properties. One approach for strongly coupled, two-temperature plasmas is to extend integral equation theories for the equilibrium pair distribution functions to allow multiple temperatures. In this work, we present the first direct comparisons of three such extensions and evaluate their accuracy against molecular dynamics (MD) simulations.

Pair distribution functions are normally applied in the context of thermal equilibrium, where they can be used to evaluate the pressure, internal energy, and other thermodynamic state variables using exact formulas from equilibrium statistical mechanics. In non-equilibrium plasmas—especially those far from equilibrium—they are used in extensions of the ideal gas kinetic theory to treat strongly coupled plasmas. The pair distributions enter into these models in the form of effective scattering potentials5,6 or local field corrections,7–11 designed to take approximate account of how the collisional transfer of momentum and energy in the plasma is affected by the many-body physics of strong coupling. The pair distribution functions of two-temperature systems are also useful for testing the range of validity of approximate one-component models of strongly coupled plasmas. Most importantly, by treating electrons and ions on equal footing, two-component plasma models grant access to electron-ion transport physics that lie beyond the scope of a one-component treatment.

The present work makes use of a model plasma consisting of ions and positively charged electrons. This approach is useful in both modeling and simulation (e.g., Ref. 9) to circumvent the collapse (recombination) of classical electron-ion plasmas, which is typically treated using softened electron-ion pseudopotentials. By instead using positively charged electrons, we are able to isolate the relevant two-temperature physics, which should not depend on the sign of the charge. Future work will use a recently developed method for modeling strongly coupled electron-ion plasmas12 to explore the effect of negatively charged electrons on pair correlations and transport. Notwithstanding, the results shown here are immediately applicable to ionic mixtures with unequal temperatures.6 

At weak coupling, the pair distribution functions are accurately described by the Debye–Hückel theory of electrolytes.13 For strongly coupled plasmas, however, the triplet and higher-order correlations ignored in the Debye–Hückel approach become important. At thermal equilibrium, these correlations are well approximated by integral equation methods developed from equilibrium statistical mechanics. The most successful of these involve solving the Ornstein–Zernike (OZ) relations together with an approximate closure, e.g., the hypernetted-chain (HNC) approximation.14 When the plasma has more than one temperature, two methodologies have been explored: (a) to map the multi-temperature plasma to an effective one-temperature plasma or (b) to extend equilibrium integral equation theories to allow multiple temperatures.

The canonical example of the mapping approach is the Yukawa one-component plasma (YOCP) model. In the YOCP model, the plasma is partitioned into a strongly coupled component of interest and a weakly coupled background. All the physics of this background are condensed into a constant screening parameter that modifies the interaction between the strongly coupled particles. These conditions are realized in dusty plasmas and in present-day ultracold neutral plasma experiments, where the YOCP model has been successfully applied to study the dust and ions, respectively.15–17 However, the YOCP cannot be used to describe processes that involve electrons, e.g., ambipolar diffusion or electron-ion temperature relaxation.

The other class of approaches extends the theory of equilibrium density correlations to the case of a plasma with two distinct temperatures. In an early investigation, Salpeter derived pair correlation functions for a weakly coupled electron-ion plasma using arguments in the vein of the Debye–Hückel theory.18 Boercker and More extended Salpeter's results to strong ion coupling using an ansatz for a two-temperature partition function, but still required that the electron-ion Coulomb coupling be weak.19 The extension to arbitrary coupling involves a generalization of the theory of pair correlations in equilibrium liquids. The approaches considered here introduce the notion of a “cross” temperature Tab that serves as the kinetic energy scale for inter-species correlations. One must also determine if the OZ equations themselves should be modified. An attractive feature of such models is that all species are treated on equal footing, in contrast to the YOCP. This permits the direct calculation of all pair correlation functions and further allows for the possibility of studying two-temperature physics when both species are strongly coupled.

Our main goal is to determine the most accurate approximation available to extend equilibrium integral-equation theories of density correlations to two-temperature plasmas. There seems to be no consensus at present regarding the form of the cross temperatures, Tab, or whether it is necessary to modify the Ornstein–Zernike relations. This work considers three formulations20–22 that have appeared in recent work on two-temperature strongly coupled plasmas.11,21–24 Our main finding is that the model proposed by Seuferling et al. (“SVT”) in Ref. 20 predicts pair distribution functions that agree with MD over a range of coupling strengths similar to what is seen for the usual equilibrium HNC theory.

We restrict our scope to a plasma with two species of classical point charges, labeled i and e, with distinct masses and temperatures. We focus on testing cases where TeTi and memi, i.e., the lighter “electrons” are warmer than the massive “ions.” This is the parameter regime of greatest importance in current strongly coupled plasma contexts. We also take both species to have equal number density (ne=ni=n/2) and unit charge Zi=Ze=1, so that the interaction potential for all particles is the repulsive Coulomb potential

vab(r)=e2r,
(1)

and the Coulomb coupling strength of each species is

Γs=e2/askBTs,
(2)

where as=(3/4πns)1/3 is the mean spacing between particles of species s, Ts is their temperature, e is the elementary charge, and kB is the Boltzmann constant.

Another basic assumption of this work is the existence of a two-temperature steady state, or “quasi-equilibrium.” In a plasma, the collisional exchange of energy tends to be most efficient between particles of the same mass and least efficient between particles of very different mass. It is frequently the case that particles of each species equilibrate among themselves before the system as a whole relaxes to thermal equilibrium. On timescales longer than the intraspecies thermal relaxation time but shorter than the interspecies thermal relaxation time, it is often accurate to take the velocity distributions to be Maxwellian with temperatures Ti and Te.

Section II introduces the three theories and discusses some of their asymptotic limits. Section III provides details on the MD techniques used to simulate a two-temperature quasi-equilibrium plasma. Section IV compares the pair distribution functions of the theoretical models with the MD results. Section V uses the SVT model to study when the YOCP model for ion-ion correlations breaks down as the electron coupling strength increases. Section VI offers some concluding remarks and describes how the present results will be of use to future studies of two-temperature plasmas.

At thermal equilibrium, the Ornstein–Zernike (OZ) relations are25 

ĥab(k)=ĉab(k)+s=i,ensĥas(k)ĉsb(k),
(3)

where ĥab(k) and ĉab(k) are the Fourier transformed total and direct correlation functions, respectively, and k is the wavenumber. The OZ equations must be solved in conjunction with approximate closure relations. The hypernetted-chain (HNC) closure is given by26 

gab(r)=exp[vab(r)kBT+hab(r)cab(r)],
(4)

where gab(r)=1+hab(r) are the radial distribution functions (RDFs). The HNC closure is very accurate when vab(r) is long-ranged, as is the case for the Coulomb potential. However, it is reasonable to expect that bridge functions (which are neglected in HNC) contribute non-negligibly to the RDFs when max(Γi,Γe)10, based on knowledge of the OCP RDFs. A number of proposed improvements model the neglected bridge functions; see for example, Refs. 27–30. We revisit the approximate nature of the HNC closure when comparing with MD results in Sec. IV.

To extend the theory of spatial correlations at equilibrium to multi-temperature systems one must address two points: (a) how to characterize the “cross-temperatures” Tab that set the kinetic energy scale for inter-species correlations and (b) whether the OZ relations should be modified. Most investigations to date have extended the HNC–OZ system of equations in one of three ways. We will refer to them in this work as the SQRT, MASS, and SVT models.

In the SQRT21 and MASS22 models, the OZ relations are taken to be the same as at equilibrium, and the HNC closures are assumed to be

gab(r)=exp[vab(r)kBTab+hab(r)cab(r)].
(5)

The models are distinguished by different ansatzes for the cross-temperatures

TabSQRT=TaTb,
(6a)
TabMASS=maTb+mbTama+mb,
(6b)

respectively. A distinguishing feature of the SQRT model is that it is mass-independent. In the parameter space of this work (mime,TiTe), it follows that TeiSQRTTeiMASS. Consequently, the SQRT model should be expected to result in stronger interspecies correlations.

In the SVT20 model, the cross-temperature is the same as in the MASS model from Eq. (6b), but the OZ relations are modified to be

ĥab=ĉab+s=i,ens(mabTasmaTabĉasĥsb+mabTsbmbTabĥasĉsb),
(7)

which we will call the SVT–OZ equations.31 Here, mab=mamb/(ma+mb) is the reduced mass of an a, b pair.

The SVT model is based on an ansatz for the two- and three-particle phase-space distribution functions,

Fab(2)=fa(p1)fb(p2)gab(r12),
(8a)
Fabc(3)=fa(p1)fb(p2)fc(p3)gabc(r12,r13,r23),
(8b)

where

fs(p)=(2πmskBTs)32exp(p22mskBTs)
(9)

is the Maxwell–Boltzmann distribution with temperature Ts normalized to unity, and gabc is the triplet distribution function. The cross-temperature, TabMASS naturally arises after integrating the two-particle BBGKY equation over momenta, which gives an Yvon–Born–Green-like equation for the RDFs gab(r12) in terms of the triplet functions gabc(r12,r13,r23).20,23 The SVT–OZ relations are derived by assuming both the superposition approximation for the triplet functions, gabcgabgacgbc, and the HNC approximation from Eq. (5) for the direct correlation functions. Several steps from this point onward are missing from the derivation in Ref. 20. These steps are written out in full in  Appendix A.

In the MASS and SVT models, the interplay between the mass and temperature dependence of Tei is important. From Eq. (6b), one sees that the mass dependence is dominant, causing Tei to rapidly converge to Te for mi20me. From this, one expects the strength of electron-ion correlations in the MASS and SVT models to be similar to that of the electron-electron correlations when the masses are sufficiently different.

The basic screening physics of each model can be understood through the weakly coupled limit. In this limit, ĉabv̂ab/kBTab, and the OZ (or SVT–OZ) equations can be explicitly solved for the partial static structure factors

Sab(k)=δab+nanbĥab(k),
(10)

where δab is the Kronecker delta. The expressions for each model are written in  Appendix B, from which one can compare the models in both k-space and in real space.

First, each model shows qualitative differences in the long-wavelength (k0) limit. The values of Sab(0) from each model are tabulated in Table I. In the long-wavelength limit, the SQRT and SVT model structure factors take finite values, as one would expect of a plasma that exhibits Debye screening. In fact, when memi, both SQRT and SVT give the Sii(k) of a weakly coupled one-component plasma screened by a background species. However, in the MASS model, Sab(0)=0, characteristic of the OCP.32 The physical content of these differences is further elucidated by examining the charge density structure factor, SZZ(k)=12(Sii+2Sei+See). (Note that at weak coupling SZZ is the same whether the electrons are positively or negatively charged due to the leading sign dependence in Sei.) In the long-wavelength limit, SZZ(k) describes variations in the total charge density; the condition that the plasma be quasineutral is equivalent to having SZZ(0)=0. On the other hand, the long-wavelength limit of the partial structure factors Sab(0) describes the screening. Of the models studied here, only SVT satisfies SZZ(0)=0 with nonzero Sab(0). Obviously, the MASS model is quasineutral as well, though Sab(0)=0 suggests that it does so not by self-consistent screening, but by not allowing long-wavelength density variations of any kind. The SQRT model is interesting in that its nonzero Sii(0) implies Debye screening of ions by electrons is included, yet it is not quasineutral (SQRT SZZ(0)0), suggesting that the effect of the ions on the electrons is not consistently treated.

TABLE I.

Long-wavelength limits of the static structure factors for the weakly coupled limit of the SQRT and SVT models, as well as the SVT model when the mass difference is large. In the MASS model (not listed), all Sab(0)=0. The various inverse screening lengths are defined in  Appendix B.

Sab(0)SQRTSVTSVT, memi
ii κe2/κ2 (κ2κei2κe2κi2)/κ2κei2 κe2/κ2 
ei κeκi/κ2 miTi+meTemiTe+meTiκi2κe2/κ2κei2 κe2/κ2 
ee κi2/κ2 (κ2κei2κe2κi2)/κ2κei2 κe2/κ2 
Sab(0)SQRTSVTSVT, memi
ii κe2/κ2 (κ2κei2κe2κi2)/κ2κei2 κe2/κ2 
ei κeκi/κ2 miTi+meTemiTe+meTiκi2κe2/κ2κei2 κe2/κ2 
ee κi2/κ2 (κ2κei2κe2κi2)/κ2κei2 κe2/κ2 

Second, the pole structure of the static structure factors in each model gives rise to different functional forms for the RDFs. In the SQRT model, Sab(k) has a single imaginary pole, which leads to an exponentially screened potential after inverse Fourier transformation. In contrast, the MASS and SVT structure factors have two imaginary poles, so that the screening comes from the difference of two exponentials

gab(r)exp{A1eK1r+A2eK2r4πnanbr}.
(11)

Here, A1, A2, K1, and K2 are constant coefficients. Their values for each model are listed in Table II. Observe that in the limit where memi, the SVT gii(r) is that of a weakly coupled YOCP screened by electrons. The relationship between the SVT model and the screened OCP is expounded upon in Sec. V.

TABLE II.

Coefficients appearing in Eq. (11) for the weak-coupling form for the RDFs for each model, as well as for the SVT model when memi. The various inverse screening lengths are defined in  Appendix B, and c=meTe+miTimiTe+meTi in the SVT column.

SQRTMASSSVTSVT, memi
iieieeiieieeiieieeiieiee
A1 κi2 κiκe κe2 κ+2(κi2κ2)κ+2κ2 κei2κ+2κ+2κ2 κ+2(κe2κ2)κ+2κ2 κi4κ2κei2 κei2κ2cκi2κe2κ2κei2 κe4κ2κei2 κi2 κe2 κe4κi2 
K1 κ κ κ κ+ κ+ κ+ κ κ κ κ κ κ 
A2 κ2(κi2κ+2)κ+2κ2 κei2κ2κ+2κ2 κ2(κe2κ+2)κ+2κ2 κi2(κei2κe2)κ2κei2 κei4cκi2κe2κ2κei2 κe2(κei2κi2)κ2κei2 κe2(κe2κi2)κi2 
K2 … … … κ κ κ κei κei κei … … κe 
SQRTMASSSVTSVT, memi
iieieeiieieeiieieeiieiee
A1 κi2 κiκe κe2 κ+2(κi2κ2)κ+2κ2 κei2κ+2κ+2κ2 κ+2(κe2κ2)κ+2κ2 κi4κ2κei2 κei2κ2cκi2κe2κ2κei2 κe4κ2κei2 κi2 κe2 κe4κi2 
K1 κ κ κ κ+ κ+ κ+ κ κ κ κ κ κ 
A2 κ2(κi2κ+2)κ+2κ2 κei2κ2κ+2κ2 κ2(κe2κ+2)κ+2κ2 κi2(κei2κe2)κ2κei2 κei4cκi2κe2κ2κei2 κe2(κei2κi2)κ2κei2 κe2(κe2κi2)κi2 
K2 … … … κ κ κ κei κei κei … … κe 

Classical molecular dynamics simulations were carried out using the open-source code LAMMPS.33 A two-component, two-temperature plasma was created in a three-dimensional periodic box. The charged particles were made to interact through the repulsive Coulomb potential, Eq. (1), and the long-range part of the Coulomb interaction was accounted for using the particle-particle, particle-mesh method.34 

Every simulation system consisted of 104 particles of each species, each singly charged. The time step for numerical integration was chosen based on the inverse electron plasma frequency, ωpe1=me/4πe2ne. All simulations used time steps in the range δt=0.0050.01ωpe1, which was sufficient to resolve the dynamics of both species.

The equilibration of a two-species system to two different temperatures remains a nontrivial issue from a numerical point of view.35 For the present simulations, each species was coupled to its own Langevin thermostat. The Langevin collision frequencies were chosen such that both species attained their target temperatures within 1% statistical fluctuations. Figure 1 shows how if the thermostat collision frequency was too weak, the ions thermalized to a temperature that was higher than the target temperature. One can see that even a 1% drift from the requested Ti is large enough to make a discernible difference in the ion-ion RDF. We attribute this effect to the fact that in a two-temperature simulation, the thermostats must work against the plasma's natural inclination to thermally relax, which requires that the thermostat collision frequency be greater than the electron-ion collision frequency. If these two rates are comparable, however, then one expects the ions (which couple to the thermostat inefficiently when their mass is large) to thermalize to a temperature greater than the thermostat temperature but less than the temperature they would attain if allowed to relax.

FIG. 1.

Effects of varying the thermostat Langevin collision frequency on the RDFs and temperature fluctuations. For the simulations shown, Γi=50,Γe=1, and mi=30me. The lines show different values of the inverse Langevin collision frequency: ν1=5δt (solid red), 10δt (dashed blue), and 20δt (dash-dotted green).

FIG. 1.

Effects of varying the thermostat Langevin collision frequency on the RDFs and temperature fluctuations. For the simulations shown, Γi=50,Γe=1, and mi=30me. The lines show different values of the inverse Langevin collision frequency: ν1=5δt (solid red), 10δt (dashed blue), and 20δt (dash-dotted green).

Close modal

It was also observed that for high mass ratios, the ion-ion RDF takes much longer to stabilize than the ion temperature. Even after the ions acquire the temperature of their heat bath, Ti, spatial correlations between ions continue to develop for hundreds to thousands of ωpe1 of simulation time. In comparison, gee and gei stabilize on the same timescale as Ti, though small variations thereafter occur in response to the evolution of gii. For reference, the OCP typically requires only a few plasma periods of averaging time for well-resolved RDFs. Since the case of the large mass ratio is of particular experimental importance, the computational burden of simulating such plasmas underscores the need for a reliable theoretical model of correlations in two-temperature plasmas.

We compared the RDFs obtained from a system under Langevin thermostats with those of a system equilibrated using two simultaneous Nosé–Hoover thermostats. At higher mass ratios, the results remain identical irrespective of the choice of thermostat. At lower mass ratios (mi/me5), the system under Nosé–Hoover thermostats displayed the “flying ice cube effect,” in which the system accumulated a spurious net momentum, leading to incorrect RDFs.36 The Langevin thermostats, however, were found to give consistent RDFs for all mass ratios.

The simulations were carried out in three stages. First, we performed an initial thermostatting stage until each species reached its target temperature. The required length of this phase depended on the mass ratio. It was found that for mi = me, 400 electron plasma periods were sufficient and that this number scaled with increased ion mass as mi/me. Second, the evolution of the RDFs was monitored until it was seen that the ion-ion correlations had fully developed. Third, time-averaged RDFs were computed while keeping both the thermostats on. The thermostats were kept active to prevent electron-ion temperature relaxation over the timescales necessary to accurately sample the RDFs. Because the thermostats were left on during the entire simulation period, the total energy was not conserved.

We have evaluated each of the three HNC extensions described in Sec. II and conducted MD simulations as described in Sec. III for several combinations of coupling strengths and mass ratios. Here we present an illustrative subset of the comparisons made, shown in Fig. 2. Plots for other parameter combinations can be found in the supplementary material.

FIG. 2.

Model RDFs compared with molecular dynamics simulation results. The connected black circles are MD, the solid red lines are the SVT model, the dotted orange lines are the MASS model, and the dash-dotted blue lines are the SQRT model.

FIG. 2.

Model RDFs compared with molecular dynamics simulation results. The connected black circles are MD, the solid red lines are the SVT model, the dotted orange lines are the MASS model, and the dash-dotted blue lines are the SQRT model.

Close modal

Figures 2(a) and 2(b) show the radial distribution functions for a plasma of strongly coupled ions and weakly coupled electrons, with mi = me and mi=30me. The first observation to make is that the strength of electron-ion correlations is clearly set by TeiMASS, not by TeiSQRT. The too-wide Coulomb hole in gei(r) shows that the SQRT model overestimates the strength of electron-ion coupling. Furthermore, the SQRT model predicts electron-electron correlation functions that qualitatively differ from MASS, SVT, and MD. The physical reason is most clearly illustrated by examining the weakly coupled limit of the SVT gee(r) when memi. Identifying the potential of mean force as ϕee=kBTelngee, one can write [see Eq. (11) and Table II]

ϕeeSVT(r)e2r[eκerκe2κi2(eκereκr)].
(12)

The first term is the screened repulsion that electrons would experience from one another if they were an OCP, while the “attractive” second term results from the tendency for electrons to cluster when they form screening clouds around ions. These two processes compete, giving rise to the slow decay in the SVT, MASS, and MD gee(r) compared to the SQRT model, which lacks this second “attractive” part. These deficiencies in the SQRT gei(r) and gee(r) were present at all coupling strengths and mass ratios investigated. The errors between SQRT and MD worsen at stronger coupling strengths, as can be seen in the supplementary material.

The remaining comparison of the MASS and SVT models highlights the question of whether the OZ equations require modification to describe a two-temperature system. In all cases studied, the SVT radial distribution functions more closely agree with MD, though the differences between the MASS and SVT RDFs often appear small. In fact, in Ref. 22, the MASS model's apparent accuracy is cited as evidence that SVT's modified OZ equations are unnecessary. Important differences in favor of the SVT approach surface when comparing the structure factors. An example is shown in Fig. 3. The ion-ion structure factor vanishes in the MASS model as k0, indicating that the ions are thermodynamically similar to an unscreened OCP, despite the presence of a screening electron background. In contrast, the SVT model gives a finite value, in line with both MD and the YOCP model. This behavior is demonstrated analytically in Sec. V.

FIG. 3.

Model ion-ion static structure factors compared with molecular dynamics simulation for Γi=4,Γe=0.1, and mi=30me. The inset shows Sii(k) near k = 0, including the YOCP model (green squares).

FIG. 3.

Model ion-ion static structure factors compared with molecular dynamics simulation for Γi=4,Γe=0.1, and mi=30me. The inset shows Sii(k) near k = 0, including the YOCP model (green squares).

Close modal

Since all the models considered are variants of the HNC approximation, it should be expected that they will all suffer inaccuracies at higher coupling due to the lack of bridge functions. In the OCP, bridge functions primarily correct the RDF oscillation amplitudes, which are somewhat too small without the bridge functions. Other differences such as the size of the Coulomb hole and the oscillation phase are relatively minor, so if these features are a point of disagreement between the models and MD, it is more likely due to the two-temperature modeling than due to the lack of bridge functions.

Figure 2(c) shows the RDFs when both species are strongly coupled. As expected, the SVT model underestimates the peak of gii(r) but otherwise agrees well with MD. In contrast, the MASS model appears to break down entirely in this regime of strong electron coupling. An unexpected feature of the MD RDFs is that at high mass ratio, the height of the first peak of gee exceeds that of gei. Ordinarily, one expects the height of this peak to correlate with the strength of the bare interaction compared to the kinetic energy, so that since Ti<Tei<Te, one anticipates max(gii)>max(gei)>max(gee). For low mass ratios, both MD and the HNC models bear out this trend at all coupling strengths, while at higher mass ratios, the HNC models do not capture the augmented first correlation peak in gee observed in MD.

Figure 2(d) shows the breakdown of the two-temperature HNC models at higher ion coupling strength. All three two-temperature models overestimate the strength of correlations in the plasma, exhibiting Coulomb holes and RDF oscillations that are larger than those seen in the MD simulations. This is in contrast to the usual equilibrium HNC theory fails, which underpredicts the peaks. For higher mass ratios and/or lower temperature ratios (see the supplementary material), the SVT RDFs are in surprisingly good agreement with MD even at such strong coupling. These are cases that happen to lie in the transitional regime where SVT goes from underpredicting to overpredicting the RDF peaks.

We now compare the ion-ion correlations of the SVT model to the Yukawa OCP to test the YOCP's limitations as Γe increases. In the classical YOCP model, the electrons are an ideal background that screens the ions. The ions then interact through a Debye-screened potential

viiY(r)=e2reκer,
(13)

where κe=3Γeae1 is the inverse electron Debye length. The YOCP model is valid only when the screening background is weakly coupled, while the SVT model predicts accurate ion-ion RDFs even when Γe exceeds unity. By comparing the YOCP ion-ion RDF gY(r) with gii(r) from two-temperature SVT calculations, we can quantitatively assess at what Γe the YOCP model fails.

A result of the weak-coupling approximation is that κe does not depend on the sign of the electron charge. For this reason, the weak-coupling assumption of the YOCP can be tested using positively charged electrons in the SVT calculations; however, an important caveat must be made. As the electron coupling strength increases, the nature of how they screen the ions is expected to become increasingly dependent on the sign of their charge. It is reasonable to expect, though, that the Γe at which the exponential screening approximation fails is about the same value at which the sign of the electron charge becomes important, since they are both tied to the weak-coupling assumption.

For a given Γi, we solve the HNC–SVT–OZ equations for gii(r) at several Γe and solve the ordinary HNC–OZ equations for gY(r) at several κe. For each Γe, the best-fit κe was chosen to be the one that minimizes the integrated absolute difference between gY and gii from HNC

Δ=dr|gY(r;κe)gii(r)|.
(14)

Figure 4 shows the best-fit YOCP κe over a wide range in Γi and Γe with the mass ratio fixed at mi=1836me. Immediately, one sees that when the electrons are weakly coupled, the best-fit κe is independent of the ion coupling strength and furthermore is essentially the inverse electron Debye length, plotted in black in the figure. The reason becomes clear upon investigating the SVT–OZ equations at weak electron coupling.

FIG. 4.

Inverse screening length of the YOCP whose gY(r) best matches the SVT gii(r) at the same ion coupling strength. The fit criterion is given by Eq. (14). Multiple points at the same Γe are for different values of Γi, indicated by color.

FIG. 4.

Inverse screening length of the YOCP whose gY(r) best matches the SVT gii(r) at the same ion coupling strength. The fit criterion is given by Eq. (14). Multiple points at the same Γe are for different values of Γi, indicated by color.

Close modal

In the limit of weak electron coupling, the Debye–Hückel approximation should be excellent for the electron-electron direct correlation function. Since TeiTe, the same should be true of the electron-ion direct correlation function, giving

ĉee(k)Zi1ĉei(k)4πe2kBTe1k2.
(15)

Due to the large mass ratio, the SVT–OZ equations from Eq. (7) become

ĥii=ĉii+niĥiiĉii+neTeTiĥeiĉei,
(16a)
ĥei=ĉei+niĥiiĉei+neĥeiĉee,
(16b)
ĥee=ĉee+niĥeiĉei+neĥeeĉee.
(16c)

Since ĉei and ĉee are known, ĥei can be eliminated from the first equation to find

ĥii=(ĉii+neĉee1neĉeeTeTiĉei)(1+niĥii).
(17)

If we introduce the notion of the “screened” ion-ion direct correlation function

ĉscr=ĉii+neĉee1neĉeeTeTiĉei,
(18)

then the ion structure factor is given by

Sii(k)=11niĉscr(k),
(19)

meaning that ĉscr mediates a one-to-one mapping between the ion structure of the two-component plasma and that of an equivalent screened one-component plasma.

In the Debye–Hückel approximation for the electrons

niĉscr=niĉii+TiTe1λDi2k211+λDe2k2.
(20)

Now if we decompose ĉii(k) into its singular Coulombic part and a remainder ĉiiR=ĉii+v̂ii/kBTi that is regular as k0,32,37 we find

niĉscr(k)=niĉiiR(k)λDe2λDi211+λDe2k2.
(21)

Thus the long wavelength limit of the ion structure factor is

limk0Sii(k)=11niĉiiR(0)+(λDe/λDi)2.
(22)

Repeating these steps using the ordinary OZ equations results in Eq. (20), but without the factor of Ti/Te, which causes ĉscr to remain singular. This is the reason why the MASS model structure factor is zero in the k0 limit, while the same limit in the SVT model is YOCP-like (nonzero). In passing, it is interesting to note that inserting Eqs. (19) and (20) back into Eq. (16) give the same structure factors found by Boercker and More.19 

Figure 5 demonstrates the breakdown of YOCP behavior when the electrons become strongly coupled. Interestingly, even when Γe1, both the fitted YOCP model and the Debye–Hückel model are in fair agreement with the full two-component SVT calculation. However, further increases to Γe result in ion-ion RDFs that rapidly become non-YOCP-like; even the fitted YOCP underpredicts the ion-ion correlation strength. In other words, the mapping between the two-component system and effective one-component system given by Eq. (18) can no longer be reproduced by an effective Yukawa potential.

FIG. 5.

Comparison of ion-ion RDFs obtained from SVT (solid red), YOCP with fitted κe (dashed blue), and YOCP with κe equal to the inverse electron Debye length (dotted black).

FIG. 5.

Comparison of ion-ion RDFs obtained from SVT (solid red), YOCP with fitted κe (dashed blue), and YOCP with κe equal to the inverse electron Debye length (dotted black).

Close modal

By comparison with molecular dynamics simulations, it has been demonstrated that the model proposed by Seuferling et al.20 accurately extends the Ornstein–Zernike theory of pair correlations to two-temperature plasmas up to and slightly beyond the coupling strengths achieved by present-day ultracold neutral plasma experiments. The assumption of a mass-weighted “cross-temperature” correctly predicts the suppression of electron-ion correlations when the mass ratio is large. Further, we have shown that the modifications made by SVT to the Ornstein–Zernike equations are necessary to give nonzero long-wavelength limits of the static structure factors, which correctly reflects the self-consistent screening of ions by electrons and vice-versa. These findings are given additional weight by our direct comparisons of the ion-ion correlation functions in the SVT and Yukawa OCP models, which indicate that the Yukawa OCP model will become unsuitable even for modeling ion correlations once Γe1.

The present work marks important progress towards a fully two-component description of correlations in classical strongly coupled plasmas. In particular, it suggests that the SVT model can be used to obtain accurate effective scattering potentials or static local field corrections needed in quasi-static descriptions of transport and relaxation processes of strongly coupled plasmas.5,7–11 However there remain interesting physical challenges to overcome. Future work will address the issue of the electron charge, which was taken to be positive in this work to decouple the relevant two-temperature physics from the physics of classical recombination. There is also the question of how to best simulate a two-temperature steady state, both in terms of technical choices regarding thermostats and in terms of the basic statistical mechanics of the simulated ensemble.

See supplementary material for plots of the radial distribution functions and static structure factors for all coupling strengths and mass ratios investigated in this work.

This material is based upon work supported by the National Science Foundation under Grant No. PHY-1453736 and by the Air Force Office of Scientific Research under Award No. FA9550-16-1-0221. 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.

We begin from Eq. (5) of Ref. 23, which after applying the superposition approximation, gabcgabgacgbc, can be written as

r1[kBTablngab+vab]=cncdr3[mabmavacr1mabmbvbcr2]gacgbc,
(A1)

where particles 1 and 2 are of species a and b, respectively (which could be the same or different), and the sum runs over all species labels. Eq. (A1) is a closed set of equations for the RDFs, but it is not suitable for strongly coupled systems because of the use of the superposition approximation. One introduces the direct correlation functions through an HNC-like approximation,

lngab=vabkBTab+habcab,
(A2)

in the hope that the errors from HNC will cancel somewhat the errors made by the superposition approximation. With this and the fact that the lack of external forces implies

cncdr3gacvacr1=cncdr3gbcvbcr2=0.
(A3)

Eq. (A1) becomes

r1[hab(r12)cab(r12)]=cncmabmaTacTabr1[cachbc](r12)+cncmabmbTbcTabr1[haccbc](r12)cncmabmaTacTabdr3hac(r13)γac(r13)r1hbc(r23)+cncmabmbTbcTabdr3hbc(r23)γbc(r23)r2hac(r13),
(A4)

where γ=hc is the indirect correlation function, and the operation denotes convolution. After Fourier transforming r1k,r2k, integrating over k, and dotting k on both sides, one obtains

k2[ĥab(k)ĉab(k)]=k2cncmabmaTacTabĉac(k)ĥbc(k)+k2cncmabmbTbcTabĥac(k)ĉbc(k)cncmabmaTacTabĥbc(k)d(k·)γ̂ac()ĥac(k)+cncmabmbTbcTabĥac(k)d(k·)γ̂bc()ĥbc(k),
(A5)

where is a dummy wavenumber arising from the Fourier transform of a real-space product. Since all the correlation functions must be isotropic in their arguments,

k2[ĥab(k)ĉab(k)]=k2cncmabmaTacTabĉac(k)ĥbc(k)+k2cncmabmbTbcTabĥac(k)ĉbc(k)cncmabmaTacTabĥbc(k)d(k·)γ̂ac()ĥac(|k|)+cncmabmbTbcTabĥac(k)d(k·)γ̂bc()ĥbc(|k|),
(A6)

where we have taken in the last line. The first three lines together form k2 times the SVT–OZ equations as written in Eq. (7), so the remaining two terms must vanish. We abbreviate

zab(k)=k2d(k·)γ̂ab()γ̂ab(|k|),
(A7)

and call the last two terms of Eq. (A6) the “remainder,” Rab, so that Eq. (A6) may be written as

ĥab=ĉab+cncmabmaTacTabĉacĥbc+cncmabmbTbcTabĥacĉbc+Rab
(A8)

with

Rab=cncmabmaTacTabĥbc(k)ẑac(k)+cncmabmaTacTabĥac(k)ẑbc(k).
(A9)

For the like-species equation (a = b), Rab vanishes trivially. For the cross-species equation (ab), observe that Rab changes sign upon interchange of species labels (ab), while the other terms of Eq. (A8) do not. Therefore, Rab = 0 for all combinations of a and b, giving Eq. (7).

In the limit of weak coupling, the direct correlation functions may be approximated

ĉab(k)4πe2kBTab1k2,
(B1)

and it is straightforward to solve each model for the static stucture factors, Sab(k)=δab+nanbĥab(k), in terms of various characteristic screening lengths. Using the notation

κi2=4πe2ni/kBTiκe2=4πe2ne/kBTeκei2=4πe2nine/kBTeimassκ2=κe2+κi2κ±2=κ22±κ44κi2κe2+κei4,

one finds for the SQRT model

Sii=k2+κe2k2+κ2,
(B2a)
Sei=κeκik2+κ2,
(B2b)
See=k2+κi2k2+κ2
(B2c)

for the MASS model

Sii=k4+κe2k2(k2+κ+2)(k2+κ2),
(B3a)
Sei=κei2k2(k2+κ+2)(k2+κ2),
(B3b)
See=k4+κi2k2(k2+κ+2)(k2+κ2)
(B3c)

and for the SVT model

Sii=(k2+κei2)(k2+κe2)+κi2(κei2κe2)(k2+κ2)(k2+κei2),
(B4a)
Sei=κei2k2miTi+meTemiTe+meTiκe2κi2(k2+κ2)(k2+κei2),
(B4b)
See=(k2+κei2)(k2+κi2)+κe2(κei2κi2)(k2+κ2)(k2+κei2).
(B4c)
1.
T. C.
Killian
,
S.
Kulin
,
S. D.
Bergeson
,
L. A.
Orozco
,
C.
Orzel
, and
S. L.
Rolston
,
Phys. Rev. Lett.
83
,
4776
(
1999
). physics/9908051.
2.
S. J.
Putterman
and
K. R.
Weninger
,
Annu. Rev. Fluid Mech.
32
,
445
(
2000
).
3.
S.
Bergeson
and
M.
Kleinert
, in
APS Meeting Abstracts
(
2016
).
4.
H. G.
Rinderknecht
,
M. J.
Rosenberg
,
C. K.
Li
,
N. M.
Hoffmann
,
G.
Kagan
,
A. B.
Zylstra
,
H.
Sio
,
J. A.
Frenje
,
M. G.
Johnson
,
F. H.
Séguin
,
R. D.
Petrasso
,
P.
Amendt
,
C.
Bellei
,
S.
Wilks
,
J.
Delettrez
,
V. Y.
Glebov
,
C.
Stoeckl
,
T. C.
Sangster
,
D. D.
Meyerhofer
, and
A.
Nikroo
,
Phys. Rev. Lett.
114
,
025001
(
2015
).
5.
S. D.
Baalrud
and
J.
Daligault
,
Phys. Rev. Lett.
110
,
235001
(
2013
).
6.
N. R.
Shaffer
,
S. D.
Baalrud
, and
J.
Daligault
,
Phys. Rev. E
95
,
013206
(
2017
).
7.
S.
Ichimaru
,
S.
Mitake
,
S.
Tanaka
, and
X.-Z.
Yan
,
Phys. Rev. A
32
,
1768
(
1985
).
8.
M. W. C.
Dharma-Wardana
and
F.
Perrot
,
Phys. Rev. E
58
,
3705
(
1998
).
9.
J.
Daligault
and
G.
Dimonte
,
Phys. Rev. E
79
,
056403
(
2009
).
10.
J.
Vorberger
,
D. O.
Gericke
,
T.
Bornath
, and
M.
Schlanges
,
Phys. Rev. E
81
,
046404
(
2010
).
11.
L. X.
Benedict
,
M. P.
Suhr
,
L. G.
Stanton
,
C. R.
Scullard
,
A. A.
Correa
,
J. I.
Casto
,
F. R.
Graziani
,
L. A.
Collins
,
O.
Čertík
,
J. D.
Kress
, and
M. S.
Murillo
,
Phys. Rev. E
95
,
043202
(
2017
).
12.
S. K.
Tiwari
,
N. R.
Shaffer
, and
S. D.
Baalrud
,
Phys. Rev. E
95
,
043204
(
2017
).
13.
P.
Debye
and
E.
Hückel
,
Phys. Z.
24
,
185
(
1923
).
14.
J. P.
Hansen
and
I. R.
MacDonald
,
Theory of Simple Liquids
, 1st ed. (
Academic Press
,
1976
).
15.
S. D.
Bergeson
,
A.
Denning
,
M.
Lyon
, and
F.
Robicheaux
,
Phys. Rev. A
83
,
023409
(
2011
).
16.
T. S.
Strickler
,
T. K.
Langin
,
P.
McQuillen
,
J.
Daligault
, and
T. C.
Killian
,
Phys. Rev. X
6
,
021021
(
2016
).
17.
A.
Melzer
,
S.
Nunomura
,
D.
Samsonov
,
Z. W.
Ma
, and
J.
Goree
,
Phys. Rev. E
62
,
4162
(
2000
).
18.
E. E.
Salpeter
,
J. Geophys. Res.
68
,
1321
, doi: (
1963
).
19.
D. B.
Boercker
and
R. M.
More
,
Phys. Rev. A
33
,
1859
(
1986
).
20.
P.
Seuferling
,
J.
Vogel
, and
C.
Toepffer
,
Phys. Rev. A
40
,
323
(
1989
).
21.
R.
Bredow
,
T.
Bornath
,
W.-D.
Kraeft
, and
R.
Redmer
,
Contrib. Plasma Phys.
53
,
276
(
2013
).
22.
M. W. C.
Dharma-Wardana
and
M. S.
Murillo
,
Phys. Rev. E
77
,
026401
(
2008
).
23.
V.
Schwartz
,
T.
Bornath
,
W.-D.
Kraeft
,
S. H.
Glenzer
,
A.
Höll
, and
R.
Redmer
,
Contrib. Plasma Phys.
47
,
324
(
2007
).
24.
D. V.
Rose
,
T. C.
Genoni
,
D. R.
Welch
,
R. E.
Clark
,
R. B.
Campbell
,
T. A.
Mehlhorn
, and
D. G.
Flicker
,
Phys. Plasmas
16
,
102105
(
2009
).
25.
L. S.
Ornstein
and
F.
Zernike
,
Proc. K. Ned. Akad. Wet.
17
,
793
807
(
1914
).
26.
J. M. J.
van Leeuwen
,
J.
Groeneveld
, and
J.
de Boer
,
Physica
25
,
792
(
1959
).
27.
Y.
Rosenfeld
and
N. W.
Ashcroft
,
Phys. Rev. A
20
,
1208
(
1979
).
28.
H.
Iyetomi
and
S.
Ichimaru
,
Phys. Rev. A
27
,
3241
(
1983
).
29.
H.
Iyetomi
,
S.
Ogata
, and
S.
Ichimaru
,
Phys. Rev. A
46
,
1051
(
1992
).
30.
G.
Kahl
,
B.
Bildstein
, and
Y.
Rosenfeld
,
Phys. Rev. E
54
,
5391
(
1996
).
31.

The original formulas in Eq. (38) of Ref. 20 contain some typographical errors.

32.
M.
Baus
and
J.-P.
Hansen
,
Phys. Rep.
59
,
1
(
1980
).
33.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
34.
S.
Plimpton
,
R.
Pollock
, and
M.
Stevens
, in
Proceedings of the 8th SIAM Conference on Parallel Processing for Scientific Computing
(
1997
).
35.
D.
Fukushi
,
K.
Mae
, and
T.
Honda
,
Jpn. J. Appl. Phys.
39
,
5014
(
2000
).
36.
S. C.
Harvey
,
R. K.-Z.
Tan
, and
T. E.
Cheatham
,
J. Comput. Chem.
19
,
726
(
1998
).
37.
M.
Baus
,
J. Phys. A
11
,
2451
(
1978
).

Supplementary Material