Pair correlation functions of strongly coupled two-temperature plasma

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 potentials^5,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 plasmas¹² 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.⁶ 

At weak coupling, the pair distribution functions are accurately described by the Debye–Hückel theory of electrolytes.¹³ 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.¹⁴ 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.¹⁸ 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.¹⁹ 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 T_ab 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, T_ab, or whether it is necessary to modify the Ornstein–Zernike relations. This work considers three formulations^20–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 $Te≥Ti$ and $me≤mi$⁠, 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

and the Coulomb coupling strength of each species is

where $as=(3/4πns)1/3$ is the mean spacing between particles of species s, T_s is their temperature, e is the elementary charge, and k_B 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 T_i and T_e.

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 are²⁵ 

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 by²⁶ 

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” T_ab 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 SQRT²¹ and MASS²² models, the OZ relations are taken to be the same as at equilibrium, and the HNC closures are assumed to be

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

respectively. A distinguishing feature of the SQRT model is that it is mass-independent. In the parameter space of this work (⁠$mi≥me, Ti≤Te$⁠), it follows that $TeiSQRT≤TeiMASS$⁠. Consequently, the SQRT model should be expected to result in stronger interspecies correlations.

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

which we will call the SVT–OZ equations.³¹ 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,

where

is the Maxwell–Boltzmann distribution with temperature T_s normalized to unity, and g_abc 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, $gabc≈gabgacgbc$⁠, 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 T_ei is important. From Eq. (6b), one sees that the mass dependence is dominant, causing T_ei to rapidly converge to T_e for $mi≳20me$⁠. 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, $ĉab≈−v̂ab/kBTab$⁠, and the OZ (or SVT–OZ) equations can be explicitly solved for the partial static structure factors

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 (⁠$k→0$⁠) 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 $me≪mi$⁠, 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.³² 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 S_ZZ is the same whether the electrons are positively or negatively charged due to the leading sign dependence in S_ei.) 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.

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

Here, A₁, A₂, K₁, and K₂ are constant coefficients. Their values for each model are listed in Table II. Observe that in the limit where $me≪mi$⁠, 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.

Classical molecular dynamics simulations were carried out using the open-source code LAMMPS.³³ 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.³⁴ 

Every simulation system consisted of 10⁴ particles of each species, each singly charged. The time step for numerical integration was chosen based on the inverse electron plasma frequency, $ωpe−1=me/4πe2ne$⁠. All simulations used time steps in the range $δt=0.005–0.01ωpe−1$⁠, 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.³⁵ 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 T_i 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.

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, T_i, spatial correlations between ions continue to develop for hundreds to thousands of $ωpe−1$ of simulation time. In comparison, g_ee and g_ei stabilize on the same timescale as T_i, though small variations thereafter occur in response to the evolution of g_ii. 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/me≲5$⁠), 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.³⁶ 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 m_i = m_e, 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.

Figures 2(a) and 2(b) show the radial distribution functions for a plasma of strongly coupled ions and weakly coupled electrons, with m_i = m_e 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 $me≪mi$⁠. Identifying the potential of mean force as $ϕee=−kBTe ln gee$⁠, one can write [see Eq. (11) and Table II]

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 $k→0$⁠, 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.

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 g_ee exceeds that of g_ei. 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 g_ee 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

where $κe=3Γeae−1$ 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 g_ii from HNC

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.

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

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

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

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

then the ion structure factor is given by

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

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

Thus the long wavelength limit of the ion structure factor is

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 $k→0$ 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.¹⁹ 

Figure 5 demonstrates the breakdown of YOCP behavior when the electrons become strongly coupled. Interestingly, even when $Γe≃1$⁠, 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.

By comparison with molecular dynamics simulations, it has been demonstrated that the model proposed by Seuferling et al.²⁰ 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 $Γe≳1$⁠.

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, $gabc≈gabgacgbc$⁠, can be written as

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,

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

Eq. (A1) becomes

where $γ=h−c$ is the indirect correlation function, and the $⋆$ operation denotes convolution. After Fourier transforming $r1→k,r2→k′$⁠, integrating over $k′$⁠, and dotting $k$ on both sides, one obtains

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,

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

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

with

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

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

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

one finds for the SQRT model

for the MASS model

and for the SVT model