Isomorph invariance and thermodynamics of repulsive dense bi-Yukawa one-component plasmas

In numerous realizations of complex plasmas, dust-dust interactions are characterized by two screening lengths and are thus better described by a combination of Yukawa potentials. The present work investigates the static correlations and the thermodynamics of repulsive dense bi-Yukawa fluids based on the fact that such strongly coupled systems exhibit isomorph invariance. The strong virial-potential energy correlations are demonstrated with the aid of molecular dynamics simulations, an accurate analytical expression for the isomorph family of curves is obtained and an empirical expression for the fluid-solid phase-coexistence line is proposed. The isomorph-based empirically modified hypernetted-chain approach, grounded on the ansatz of isomorph invariant bridge functions, is then extended to such systems and the resulting structural properties show an excellent agreement with the results of computer simulations. A simple and accurate closed-form expression is obtained for the excess internal energy of dense bi-Yukawa fluids by capitalizing on the compact parameterization offered by the Rosenfeld-Tarazona decomposition in combination with the Rosenfeld scaling, which opens up the energy route to thermodynamics.


I. INTRODUCTION
Complex (or dusty) plasmas can be defined as weakly ionized plasmas that are dispersed with solid charged particulates of mesoscopic size. 1,2 For large dust densities and high charges, the system can be considered as one-component, since the dust species is dominant in terms of energy and momentum transfer. 1,2 In the spatiotemporal scales relevant for dust dynamics, plasma assumes the role of the neutralizing background responsible for the quasi-stationary dust charge and the shielding of the dust-dust interaction. Such systems are typically engineered in the laboratory and, with proper manipulation, can virtually attain all states of matter. [2][3][4] In the dense fluid region of the complex plasma phase diagram, the dust interaction energy is much higher than the dust kinetic energy but not high enough to drive crystallization.
In strongly coupled complex plasmas, regardless of dimensionality, the dust-dust interaction energy has been generally considered to conform to the Debye-H€ uckel or Yukawa type, 2-5 uðrÞ ¼ ðQ 2 =rÞ exp ðÀr=k D Þ, where Q is the dust charge and k D is the plasma screening length. This was supported by early measurements of the interaction between dust particles levitating in the sheath of a radio-frequency discharge 6 and corresponds to the paradigm of the Yukawa onecomponent plasma (YOCP) which reduces to the classical onecomponent plasma (OCP) in the absence of screening (k D ! 1).
Yukawa interactions result from the linearized Boltzmann response of isotropic plasmas around a test point charge. 2 However, even in the absence of plasma drifts, the continuous absorption of plasma fluxes on the particulate surface leads to non-Boltzmann densities for isolated dust [7][8][9] and brings forth the issue of ionizationrecombination balance for dense dust clouds. 10,11 In the latter, simple hydrodynamic models have predicted that the interaction potential acquires a double Yukawa repulsive form, 12,13 where the short-range screening length is controlled by plasma shielding and the long-range screening length is controlled by plasma sink-source competition. The possibility of long-range attraction has also been discussed in the literature. [14][15][16] We shall refer to such strongly-interacting systems as bi-Yukawa one-component plasmas (biYOCPs). It should be emphasized that the biYOCP interaction potential differs from the prototype two-Yukawa or double-Yukawa potential used to describe pair interactions in colloid-polymer mixtures, which also involves a hard-core repulsion. [17][18][19] The general form of the bi-Yukawa pair-interaction potential is uðrÞ ¼ ð 1 =rÞ exp ðÀr=k 1 Þ þ ð 2 =rÞ exp ðÀr=k 2 Þ. We shall limit ourselves to purely repulsive interactions ( 1 , 2 > 0) with a pure Coulomb short-range asymptote ( 1 þ 2 ¼ Q 2 ) and a dominant short-range screening term (k 2 ! k 1 , 2 1 ). We will focus exclusively on extended three-dimensional dense biYOCP liquids. Expressed in reduced coordinates, x ¼ r/d, where d ¼ (4pn/3) À1/3 is the Wigner-Seitz radius with n being the particle density, the biYOCP interaction potential becomes In the above b ¼ 1/(k B T), where T is the particle temperature and k B is the Boltzmann constant, the dimensionless thermodynamic state variables of the biYOCP are the coupling parameter defined by C ¼ bQ 2 /d and the normalized screening parameter defined by j ¼ d/k 1 . The external potential parameters of the biYOCP are the relative strength r ¼ 2 /Q 2 and the relative screening length l ¼ k 1 /k 2 subject to the constraints 0 r 0.5 and 0 l 1. We shall assume that, in contrast to the state variables (C, j), the external parameters (r, l) depend on neither the density nor the temperature of the dust component. These parameters should vary in different experimental realizations, since they depend on the discharge conditions as well as on the dust size and composition. It is evident that the biYOCP system reduces to the YOCP for r ¼ 0 and that it reduces to the OCP for j ¼ 0 In spite of the abundance of theoretical studies that predict a bi-Yukawa interaction potential for many laboratory realizations of complex plasmas, [12][13][14][15][16] there have been few investigations of the thermodynamic, static, and dynamic properties of the strongly coupled biYOCP. [20][21][22] This probably stems from the fact that the two state variables and two external potential parameters form a four-dimensional (C, j; r, l) space that is difficult to cover systematically. However, given the exhaustive study of the YOCP during the last decades as well as the anticipation that the dense biYOCP exhibits approximate hidden scale invariance and possibly obeys crucial additivity theorems, this number of independent parameters does not preclude the existence of simple closed-form relations for some fundamental aspects of dense biYOCP liquids.
The aim of the present paper is to provide simple analytical expressions that accurately describe the excess internal energy of dense repulsive bi-YOCP liquids for arbitrary (C, j; r, l) parameters, which also opens up the energy route to all thermodynamic properties. The treatment is based on calculating the excess internal energy with the recently proposed isomorph-based empirically modified hypernettedchain (IEMHNC) approximation 23 and on taking advantage of the compact parameterization offered by the Rosenfeld-Tarazona decomposition 24 in combination with the Rosenfeld scaling. 25 The methodology should be highly accurate provided that dense repulsive biYOCP liquids are Roskilde simple systems, i.e., that they possess continuous paths of constant excess entropy in their phase diagram along which a large number of properties in reduced units is invariant to a good approximation. [26][27][28][29][30][31][32][33][34][35] Hence, the strong virialpotential energy correlations of the dense biYOCP are first demonstrated with the aid of molecular dynamics (MD) simulations, an analytical expression for the biYOCP isomorph family of curves is obtained and an empirical expression for the biYOCP melting line is introduced. These characteristics are also justified in light of their YOCP counterparts with the aid of approximate additivity theorems. The IEMHNC integral approach, based on the conjecture of isomorph invariant bridge functions, is then extended to biYOCP systems and the resulting structural properties are compared to molecular dynamics results. The Rosenfeld-Tarazona decomposition of the reduced excess internal energy 24 is finally generalized to dense biYOCP liquids; the static component is computed from the asymptotic high-density limit with hard-sphere Percus-Yevick radial distribution functions, while the thermal component is computed from the IEMHNC approach and fitted to the Rosenfeld scaling. 25

II. ISOMORPH THEORY
Isomorphs are curves in the phase diagram along which a large set of the structural but also dynamic properties are approximately invariant when expressed in properly reduced units, [26][27][28][29][30][31][32][33][34][35] in which the length is normalized to the mean-cubic interparticle distance a ¼ n À1/3 , the energy to k B T and the time to s ¼ n À1=3 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m=ðk B TÞ p , where m is the particle mass. For one-component systems, isomorphs are of particular interest because they effectively transform the phase diagram from two-dimensional to one-dimensional. Such curves exist for a rather sizable class of condensed matter systems, which are known as Roskilde-simple or R-simple.

A. Definition of R-simple systems
We first introduce the collective N-particle position vector configuration R ¼ ðr 1 ; r 2 ; …; r N Þ, the total potential energy UðRÞ, and the microscopic virial WðRÞ ¼ ð1=3Þ P i r i Á F i ¼ Àð1=3ÞR Á rUðRÞ. R-simple systems are systems that are characterized by strong correlations between their virial and their potential energy fluctuations in constant-volume (NVT) equilibrium, which are quantified via the standard Pearson coefficient defined as R WU ¼ hDWDUi= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hðDWÞ 2 ihðDUÞ 2 i q , where h…i denotes canonical ensemble averaging and DA ¼ A À hAi statistical fluctuations around the thermodynamic mean. The strong W-U correlations are empirically delimited by the practical condition R WU տ 0:9. 26-29 Such a definition enables a straightforward identification of R-simple systems with the aid of NVT Molecular Dynamics (MD) simulations, 27,28 but allows neither the rigorous derivation of their fundamental properties nor the numerical determination of the isomorph curves. We shall refer to it as the "practical definition" of R-simple systems.
For a more concrete characterization, the so-called "axiomatic definition" of R-simple systems has been formulated, 31,32 according to which an R-simple system is defined by the relation where the equilibrium configurations R a ; R b correspond to the same density (not necessarily configurations of the same temperature) with k being an arbitrary positive uniform scaling factor (though in practice of the order of unity as we will discuss in what follows). The validity of Eq. (4) for all possible equal density configurations implies perfect Physics of Plasmas ARTICLE scitation.org/journal/php W-U correlations (R WU 1) which is strictly true only for Euler homogeneous pair-interactions such as those realized in inversepower law (IPL) systems. Consequently, for typical R-simple systems, Eq. (4) should always be understood to be valid in a more restricted sense, i.e., for most physically relevant configurations. It is evident that both definitions depend on the thermodynamic state and should be satisfied in an extended region of the system's phase diagram. Overall, R-simple systems are usually identified according to the practical definition and it is accepted that all properties derived within the axiomatic definition will be followed in some approximate fashion. Three basic properties of R-simple systems are the following: [31][32][33][34] (1) Single phases possess isomorphs and the isomorph curves constitute configurational adiabats, lines of constant excess entropy, of the respective phase region. (2) Phase coexistence curves are themselves isomorphs to the first-order; small deviations from this have recently been predicted using the isomorph theory. 34 (3) Two arbitrary configurations R 1 ; R 2 that are identical in reduced units belong to thermodynamic state points on the same isomorph and have approximately equal probabilities in the canonical ensemble, i.e., where the C 12 factor depends on the two thermodynamic state points but not on the microscopic configurations.

B. Direct isomorph check
Different methods have been developed for the tracing of the isomorph curves of R-simple systems via MD simulations. 27,30,33 The socalled direct isomorph check is based on the equal canonical probabilities and works as follows. 27 Equation (5) can be rewritten as Since these two configurations should be identical in reduced units, we also have that n Let us suppose a (n 1 , T 1 ) reference state and an isomorphic state (n 2 , T 2 ) of rescaled density n 2 ¼ (1/k)n 1 . The goal is to determine T 2 . The potential energy UðR 1 Þ is first extracted by an equilibrium (n 1 , T 1 ) simulation, the configuration is then rescaled to R 2 ¼ kR 1 and the potential energy UðR 2 Þ is obtained. This procedure is repeated for several configurations and a UðR 1 Þ; UðR 2 Þ scatter plot is obtained. As a consequence of Eq. (6), the points should be well-approximated by a straight line whose slope h can be obtained with linear regression. This allows for the determination of T 2 from T 2 ¼ hT 1 . It is evident that an isomorph line consisting of M points can be simply traced by applying the above procedure M times, always starting from the newly obtained isomorphic point.
Theoretically, the direct isomorph check should work for arbitrary density jumps, i.e., for any positive scaling factor k. However, the approximate nature of isomorph invariance for non-IPL systems and thus of Eq. (6) should be taken into account. This imposes a k Շ 1:1 restriction, which indicates that tracing is usually performed with density increments ðn 1 À n 2 Þ=n 1 Շ 10%. 33 We note that, for some interaction potentials, it is more practical to carry out MD simulations directly in reduced units by setting the temperature and density equal to unity. The direct isomorph check is then performed by controlling the length and energy parameters of the potential instead of the density and temperature. 33

III. ISOMORPHIC ASPECTS OF THE YOCP STRUCTURE AND THERMODYNAMICS
We recapitulate basic isomorphic aspects concerning key structural and thermodynamic YOCP properties, which will later be generalized to biYOCP systems. The strong W-U equilibrium correlations of the YOCP are documented, the analytical expression of the YOCP isomorph family of curves is presented, the isomorph-based empirically modified hypernetted-chain approximation is introduced and the Rosenfeld-Tarazona decomposition of the internal energy is discussed.

A. Strong correlations and isomorph curves
A recent dedicated MD investigation has revealed that Yukawa liquids are R-simple systems providing a general framework for the understanding of earlier discovered invariant YOCP aspects. 33 The YOCP exhibits a high correlation coefficient R WU > 0.98 for an extended region of its liquid phase and YOCP isomorphs are accurately described by an analytical expression of the form 33 where a ¼ a/d ¼ (4p/3) 1=3 is the ratio of the cubic mean interparticle distance to the Wigner-Seitz radius. This expression is equivalent to the semiempirical relation earlier proposed for the YOCP melting line 36,37 that well describes the extensive MD liquid-solid phase coexistence data available in the literature 38,39 In the above, C OCP m ¼ 171:8 is the coupling parameter at the OCP melting point. 38

B. Isomorph-based integral theory approach
For isotropic one-component systems, the integral equation theory of liquid structure comprises of the Ornstein-Zernike equation 40 hðrÞ ¼ cðrÞ þ n ð cðr 0 Þhðjr À r 0 jÞd 3 r 0 ; combined with the exact closure equation 40 gðrÞ ¼ exp ÀbuðrÞ þ hðrÞ À cðrÞ þ BðrÞ ½ ; where h(r) ¼ g(r) À 1 is the total correlation function, g(r) is the radial distribution function, c(r) is the direct correlation function, and B(r) is the unknown bridge function. Different theoretical approaches employ different approximations for the bridge function, for instance B(r) ¼ 0 corresponds to the hypernetted-chain (HNC) approximation and B(r) ¼ ln[1 þ c(r)] À c(r) corresponds to the Percus-Yevick (PY) approximation, 40 where we introduced the indirect correlation function c(r) ¼ h(r) À c(r) for convenience. The isomorph-based empirically modified hypernetted-chain (IEMHNC) approximation is based on the conjecture that bridge functions remain completely invariant along isomorph curves. 23 For Physics of Plasmas ARTICLE scitation.org/journal/php strongly coupled Yukawa systems, this recently developed approach assumes that the YOCP bridge function B YOCP (r, C, j) is an isomorph invariant up to and including the OCP limit and takes advantage of the availability of simulation-extracted OCP bridge functions B OCP (r, C). From Eq. (7), we have Accurate closed-form expressions are available for the OCP bridge function B OCP (r, C) that have been obtained with the aid of Monte Carlo simulations and read as 41 where the coefficients b 0 (C), b 1 (C), c 1 (C), c 2 (C), and c 3 (C) are described by Eqs. (21) and (23)

C. The Rosenfeld-Tarazona decomposition of the internal energy
On the basis of density functional theory and thermodynamic perturbation theory around the hard sphere reference system, by retaining the leading order terms of an asymptotic high-density expansion, Rosenfeld and Tarazona obtained a simple general decomposition of the excess internal energy of fluids valid at high densities and for predominant repulsive interactions. 24 In standard (n, T) state variables, the RT decomposition reads as which leads to the well-known T À2/5 RT scaling for the excess isochoric heat capacity C ex V ¼ ð@U ex =@TÞ V . The RT scaling has been demonstrated to be closely followed by R-simple liquids. 46 It is worth pointing out that the original formulation of isomorph theory proved that C ex V is isomorph invariant, 27 but this does not apply in the more generic reformulation 31 and has been confirmed not to apply for all Rsimple liquids. 47,48 For the YOCP system, when expressed in (C, j) state variables and via reduced excess internal energies u ex ¼ U ex /(Nk B T), the RT decomposition reads as 24 u ex ðC; jÞ ¼ u st ðC; jÞ þ u th ðC; jÞ ; where u ex ðC; jÞ ¼ ð3=2Þ The "static component" u st refers to the excess internal energy of the system with its constituents frozen in a regular structure associated with the infinite-coupling limit. In this unphysical limit, the static correlations for repulsive pair additive systems coincide with the unitary packing fraction limit of the PY approximation for hard spheres. 24 For the YOCP, we have with g ¼ ðp=6Þnr 3 hs is the packing fraction, r hs is the hard sphere radius, and y ¼ r/r hs . The integral involving the PY hard-sphere radial distribution function g PY (y; g) has an analytical form 49,50 that, in the limit g ! 1, yields which has a removable divergence in the OCP j ! 0 limit u YOCP st ðC; jÞ ¼ ð3=2ÞC=j 2 À 9C=10 þ OðjÞ that cancels out exactly with the background contribution leading to the well-known OCP result u OCP st ¼ À9C=10. 51,52 Note that Eq. (16) can also be derived within the Ion-Sphere Model (ISM) 51,52 extended to the YOCP. 53,54 The ISM approach is based on purely electrostatic arguments and offers a physically transparent proof but does not allow for straightforward generalizations to more complicated interaction potentials.
The "thermal component" u th accounts for deviations from fixed positions as caused by thermal motion. It can be acquired from Eq. (13), that for IPL and Yukawa potentials, we have u th / C 2=5 , 24 We point out that the j-dependence of the thermal component stems exclusively from the melting line C YOCP m ðjÞ. In view of Eqs. (7) and (8)

IV. ISOMORPHIC ASPECTS OF THE BIYOCP STRUCTURE AND THERMODYNAMICS A. Strong correlations and isomorph curves
It can be argued that if u 1 (r), u 2 (r) are the pair potentials of two R-simple liquids, then the pair potential u 1 (r) þ u 2 (r) should also correspond to a R-simple liquid. 29,55 This additivity theorem, based on the quasi-universality of simple liquids, 32 clearly suggests that the dense biYOCP system is R-simple. Furthermore, Eq. (5) suggests that the isomorph curves are also additive when expressed in (C, j) state variables. Therefore, the biYOCP isomorph curves will be given by C biYOCP iso ðC; j; r; lÞ ¼ C YOCP iso ½Cð1 À rÞ; j þ C YOCP iso ðCr; ljÞ. Courtesy of Eq. (7), this leads to the analytical expression

Physics of Plasmas
As emphasized in Sec. II A, the properties of R-simple systems are satisfied in an approximate manner for non-IPL potentials. Hence, MD simulations should be carried out in order to verify that biYOCP systems are indeed R-simple and to quantify the accuracy of the analytical expression for the biYOCP isomorphs. The MD simulations were carried out on graphics cards with the RUMD open-source software. 58 The simulations were performed directly in reduced units with N ¼ 8192 particles and a time-step of Dt ¼ 2.5 Â 10 À3 s. The interaction potential was truncated at the cutoff distance r cut ¼ 10a with the shifted-force cutoff method. 59 Since the Ewald decomposition is not implemented in RUMD, long-range interactions close to the unscreened Coulomb limit cannot be handled accurately 60,61 and we confined ourselves to j ! 1. We also focused on the dense fluid region of the phase diagram, 0:1 < C=C biYOCP m < 1.
The biYOCP system exhibited a very high correlation coefficient for all state points studied. The strong correlations between the virial and potential energy fluctuations are documented in Fig. 1 for r ¼ 0.10, l ¼ 0.15 and two thermodynamic states, both characterized by R WU > 0.99. Overall, it is concluded that biYOCP liquids are Rsimple for any r < 0.5, l < 1 at least in the phase region defined by 0:1 < C=C biYOCP m < 1 and j ! 1. The direct isomorph check with a screening parameter increment Dj/j ¼ 3.5% was employed in order to trace isomorphs in the interval 1.0 j 5.0. The screening parameters of the state points that are generated by the direct isomorph check belong to a geometric progression with a common ratio of 1.035. The initial step of the numerical procedure is illustrated in Fig. 2 and detailed in the respective caption, see also Sec. II B. The numerical isomorph curves are illustrated in Fig.  3 together with the analytical expression of Eq. (18). We considered three reduced coupling parameters C=C biYOCP m ¼ 0:1; 0:5; 0:8 f gand two (r, l) combinations: r ¼ 0.10, l ¼ 0.15 in Fig. 3(a) and r ¼ 0.20, l ¼ 0.25 in Fig. 3(b).
The closed-form description of the isomorph lines, see Eq. (18), displays a remarkable accuracy when compared with the discrete numerical isomorph points. The mean absolute relative deviations between the analytical and numerical isomorphs never exceeded 3.3%, thus verifying the additivity theorems. In the same figure, the melting line for the biYOCP has been compared to the melting line for the YOCP as obtained from MD simulations 39 and from the analytical expression of Eq. (8). As expected from the longer range of the FIG. 1. Scatter plot demonstrating the strong W-U correlations of the biYOCP system for r ¼ 0.10 and l ¼ 0.15 at the state points j ¼ 2, C ¼ 0.1C m and j ¼ 2, C ¼ 0.8C m . The virial and potential energy equilibrium fluctuations have been normalized by the particle number N; j ¼ 2, C ¼ 0.1C m is characterized by R WU ¼ 0.996 (green circles) and j ¼ 2, C ¼ 0.8C m by R WU ¼ 0.998 (red squares) revealing that the biYOCP is R-simple at both thermodynamic states.

FIG. 2.
Direct isomorph check for the biYOCP system when r ¼ 0.10 and l ¼ 0.15. Scatter plot of the potential energy of two initial configurations vs the potential energy of two respective scaled configurations. The two initial configurations correspond to the state points j 1 ¼ 1.0, C 1 ¼ 0.5C m (cyan circles) & j 1 ¼ 1.0, C 1 ¼ 0.8C m (orange squares); for both cases, the scaled screening parameter is j 2 ¼ 1.035. The temperature T 2 of the scaled configurations is evaluated by multiplying the linear regression slope h with the temperature T 1 of the initial configurations, i.e., T 2 ¼ hT 1 . The coupling parameter is then obtained from C 2 ¼ (j 1 /j 2 )(C 1 /h).

ARTICLE
scitation.org/journal/php interactions, for the same screening parameter, the melting point of the biYOCP is always smaller than the melting point of the YOCP, with the obvious exception of j ¼ 0 where both systems reduce to the OCP.

B. Isomorph-based integral theory approach
The results of Sec. IV A allow us to extend the IEMHNC approach to biYOCP systems. Taking into account the fact that the biYOCP reduces to the OCP for j ! 0 and invoking the ansatz of isomorph-invariant bridge functions, we can employ the biYOCP isomorph mapping introduced in Eq. (18) in order to acquire a correspondence relation between the unknown biYOCP bridge function and the known OCP bridge function, 23 The range of validity of the simulation-extracted OCP bridge function, Eq. (12), dictates the validity range of the constructed biYOCP bridge function, Eq. (20), which becomes 5:25 C biYOCP iso ðC; j; r; lÞ 171:8. The biYOCP radial distribution functions obtained from the IEMHNC approach and MD simulations will be compared to quantify the accuracy of the integral theory. The numerical procedure followed for the solution of the system of Eqs. (9), (10), (12), and (20) is based on Picard iterations in Fourier space combined with standard mixing algorithms and established long-range decomposition schemes (when necessary) in order to improve convergence. 23 The upper range cut-off was selected to be R max ¼ 20 d, the real space resolution was Dr ¼ 10 À3 d and the Fourier-space convergence criterion in terms of the indirect correlation function read as jc n ðkÞ À c nÀ1 ðkÞj < 10 À5 8k. The MD simulations were performed with the RUMD software as described in Sec. IV A.
The IEMHNC and MD radial distribution functions for the biYOCP system have been compared in Fig. 4 for six different combinations of external parameters (r, l) and phase variables (C, j). Despite the lack of adjustable parameters in the IEMHNC approach, the agreement between the theory and MD simulations is excellent with the resulting g(r) being almost indistinguishable especially within

ARTICLE
scitation.org/journal/php the first coordination cell. For a more tangible comparison, the correlation void, i.e., the distance where g(r) ¼ 1/2, and the magnitude of the first maximum were extracted. For all the state points investigated, the IEMHNC deviations from MD never exceeded 1% for the correlation void, 2% for the first maximum magnitude. The results indicate that the IEMHNC approach remains remarkably accurate also for biYOCP liquids and suggest that the IEMHNC structural and thermodynamic properties can be considered as nearly exact.

C. The Rosenfeld-Tarazona decomposition of the internal energy
The Rosenfeld-Tarazona decomposition of the reduced excess internal energy for biYOCP liquids is naturally We should derive closed-form expressions for the static component u biYOCP st ðC; j; r; lÞ and the thermal component u biYOCP th ðC; j; r; lÞ. As per usual, for the static part we need to consider the unitary packing fraction limit of the PY approximation, whereas for the thermal part, we will consider the result of the IEMHNC approximation.
We begin with the static component. For the biYOCP interaction potential, the infinite-coupling limit can be expressed as Owing to the presence of the PY radial distribution function for hard spheres at g ¼ 1 packing, there is an inherent additivity in the static component with each integral representing an independent YOCP contribution. Therefore, the static part for the biYOCP is simply given by u biYOCP st ðC; j; r; lÞ ¼ Cð1 À rÞ jðj þ 1Þ ðj þ 1Þ þ ðj À 1Þe 2j þ Cr ljðlj þ 1Þ ðlj þ 1Þ þ ðlj À 1Þe 2lj : We continue with the thermal component. On the basis of Rosenfeld's theoretical and numerical analysis for the YOCP, 25 it is expected that, for C=C biYOCP m ! 0:1, the thermal part will be described by a relation of the form u biYOCP th ðC; j; r; lÞ ¼ dðr; lÞ C C biYOCP m ðj; r; lÞ with d(r, l) a function of the external potential parameters. The IEMHNC approach will be utilized to quantify the accuracy of this relation and to determine the unknown function d(r, l). A systematic parametric study was performed, where the IEMHNC approximation was numerically solved over the whole range of physically meaningful external potential parameters (r, l), within the dense liquid phase 0:1C biYOCP m ðj; r; lÞ C C biYOCP m ðj; r; lÞ and for typical normalized screening parameters 0 j 4.0. In particular, we probed five equally stepped r values ranging from 0.05 to 0.45 and eleven l values equally stepped in the 0.05-0.95 range including the Coulomb limit l ¼ 0. For each of the 55 (r, l) combinations, we studied eight equally stepped screening parameters ranging from 0.5 to 4.0 and ten equally stepped coupling parameters. In total, 4400 (r, l, C, j) combinations were investigated, which are summarized in Table I.
For each case, the IEMHNC radial distribution function was numerically computed and the reduced excess internal energy was calculated from u biYOCP ex ðC; j; r; lÞ ¼ ð3=2ÞC Ð 1 0 x½re Àjx þ ð1 ÀrÞe Àljx g ðx; C; j; r; lÞdx, where x ¼ r/d. Afterward, the static part was calculated from Eq. (23) and the thermal part was obtained from the decomposition u biYOCP th ðC; j; r; lÞ ¼ u biYOCP ex ðC; j; r; lÞ Àu biYOCP st ðC; j; r; lÞ. Then, a two-step least-square fitting procedure was followed; (1) for each fixed (r i , l i ), the corresponding d value was determined by least-square fitting the u biYOCP th ðC; j; r i ; l i Þ values obtained for all (C, j) combinations to Eq. (24). This was repeated for all (r, l) combinations and a d(r i , l i ) dataset was obtained. (2) The function d(r, l) was eventually extracted with a two-dimensional least-squares fit to the d(r i , l i ) dataset.
The results revealed that Rosenfeld's expression u th / (C/C m ) 2/5 provides an excellent description of the thermal internal energy component of dense biYOCP liquids regardless of the value of the external parameters (r, l). To be more quantitative, for any of the 55 (r i , l i ) combinations, the mean absolute relative deviations between the IEMHNC dataset (consisting of 80 points) and the Eq. (24) expression with a least-square fitted prefactor d(r i , l i ) never exceeded 2.0%. In fact, the highest deviations were recorded for l ¼ 0.00, i.e., the Yukawa-plus-Coulomb special case, which is an expected result in view of the well documented observation that the 2/5 exponent in Rosenfeld's expression becomes less accurate for Coulomb potentials. 52,54,62 The collapse of the IEMHNC results to the analytical curve has been illustrated in Fig. 5 for two (r, l) combinations.
Concerning the d(r, l) function, no systematic dependencies on either r or l were observed. As illustrated in Fig. 6, the least-square fitted d prefactor remained nearly constant as r, l vary. The largest deviations from the average value were attained for l ¼ 0.00. Excluding this case, we obtained d(r, l) ¼ 3.128 for the best-fit solution with a merely 0.23% mean relative error between the constant value and the discrete points d(r i , l i ). This result is reminiscent of previous results obtained for the YOCP case (r ¼ 0) which predicted d to lie in the range 3.0 d 3.2. 25,54,63 Furthermore, the result is also accurate for l ¼ 0.00, for which it exhibits only a 0.75% mean relative error This general result is valid for arbitrary external parameters r < 0.5, l < 1. It should be noted that the j, r, l dependencies of the thermal component stem exclusively from the melting line C biYOCP m ðj; r; lÞ. Excluding the Coulomb limit l ¼ 0, for the 4000 (r, l 6 ¼ 0, C, j) combinations investigated, the mean relative deviation between Eq. (25) and the IEMHNC approach result is 1.4% whereas the maximum relative deviation is 5.7%.
Combining all the above, the reduced excess internal energy for dense repulsive biYOCP liquids is described by Eqs. (19), (21), (23), and (25). Excluding the Coulomb limit l ¼ 0, the mean relative deviation between Eqs. (19), (21), (23), and (25) and the IEMHNC approximation result is 0.1% whereas the maximum relative deviation is 3.3% for the 4000 (r, l 6 ¼ 0, C, j) combinations studied. Higher deviations are generally observed for large values of j and low values of C/C m , see also Ref. 63. Given the four parameters involved in the biYOCP interaction, such an accurate but compact representation would have been impossible had we not taken advantage of the known YOCP representations.

V. DISCUSSION AND CONCLUSIONS
Repulsive bi-Yukawa pair-interactions have been theoretically predicted to adequately describe dust-dust interactions in complex plasmas that are engineered in weakly ionized gas discharges. [12][13][14][15][16] In isotropic complex plasmas, the short-range characteristic length arises from the polarization of the plasma background and is close to the linear Debye-H€ uckel screening length, whereas the long-range characteristic length arises from the competition between plasma ionization and plasma absorption or recombination. Complex plasmas have served as the main motivation for the present investigation, but we should point out that there is only indirect evidence (obtained in fluid demixing experiments 64 ) which support such an interaction model. Repulsive bi-Yukawa pair-interactions have also been utilized as effective ion-ion potentials for the regime of warm dense matter (WDM), 65,66 where strongly coupled classical ions coexist with degenerate electrons. 67 In WDM, the long-range characteristic length arises from the polarization of the free electron cloud and is very close to the linear Thomas-Fermi screening length, whereas the short-range characteristic length arises from the overlapping bound electron wavefunctions and leads to increased ion repulsion. 65,66,68 Such effective potentials have been extracted from density functional molecular dynamics simulations considering all electrons 68 and have also been indirectly confirmed by experiments. 69,70 Overall, the paradigm of bi-Yukawa one-component plasmas can be utilized for the study of the properties of strongly coupled systems whose interparticle interactions are characterized by two fundamental length scales. Therefore, it is entirely plausible that such effective potentials emerge in a broader range of classical soft matter systems and degenerate subatomic systems than currently suspected. Theoretical investigations of dynamic, structural, and thermodynamic properties of the dense biYOCP can greatly benefit from the comprehensive understanding of the strongly coupled YOCP. In particular,   24), as a function of the external parameters r < 0.5 and l < 1. The 55 discrete points cluster around a constant value, which is best represented by d ¼ 3.128 (red solid line).

ARTICLE
scitation.org/journal/php the approximate isomorph invariance of both model systems (together with the ensuing additivity theorems) provides a general framework that allows the straightforward extension of successful theoretical approaches and even of accurate semiempirical expressions from the single to the double Yukawa potential.
Here, we focused on structural and mainly on thermodynamic aspects of dense repulsive biYOCP liquids. An accurate description for the isomorph family of phase-diagram curves has been obtained from molecular dynamics simulations, with the analytical expression constructed on the basis of approximate additivity theorems. The isomorph-based empirically modified hypernetted-chain approach, based on the isomorph invariance of the bridge functions and known to be highly accurate for the YOCP, has been extended to the biYOCP and the resulting structural properties have revealed an excellent agreement with the results of computer simulations. In spite of the involvement of two state variables and two external parameters in the pair-interaction potential, this allowed us to obtain a simple and accurate closed-form expression for the excess internal energy of dense biYOCP liquids by capitalizing on a compact parameterization procedure that was yet again initially established for the YOCP.