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

TL;DR

This paper explores the thermodynamics and structure of dense bi-Yukawa plasmas, demonstrating isomorph invariance, deriving analytical expressions, and validating results with simulations, advancing understanding of complex plasma interactions.

Contribution

It introduces an isomorph-based approach to analyze bi-Yukawa fluids, including analytical expressions and empirical models validated by simulations.

Findings

Strong virial-potential energy correlations confirmed

Analytical isomorph family of curves derived

Excellent agreement between theory and simulations

Abstract

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.

Tables

Table 1. Table 1: Parametric study for the determination of δ ​ ( σ , μ ) 𝛿 𝜎 𝜇 \delta(\sigma,\mu) from the IEMHNC approach. For each parameter; a min subscript 𝑎 min a_{\mathrm{min}} denotes the lowest value studied, a max subscript 𝑎 max a_{\mathrm{max}} the highest value probed, N a subscript 𝑁 𝑎 N_{a} the distinct number of values investigated, Δ ​ a Δ 𝑎 \Delta{a} the constant step between successive values. Note that the Yukawa plus Coulomb limit corresponding to μ = 0 𝜇 0 \mu=0 was also studied.

Parameter	$a_{\min }$	$a_{\max }$	$N_a$	$\Delta a$
$\sigma$	0.05	0.45	5	0.1
$\mu$	0.05	0.95	10	0.1
$\Gamma /{\Gamma }_m$	0.1	1.0	10	0.1
$\kappa$	0.5	4.0	8	0.5

Equations

$$\beta u(x)=\frac{\Gamma }{x}\left[(1-\sigma )\exp \left(-\kappa x\right)+\sigma \exp \left(-\mu \kappa x\right)\right].$$

$$\beta u(x)=\frac{\Gamma }{x}\exp \left(-\kappa x\right),$$

$$\beta u(x)=\Gamma /x.$$

$$U({\mathbf{R}}_a)<U({\mathbf{R}}_b)⟺U(\lambda {\mathbf{R}}_a)<U(\lambda {\mathbf{R}}_b),$$

$$\exp \left[-{\beta }_{1}U({\mathbf{R}}_1)\right]\cong C_{12}\exp \left[-{\beta }_{2}U({\mathbf{R}}_2)\right].$$

$$U({\mathbf{R}}_2)\cong \frac{T_2}{T_1}U({\mathbf{R}}_1)+D_{12}.$$

$${\Gamma }_{\mathrm{iso}}^{\mathrm{YOCP}}(\Gamma ,\kappa )=\Gamma e^{-\alpha \kappa }\left[1+\alpha \kappa +\frac{(\alpha \kappa {)}^2}{2}\right]=const.,$$

$${\Gamma }_{\mathrm{m}}^{\mathrm{YOCP}}(\kappa )={\Gamma }_{\mathrm{m}}^{\mathrm{OCP}}{e}^{\alpha \kappa }{\left[1+\alpha \kappa +\frac{(\alpha \kappa {)}^2}{2}\right]}^{-1}.$$

$$h(r)=c(r)+n\int c(r^{\prime })h(|\mathbf{r}-{\mathbf{r}}^{\prime }|)d^3{r}^{\prime },$$

$$g(r)=\exp \left[-\beta u(r)+h(r)-c(r)+B(r)\right],$$

$$B_{\mathrm{YOCP}}\left(\frac{r}{d},\Gamma ,\kappa \right)=B_{\mathrm{OCP}}\left(\frac{r}{d},\frac{\Gamma }{e^{\alpha \kappa }}\left[1+\alpha \kappa +\frac{(\alpha \kappa {)}^2}{2}\right]\right).$$

$$B_{\mathrm{OCP}}(r,\Gamma )=\Gamma [-b_0(\Gamma )+c_1(\Gamma ){\left(\frac{r}{d}\right)}^4+c_2(\Gamma ){\left(\frac{r}{d}\right)}^6$$

$$+c_3(\Gamma ){\left(\frac{r}{d}\right)}^8]\exp \left[-\frac{b_1(\Gamma )}{b_0(\Gamma )}{\left(\frac{r}{d}\right)}^2\right],$$

$$U_{\mathrm{ex}}(n,T)=A(n)+B(n)T^{3/5},$$

$$u_{\mathrm{ex}}(\Gamma ,\kappa )=u_{\mathrm{st}}(\Gamma ,\kappa )+u_{\mathrm{th}}(\Gamma ,\kappa ),$$

$$u_{\mathrm{st}}^{\mathrm{YOCP}}(\Gamma ,\kappa )=6\Gamma \underset{\eta \to 1}{\lim }{\int }_0^{\infty }y{g}_{\mathrm{PY}}(y;\eta )e^{-2{\eta }^{1/3}\kappa y}dy,$$

$$u_{\mathrm{st}}^{\mathrm{YOCP}}(\Gamma ,\kappa )=\frac{\kappa (\kappa +1)\Gamma }{(\kappa +1)+(\kappa -1)e^{2\kappa }},$$

$$u_{\mathrm{th}}^{\mathrm{YOCP}}(\Gamma ,\kappa )=3.0{\left[\frac{\Gamma }{{\Gamma }_{\mathrm{m}}^{\mathrm{YOCP}}(\kappa )}\right]}^{2/5}.$$

$${\Gamma }_{\mathrm{iso}}^{\mathrm{biYOCP}}$$

$$+\sigma e^{-\mu \alpha \kappa }\left[1+\mu \alpha \kappa +\frac{(\mu \alpha \kappa {)}^2}{2}\right]\}=\mathrm{const}$$

$${\Gamma }_{\mathrm{m}}^{\mathrm{biYOCP}}(\kappa ;\sigma ,\mu )$$

$$+\frac{\sigma }{e^{\mu \alpha \kappa }}\left[1+\mu \alpha \kappa +\frac{(\mu \alpha \kappa {)}^2}{2}\right]\}.$$

$$B_{\mathrm{biYOCP}}\left(\frac{r}{d},\Gamma ,\kappa ;\sigma ,\mu \right)=B_{\mathrm{OCP}}(\frac{r}{d},\Gamma \{\frac{(1-\sigma )}{e^{\alpha \kappa }}[1+$$

$$\alpha \kappa +\frac{(\alpha \kappa {)}^2}{2}]+\frac{\sigma }{e^{\mu \alpha \kappa }}\left[1+\mu \alpha \kappa +\frac{(\mu \alpha \kappa {)}^2}{2}\right]\}).$$

$$u_{\mathrm{ex}}^{\mathrm{biYOCP}}(\Gamma ,\kappa ;\sigma ,\mu )=$$

$$u_{\mathrm{th}}^{\mathrm{biYOCP}}(\Gamma ,\kappa ;\sigma ,\mu ).$$

$$u_{\mathrm{st}}^{\mathrm{biYOCP}}=6\Gamma (\sigma -1)\underset{\eta \to 1}{\lim }{\int }_0^{\infty }y{g}_{\mathrm{PY}}(y;\eta )e^{-2{\eta }^{1/3}\kappa y}dy$$

$$+6\sigma \Gamma \underset{\eta \to 1}{\lim }{\int }_0^{\infty }y{g}_{\mathrm{PY}}(y;\eta )e^{-2{\eta }^{1/3}\kappa \mu y}dy.$$

$$u_{\mathrm{st}}^{\mathrm{biYOCP}}$$

$$+\Gamma \sigma \frac{\mu \kappa (\mu \kappa +1)}{(\mu \kappa +1)+(\mu \kappa -1)e^{2\mu \kappa }}.$$

$$u_{\mathrm{th}}^{\mathrm{biYOCP}}(\Gamma ,\kappa ;\sigma ,\mu )=\delta (\sigma ,\mu ){\left[\frac{\Gamma }{{\Gamma }_{\mathrm{m}}^{\mathrm{biYOCP}}(\kappa ;\sigma ,\mu )}\right]}^{\frac{2}{5}}$$

$$u_{\mathrm{th}}^{\mathrm{biYOCP}}(\Gamma ,\kappa ;\sigma ,\mu )=3.128{\left[\frac{\Gamma }{{\Gamma }_{\mathrm{m}}^{\mathrm{biYOCP}}(\kappa ;\sigma ,\mu )}\right]}^{\frac{2}{5}}.$$

Full text

I Introduction

Complex (or dusty) plasmas can be defined as weakly ionized plasmas that are dispersed with solid charged particulates of mesoscopic size introd1 ; introd2 . 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 introd1 ; introd2 . 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 introd2 ; introd3 ; introd4 . 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ückel or Yukawa type introd2 ; introd3 ; introd4 ; introd5 , $u(r)=(Q^{2}/r)\exp{(-r/\lambda_{\mathrm{D}})}$ where $Q$ is the dust charge and $\lambda_{\mathrm{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 introd6 and corresponds to the paradigm of the Yukawa one-component plasma (YOCP) which reduces to the classical one-component plasma (OCP) in the absence of screening ($\lambda_{\mathrm{D}}\to\infty$).

Yukawa interactions result from the linearized Boltzmann response of isotropic plasmas around a test point charge introd2 . 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 potent0 ; potent1 ; potent2 and brings forth the issue of ionization-recombination balance for dense dust clouds potent3 ; potent4 . In the latter, simple hydrodynamic models have predicted that the interaction potential acquires a double Yukawa repulsive form potent5 ; potent6 , 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 potent7 ; potent8 ; potent9 . We shall refer to such strongly-interacting systems as bi-Yukawa one-component plasmas (biYOCP). 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 gelati1 ; gelati2 ; gelati3 .

The general form of the bi-Yukawa pair-interaction potential is $u(r)=(\epsilon_{1}/r)\exp{(-r/\lambda_{1})}+(\epsilon_{2}/r)\exp{(-r/\lambda_{2})}$. We shall limit ourselves to purely repulsive interactions ($\epsilon_{1},\epsilon_{2}>0$) with a pure Coulomb short-range asymptote ($\epsilon_{1}+\epsilon_{2}=Q^{2}$) and a dominant short-range screening term ($\lambda_{2}\geq\lambda_{1},\epsilon_{2}\leq\epsilon_{1}$). We will focus exclusively on extended three-dimensional dense biYOCP liquids. Expressed in reduced coordinates, $x=r/d$ where $d=\left(4\pi{n}/3\right)^{-1/3}$ is the Wigner-Seitz radius with $n$ the particle density, the biYOCP interaction potential becomes

[TABLE]

In the above $\beta=1/(k_{\mathrm{B}}T)$, where $T$ is the particle temperature and $k_{\mathrm{B}}$ is the Boltzmann constant. The dimensionless thermodynamic state variables of the biYOCP are the coupling parameter defined by $\Gamma=\beta{Q}^{2}/d$ and the normalized screening parameter defined by $\kappa=d/\lambda_{1}$. The external potential parameters of the biYOCP are the relative strength $\sigma=\epsilon_{2}/{Q}^{2}$ and the relative screening length $\mu=\lambda_{1}/\lambda_{2}$ subject to the constraints $0\leq\sigma\leq 0.5$ and $0\leq\mu\leq 1$. We shall assume that, in contrast to the state variables $(\Gamma,\kappa)$, the external parameters $(\sigma,\mu)$ 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 $\sigma=0$

[TABLE]

and that it reduces to the OCP for $\kappa=0$

[TABLE]

In spite of the abundance of theoretical studies that predict a bi-Yukawa interaction potential for many laboratory realizations of complex plasmas potent5 ; potent6 ; potent7 ; potent8 ; potent9 , there have been few investigations of the thermodynamic, static and dynamic properties of the strongly coupled biYOCP biYuka1 ; biYuka2 ; biYuka3 . This probably stems from the fact that the two state variables and two external potential parameters form a four-dimensional $(\Gamma,\kappa;\sigma,\mu)$ 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 $(\Gamma,\kappa;\sigma,\mu)$ 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 hypernetted-chain (IEMHNC) approximation approx3 and on taking advantage of the compact parameterization offered by the Rosenfeld-Tarazona decomposition RTdeco1 in combination with the Rosenfeld scaling DecRos1 .

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 Roski01 ; Roski02 ; Roski03 ; Roski04 ; Roski05 ; Roski06 ; Roski07 ; Roski08 ; Roski09 ; Roski13 . Hence, the strong virial-potential energy correlations of the dense biYOCP are first demonstrated with the aid of molecular dynamics 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 RTdeco1 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 DecRos1 .

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 Roski01 ; Roski02 ; Roski03 ; Roski04 ; Roski05 ; Roski06 ; Roski07 ; Roski08 ; Roski09 ; Roski13 , in which the length is normalized to the mean-cubic inter-particle distance $a=n^{-1/3}$, the energy to $k_{\mathrm{B}}T$ and the time to $\tau=n^{-1/3}\sqrt{m/(k_{\mathrm{B}}T)}$, 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.

II.1 Definition of R-simple systems

We first introduce the collective N-particle position vector configuration $\boldsymbol{R}=(\boldsymbol{r}_{1},\boldsymbol{r}_{2},...,\boldsymbol{r}_{\mathrm{N}})$, the total potential energy $U(\boldsymbol{R})$ and the microscopic virial $W(\boldsymbol{R})=(1/3)\sum_{i}\boldsymbol{r}_{i}\cdot\boldsymbol{F}_{i}=-(1/3)\boldsymbol{R}\cdot\nabla{U}(\boldsymbol{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}=\langle\Delta{W}\Delta{U}\rangle/\sqrt{\langle(\Delta{W})^{2}\rangle\langle(\Delta{U})^{2}\rangle}$ where $\langle...\rangle$ denotes canonical ensemble averaging and $\Delta{A}=A-\langle{A}\rangle$ statistical fluctuations around the thermodynamic mean. The strong $W-U$ correlations are empirically delimited by the practical condition $R_{WU}\gtrsim 0.9$ Roski01 ; Roski02 ; Roski03 ; Roski04 . Such a definition enables a straightforward identification of R-simple systems with the aid of $NVT$ Molecular Dynamics (MD) simulations Roski02 ; Roski03 , 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 Roski06 ; Roski07 , according to which an R-simple system is defined by the relation

[TABLE]

where the equilibrium configurations $\boldsymbol{R}_{\mathrm{a}},\,\boldsymbol{R}_{\mathrm{b}}$ correspond to the same density (not necessarily configurations of the same temperature) with $\lambda$ 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 W-U correlations ($R_{WU}\equiv 1$) which is strictly true only for Euler homogeneous pair-interactions such as those realized in inverse-power 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 Roski06 ; Roski07 ; Roski08 ; Roski09 : (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 isomorph theory Roski09 . (3) Two arbitrary configurations $\boldsymbol{R}_{1},\boldsymbol{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.

[TABLE]

where the $C_{12}$ factor depends on the two thermodynamic state points but not on the microscopic configurations.

II.2 Direct isomorph check

Different methods have been developed for the tracing of the isomorph curves of R-simple systems via MD simulations Roski02 ; Roski05 ; Roski08 . The so-called direct isomorph check is based on the equal canonical probabilities and works as follows Roski02 . Eq.(5) can be rewritten as

[TABLE]

Since these two configurations should be identical in reduced units, we also have that $n_{1}^{1/3}\boldsymbol{R}_{1}=n_{2}^{1/3}\boldsymbol{R}_{2}$. Let us suppose a $(n_{1},T_{1})$ reference state and an isomorphic state $(n_{2},T_{2})$ of re-scaled density $n_{2}=(1/\lambda){n}_{1}$. The goal is to determine $T_{2}$. The potential energy $U(\boldsymbol{R}_{1})$ is first extracted by an equilibrium $(n_{1},T_{1})$ simulation, the configuration is then re-scaled to $\boldsymbol{R}_{2}=\lambda\boldsymbol{R}_{1}$ and the potential energy $U(\boldsymbol{R}_{2})$ is obtained. This procedure is repeated for several configurations and a $U(\boldsymbol{R}_{1}),\,U(\boldsymbol{R}_{2})$ scatter plot is obtained. As a consequence of Eq.(6), the points should be well-approximated by a straight line whose slope $\theta$ can be obtained with linear regression. This allows for the determination of $T_{2}$ from $T_{2}=\theta T_{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 $\lambda$. However, the approximate nature of isomorph invariance for non-IPL systems and thus of Eq.(6) should be taken into account. This imposes a $\lambda\lesssim 1.1$ restriction, which indicates that tracing is usually performed with density increments $(n_{1}-n_{2})/n_{1}\lesssim 10\%$ Roski08 . 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 Roski08 .

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.

III.1 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 Roski08 . 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 Roski08

[TABLE]

where $\alpha=a/d=(4\pi/3)^{1/3}$ is the ratio of the cubic mean inter-particle distance to the Wigner-Seitz radius.

This expression is equivalent to the semi-empirical relation earlier proposed for the YOCP melting line isoKhr1 ; isoKhr2 that well describes the extensive MD liquid-solid phase coexistence data available in the literature Hamagu1 ; Hamagu2

[TABLE]

In the above $\Gamma_{\mathrm{m}}^{\mathrm{OCP}}=171.8$ is the coupling parameter at the OCP melting point Hamagu1 .

III.2 Isomorph-based integral theory approach

For isotropic one-component systems, the integral equation theory of liquid structure comprises of the Ornstein-Zernike equation bookre1

[TABLE]

combined with the exact closure equation bookre1

[TABLE]

where $h(r)=g(r)-1$ is the total correlation function, $g(r)$ the radial distribution function, $c(r)$ the direct correlation function and $B(r)$ 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+\gamma(r)]-\gamma(r)$ corresponds to the Percus-Yevick (PY) approximation bookre1 , where we introduced the indirect correlation function $\gamma(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 approx3 . For strongly coupled Yukawa systems, this recently developed approach assumes that the YOCP bridge function $B_{\mathrm{YOCP}}(r,\Gamma,\kappa)$ is an isomorph invariant up to and including the OCP limit and takes advantage of the availability of simulation-extracted OCP bridge functions $B_{\mathrm{OCP}}(r,\Gamma)$. From Eq.(7) we have

[TABLE]

Accurate closed-form expressions are available for the OCP bridge function $B_{\mathrm{OCP}}(r,\Gamma)$ that have been obtained with the aid of Monte Carlo simulations and read as OCPbri1

[TABLE]

where the coefficients $b_{0}(\Gamma)$, $b_{1}(\Gamma)$, $c_{1}(\Gamma)$, $c_{2}(\Gamma)$ and $c_{3}(\Gamma)$ are described by Eqs.(21,23) of Ref.OCPbri1 . The IEMHNC approach has been validated against extensive YOCP computer simulations approx3 performing consistently better than other rigorous integral theory approaches approx4 ; approx6 ; approx7 ; approx8 as far as thermodynamic properties and structural quantities are concerned. Overall, it has been concluded that within its range of validity, $5.25\leq\Gamma_{\mathrm{iso}}^{\mathrm{YOCP}}(\Gamma,\kappa)\leq 171.8$, the IEMHNC radial distribution functions are accurate within $1.2\%$ inside the first co-ordination cell as well as that the IEMHNC excess internal energies and pressures are accurate within $0.5\%$ approx3 .

III.3 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 RTdeco1 . In standard $(n,T)$ state variables, the RT decomposition reads as

[TABLE]

which leads to the well-known $T^{-2/5}$ RT scaling for the excess isochoric heat capacity $C_{\mathrm{V}}^{\mathrm{ex}}=\left(\partial{U}_{\mathrm{ex}}/\partial{T}\right)_{V}$. The RT scaling has been demonstrated to be closely followed by R-simple liquids RTdeco2 . It is worth pointing out that the original formulation of isomorph theory proved that $C_{\mathrm{V}}^{\mathrm{ex}}$ is isomorph invariant Roski02 , but this does not apply in the more generic reformulation Roski06 and has been confirmed not to apply for all R-simple liquids RTdeco3 ; RTdeco4 .

For the YOCP system, when expressed in $(\Gamma,\kappa)$ state variables and via reduced excess internal energies ${u}_{\mathrm{ex}}=U_{\mathrm{ex}}/(Nk_{\mathrm{B}}T)$, the RT decomposition reads as RTdeco1

[TABLE]

where $u_{\mathrm{ex}}(\Gamma,\kappa)=(3/2)\int_{0}^{\infty}x^{2}[\beta{u}(x;\Gamma,\kappa)]g(x;\Gamma,\kappa)dx$ in Wigner-Seitz coordinates $x=r/d$.

The static component $u_{\mathrm{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 un-physical limit, the static correlations for repulsive pair additive systems coincide with the unitary packing fraction limit of the PY approximation for hard spheres RTdeco1 . For the YOCP we have

[TABLE]

with $\eta=(\pi/6)nr_{\mathrm{hs}}^{3}$ the packing fraction, $r_{\mathrm{hs}}$ the hard sphere radius and $y=r/r_{\mathrm{hs}}$. The integral involving the PY hard-sphere radial distribution function $g_{\mathrm{PY}}(y;\eta)$ has an analytical form RTdeco5 ; RTdeco6 that, in the limit $\eta\rightarrow 1$, yields

[TABLE]

which has a removable divergence in the OCP $\kappa\to 0$ limit $u_{\mathrm{st}}^{\mathrm{YOCP}}(\Gamma,\kappa)=(3/2)\Gamma/\kappa^{2}-9\Gamma/10+\mathcal{O}(\kappa)$ that cancels out exactly with the background contribution leading to the well-known OCP result $u_{\mathrm{st}}^{\mathrm{OCP}}=-9\Gamma/10$ RTdeco7 ; RTdeco8 . Note that Eq.(16) can also be derived within the Ion-Sphere Model (ISM) RTdeco7 ; RTdeco8 extended to the YOCP RTdeco9 ; RTdeco0 . 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_{\mathrm{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_{\mathrm{th}}\propto\Gamma^{2/5}$ RTdeco1 , while a detailed analysis of YOCP MD simulations by Rosenfeld concluded that, when $\Gamma/\Gamma_{\mathrm{m}}^{\mathrm{YOCP}}\geq 0.1$, $u_{\mathrm{th}}$ is well approximated by DecRos1

[TABLE]

We point out that the $\kappa-$dependence of the thermal component stems exclusively from the melting line $\Gamma_{\mathrm{m}}^{\mathrm{YOCP}}(\kappa)$. In view of Eqs.(7,8), the above expression suggests that, for $\Gamma/\Gamma_{\mathrm{m}}^{\mathrm{YOCP}}\geq 0.1$, the YOCP belongs to a subclass of R-simple liquids characterized by an isomorph invariant $C_{\mathrm{V}}^{\mathrm{ex}}$. Finally, we note that an empirical expression similar to Eq.(17), reading as $u_{\mathrm{th}}^{\mathrm{YOCP}}=3.2({\Gamma}/{\Gamma_{\mathrm{m}}^{\mathrm{YOCP}}})^{2/5}-0.1$, gives a highly accurate prediction for the thermal component down to $\Gamma/\Gamma_{\mathrm{m}}^{\mathrm{YOCP}}=10^{-3}$ RTdeco0 .

IV Isomorphic aspects of the biYOCP structure and thermodynamics

IV.1 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 Roski04 ; Roski10 . This additivity theorem, based on the quasi-universality of simple liquids Roski07 , clearly suggests that the dense biYOCP system is R-simple. Furthermore, Eq.(5) suggests that the isomorph curves are also additive when expressed in $(\Gamma,\kappa)$ state variables. Therefore, the biYOCP isomorph curves will be given by $\Gamma_{\mathrm{iso}}^{\mathrm{biYOCP}}(\Gamma,\kappa;\sigma,\mu)=\Gamma_{\mathrm{iso}}^{\mathrm{YOCP}}[\Gamma(1-\sigma),\kappa]+\Gamma_{\mathrm{iso}}^{\mathrm{YOCP}}(\Gamma\sigma,\mu\kappa)$. Courtesy of Eq.(7), this leads to the analytical expression

[TABLE]

Finally, invoking Rosenfeld’s additivity of melting temperatures Rosenf1 ; Rosenf2 (explained by isomorph theory Roski10 ) and re-adjusting the coupling parameter strengths, we have $1/\Gamma_{\mathrm{m}}^{\mathrm{biYOCP}}(\kappa;\sigma,\mu)=(1-\sigma)/\Gamma_{\mathrm{m}}^{\mathrm{YOCP}}(\kappa)+\sigma/\Gamma_{\mathrm{m}}^{\mathrm{YOCP}}(\mu\kappa)$. Combining with Eq.(8), we get the biYOCP melting line

[TABLE]

As emphasized in Sec.II.1, 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 Roski11 . The simulations were performed directly in reduced units with $N=8192$ particles and a time-step of $\Delta{t}=2.5\times 10^{-3}\tau$. The interaction potential was truncated at the cutoff distance $r_{\mathrm{cut}}=10a$ with the shifted-force cutoff method Roski12 . Since the Ewald decomposition is not implemented in RUMD, long-range interactions close to the unscreened Coulomb limit cannot be handled accurately Ewald00 ; Ewald01 and we confined ourselves to $\kappa\geq 1$. We also focused on the dense fluid region of the phase diagram, $0.1<\Gamma/\Gamma_{\mathrm{m}}^{\mathrm{biYOCP}}<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 figure 1 for $\sigma=0.10,\,\mu=0.15$ and two thermodynamic states, both characterized by $R_{WU}>0.99$. Overall, it is concluded that biYOCP liquids are R-simple for any $\sigma<0.5$, $\mu<1$ at least in the phase region defined by $0.1<\Gamma/\Gamma_{\mathrm{m}}^{\mathrm{biYOCP}}<1$ and $\kappa\geq 1$.

The direct isomorph check with a screening parameter increment $\Delta\kappa/\kappa=3.5\%$ was employed in order to trace isomorphs in the interval $1.0\leq\kappa\leq 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 figure 2 and detailed in the respective caption, see also Sec.II.2. The numerical isomorph curves are illustrated in figure 3 together with the analytical expression of Eq.(18). We considered three reduced coupling parameters $\Gamma/\Gamma_{\mathrm{m}}^{\mathrm{biYOCP}}=\left\{0.1,0.5,0.8\right\}$ and two $(\sigma,\,\mu)$ combinations: $\sigma=0.10$, $\mu=0.15$ in Figure 3a and $\sigma=0.20$, $\mu=0.25$ in Figure 3b.

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 Hamagu2 and from the analytical expression of Eq.(8). As expected from the longer range of the 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 $\kappa=0$ where both systems reduce to the OCP.

IV.2 Isomorph-based integral theory approach

The results of the previous subsection allow us to extend the IEMHNC approach to biYOCP systems. Taking into account the fact that the biYOCP reduces to the OCP for $\kappa\to 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 approx3 ,

[TABLE]

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\leq\Gamma_{\mathrm{iso}}^{\mathrm{biYOCP}}(\Gamma,\kappa;\sigma,\mu)\leq 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,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 approx3 . The upper range cut-off was selected to be $R_{\mathrm{max}}=20\,d$, the real space resolution was $\Delta{r}=10^{-3}\,d$ and the Fourier-space convergence criterion in terms of the indirect correlation function read as $|\gamma_{n}(k)-\gamma_{n-1}(k)|<10^{-5}\,\forall{k}$. The MD simulations were performed with the RUMD software as described in section IV.1.

The IEMHNC and MD radial distribution functions for the biYOCP system have been compared in Figure 4 for six different combinations of external parameters $(\sigma,\mu)$ and phase variables $(\Gamma,\kappa)$. 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 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.

IV.3 The Rosenfeld-Tarazona decomposition of the internal energy

The Rosenfeld-Tarazona decomposition of the reduced excess internal energy for biYOCP liquids is naturally

[TABLE]

We should derive closed-form expressions for the static component $u_{\mathrm{st}}^{\mathrm{biYOCP}}(\Gamma,\kappa;\sigma,\mu)$ and the thermal component $u_{\mathrm{th}}^{\mathrm{biYOCP}}(\Gamma,\kappa;\sigma,\mu)$. 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

[TABLE]

Owing to the presence of the PY radial distribution function for hard spheres at $\eta=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

[TABLE]

We continue with the thermal component. On the basis of Rosenfeld’s theoretical and numerical analysis for the YOCP DecRos1 , it is expected that, for $\Gamma/\Gamma_{\mathrm{m}}^{\mathrm{biYOCP}}\geq 0.1$, the thermal part will be described by a relation of the form

[TABLE]

with $\delta(\sigma,\mu)$ 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 $\delta(\sigma,\mu)$.

A systematic parametric study was performed, where the IEMHNC approximation was numerically solved over the whole range of physically meaningful external potential parameters $(\sigma,\mu)$, within the dense liquid phase $0.1\Gamma_{\mathrm{m}}^{\mathrm{biYOCP}}(\kappa;\sigma,\mu)\leq\Gamma\leq\Gamma_{\mathrm{m}}^{\mathrm{biYOCP}}(\kappa;\sigma,\mu)$ and for typical normalized screening parameters $0\leq\kappa\leq 4.0$. In particular, we probed five equally stepped $\sigma$ values ranging from $0.05$ to $0.45$ and eleven $\mu$ values equally stepped in the $0.05$ to $0.95$ range including the Coulomb limit $\mu=0$. For each of the $55$ $(\sigma,\mu)$ combinations, we studied eight equally stepped screening parameters ranging from $0.5$ to $4.0$ and ten equally stepped coupling parameters. In total, $4400$ $(\sigma,\mu,\Gamma,\kappa)$ combinations were investigated, which are summarized in Table 1.

For each case, the IEMHNC radial distribution function was numerically computed and the reduced excess internal energy was calculated from $u_{\mathrm{ex}}^{\mathrm{biYOCP}}(\Gamma,\kappa;\sigma,\mu)=(3/2)\Gamma\int_{0}^{\infty}x[\sigma{e}^{-\kappa{x}}+(1-\sigma){e}^{-\mu\kappa{x}}]g(x;\Gamma,\kappa;\sigma,\mu)dx$ where $x=r/d$. Afterwards, the static part was calculated from Eq.(23) and the thermal part was obtained from the decomposition $u_{\mathrm{th}}^{\mathrm{biYOCP}}(\Gamma,\kappa;\sigma,\mu)=u_{\mathrm{ex}}^{\mathrm{biYOCP}}(\Gamma,\kappa;\sigma,\mu)-u_{\mathrm{st}}^{\mathrm{biYOCP}}(\Gamma,\kappa;\sigma,\mu)$. Then, a two-step least-square fitting procedure was followed; (1) For each fixed $(\sigma_{i},\mu_{i})$, the corresponding $\delta$ value was determined by least-square fitting the $u_{\mathrm{th}}^{\mathrm{biYOCP}}(\Gamma,\kappa;\sigma_{i},\mu_{i})$ values obtained for all $(\Gamma,\kappa)$ combinations to Eq.(24). This was repeated for all $(\sigma,\mu)$ combinations and a $\delta(\sigma_{i},\mu_{i})$ dataset was obtained. (2) The function $\delta(\sigma,\mu)$ was eventually extracted with a two-dimensional least-squares fit to the $\delta(\sigma_{i},\mu_{i})$ dataset.

The results revealed that Rosenfeld’s expression $u_{\mathrm{th}}\propto(\Gamma/\Gamma_{\mathrm{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 $(\sigma,\mu)$. To be more quantitative, for any of the $55$ $(\sigma_{i},\mu_{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 pre-factor $\delta(\sigma_{i},\mu_{i})$ never exceeded $2.0\,\%$. In fact, the highest deviations were recorded for $\mu=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 RTdeco0 ; RTdeco8 ; fitKhra . The collapse of the IEMHNC results to the analytical curve has been illustrated in figure 5 for two $(\sigma,\mu)$ combinations.

Concerning the $\delta(\sigma,\mu)$ function, no systematic dependencies on either $\sigma$ or $\mu$ were observed. As illustrated in figure 6, the least-square fitted $\delta$ pre-factor remained nearly constant as $\sigma,\mu$ vary. The largest deviations from the average value were attained for $\mu=0.00$. Excluding this case, we obtained $\delta(\sigma,\mu)=3.128$ for the best-fit solution with a merely $0.23\%$ mean relative error between the constant value and the discrete points $\delta(\sigma_{i},\mu_{i})$. This result is reminiscent of previous results obtained for the YOCP case ($\sigma=0$) which predicted $\delta$ to lie in the range $3.0\leq\delta\leq 3.2$ DecRos1 ; RTdeco0 ; fitKhrR . Furthermore, the result is also accurate for $\mu=0.00$, for which it exhibits only a $0.75\%$ mean relative error with respect to the discrete points. Altogether, the thermal component of the reduced excess internal energy for dense biYOCP liquids is given by

[TABLE]

This general result is valid for arbitrary external parameters $\sigma<0.5$, $\mu<1$. It should be noted that the $\kappa,\sigma,\mu$ dependencies of the thermal component stem exclusively from the melting line $\Gamma_{\mathrm{m}}^{\mathrm{biYOCP}}(\kappa;\sigma,\mu)$. Excluding the Coulomb limit $\mu=0$, for the $4000$ $(\sigma,\mu\neq 0,\Gamma,\kappa)$ 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,25). Excluding the Coulomb limit $\mu=0$, the mean relative deviation between Eqs.(19,21,23,25) and the IEMHNC approximation result is $0.1\%$ whereas the maximum relative deviation is $3.3\%$ for the $4000$ $(\sigma,\mu\neq 0,\Gamma,\kappa)$ combinations studied. Higher deviations are generally observed for large values of $\kappa$ and low values of $\Gamma/\Gamma_{\mathrm{m}}$, see also Ref.fitKhrR . 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 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 potent5 ; potent6 ; potent7 ; potent8 ; potent9 . 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ückel 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 are only indirect evidence (obtained in fluid de-mixing experiments potentO ) which support such interaction model. Repulsive bi-Yukawa pair-interactions have also been utilized as effective ion-ion potentials for the regime of warm dense matter (WDM) WDMpap1 ; WDMpap2 , where strongly coupled classical ions are coexisting with degenerate electrons WDMpapg . 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 WDMpap1 ; WDMpap2 ; WDMpap3 . Such effective potentials have been extracted from density functional molecular dynamics simulations considering all electrons WDMpap3 and have also been indirectly confirmed by experiments WDMpap4 ; WDMpap5 .

Overall, the paradigm of bi-Yukawa one-component plasmas can be utilized for the study of the properties of strongly coupled systems whose inter-particle 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, 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 semi-empirical 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.

Acknowledgments

FLC and PT acknowledge the financial support of the Swedish National Space Agency. This work was also supported by the VILLUM Foundation’s Grant No. 16515 (Matter).