# Thermodynamics of Yukawa fluids near the one-component-plasma limit

## Abstract

Thermodynamics of weakly screened (near the one-component-plasma limit) Yukawa fluids in two and three dimensions is analyzed in detail. It is shown that the thermal component of the excess internal energy of these fluids, when expressed in terms of the properly normalized coupling strength, exhibits the scaling pertinent to the corresponding one-component-plasma limit (the scalings differ considerably between the two- and three-dimensional situations). This provides us with a simple and accurate practical tool to estimate thermodynamic properties of weakly screened Yukawa fluids. Particular attention is paid to the two-dimensional fluids, for which several important thermodynamic quantities are calculated to illustrate the application of the approach.

###### pacs:

52.27.Lw, 52.25.Kn, 05.70.Ce, 64.30.-t## I Introduction

Thermodynamics of Yukawa systems (particles interacting via the pairwise Yukawa repulsive potential) is of considerable interest, in particular in the context of conventional plasmas, dusty (complex) plasmas, and colloidal suspensions. For small point-like particles the Yukawa potential (also known as screened Coulomb or Debye-Hückel potential) reads

(1) |

where is the energy scale, is the screening length, and is the distance between a pair of particles. For charged particles immersed in a neutralizing medium we have , where is the particle charge, and the screening length is equal to the corresponding Debye radius. Such system are fully characterized by two dimensionless parameters. The first is the coupling parameter, , where is the characteristic interparticle separation [in this paper we take in three dimensions and in two dimensions, where is the particle density] and is the temperature (in energy units). The second is the screening parameter, defined as .

Statical and dynamical properties of Yukawa systems in both two- and three-dimensional configurations have been extensively investigated using various theoretical and simulation techniques. Thermodynamic properties have also received considerable attention. The accurate results from Monte Carlo (MC) and molecular dynamics (MD) numerical simulations for some thermodynamic functions (in particular, for internal energy and compressibility) have been tabulated for a wide (but discrete) range of state variables and Meijer (); Tejero (); Farouki1994 (); Hamaguchi (); CG (); Totsuji04 (). Recently, simple and accurate expressions to evaluate thermodynamic properties of three-dimensional (3D) Yukawa fluids PractUP () and solids Solid () have been proposed, which are very convenient for practical applications. These expressions are based on the Rosenfeld-Tarazona (RT) scaling Rosenfeld1998 () of the excess internal energy in the fluid state (and harmonic approximation in the solid state) and deliver impressively high accuracy in a wide range of the phase diagram of 3D Yukawa systems.

The focus of the present study is on the practical estimation of the thermodynamics of Yukawa fluids near the one-component-plasma (OCP) limit. The OCP corresponds to the limit for Yukawa systems. We show that for the OCP scaling of the thermal component of the excess internal energy of 3D Yukawa fluids is superior to that of the RT scaling. This allows us to somewhat improve the accuracy of our previous practical approach to the thermodynamics of 3D Yukawa fluids PractUP () in the case of weak screening. More importantly, we show that the thermal component of the excess internal energy exhibits a quasi-universal scaling also in two-dimensional systems of soft repulsive particles (which is, however, quite different compared to the RT or OCP scaling in 3D case). Based on the 2D OCP limit we then produce a simple analytical expression for the excess energy of weakly screened 2D Yukawa systems. The accuracy of this expression is verified by comparison with the results from MD simulations. We demonstrate how this expression can be used to evaluate the specific heat, free energy, compressibility factor, and isothermal compressibility modulus of 2D Yukawa fluids in practical situations.

## Ii Yukawa fluids in three dimensions

The reduced excess energy (excess energy per particle in units of the system temperature) can be conveniently divided into static and thermal components,

(2) |

The static contribution corresponds to the value of internal energy when the particles are frozen in some regular structure (e.g. crystalline lattice) and the thermal corrections arise due the deviations from these fixed position due to the particle thermal motion. In the strongly coupled classical OCP, the simplest estimate of the static component can be obtained using the ion sphere model (ISM) Baus (); Ichimaru (). In this model each particle is restricted to the spherical cell of radius , filled with the neutralizing background. Each cell is quasineutral, the cells do not overlap, and, therefore, is simply equal to the electrostatic energy of such cell. The latter can be readily calculated, yielding Ichimaru () (the energy is negative due to particle-background interactions). In fact, this can be mathematically proven to be the exact lower bound on the excess energy of the OCP Lieb (). The ISM approach provides a rather good estimate of the actual OCP excess energy at strong coupling, in both fluid and solid phases (for example, the Madelung constant of the OCP forming a body-centered-cubic lattice is ). Another exact lower bound on the OCP excess energy is provided by the Debye-Hückel (DH) construction Mermin (). This bound is a reasonable measure of the actual energy at weak coupling. Recently a physically motivated interpolation between the DH and ISM limits has been discussed in Ref. Hybrid3D ().

When the value of the static component of the excess energy is specified (ISM approach in our case), the thermal component is also well defined. Based on the accurate Monte-Carlo (MC) simulation results from Ref. Caillol99 (), a simple two-term expression has been proposed OCP2014 ()

(3) |

This fit, along with the MC numerical data, is plotted in Fig. 1. The exponent (or close to that) in Eq. (3) provides particularly good agreement with the numerical data. This fact has been documented in a number of previous studies Caillol99 (); OCP2014 (); Stringfellow (); Farouki1994 (); Dubin1999 ().

For Yukawa systems, the reduced excess energy can also be decomposed into the static and thermal components [Eq. (2)], with the static component evaluated using the generalization of the ISM to these systems (see, in particular, Ref. ISM () for the details). Note that, in contrast to the OCP limit where the neutralizing background must be present to keep thermodynamic functions finite and well defined, for model Yukawa systems two realizations are possible. The first corresponds to the conventional system of charged particles immersed into neutralizing medium, which provides exponential screening of the particle-particle electrical interactions and keeps the overall charge neutrality (this situation can to some extent mimic real colloidal suspensions and complex plasmas IvlevBook ()). The second is the imaginary one-component system of particles interacting via the Yukawa potential without any neutralizing medium, which we will refer to as the single component Yukawa system (SCYS), and which is merely used as a useful model in condensed matter research. The difference between the excess energies of these two systems is related to the interactions between the particles and the neutralizing medium. The difference can be evaluated as ISM (); Ham94 ()

(4) |

where the first term represents the energy of the medium that, on average, neutralizes the charge of the particles (it exactly compensates the energy of the particle-particle interactions when the particles are at complete disorder). The second term gives the energy of the sheath around each particle. The static component of the internal energy in the ISM approximation is ISM ()

(5) |

In view of Eq. (4), the first term of the above equation describes the static energy of the SCYS. The coefficient of proportionality, , can be referred to as the fluid Madelung constant Rosenfeld2000 (). For SCYS we have therefore

(6) |

Note, that the same result [Eq. (6)] can be obtained using the Percus-Yevick (PY) radial distribution function of hard spheres in the limit , where is the hard sphere packing fraction Rosenfeld2000 ().

To the best of our knowledge, Rosenfeld and Tarazona were first to discuss a quasi-universal behavior of the thermal component of the internal energy () for a wide class of fluids formed by soft repulsive particle (in particular, they focused on inverse-power-law and Yukawa potentials) in three dimensions Rosenfeld1998 (); Rosenfeld2000 (). The accuracy of this scaling for various model systems has been recently investigated in extensive numerical simulations Ingebrigtsen (). It has been observed that the RT scaling is particulary useful for liquids with strong correlations between equilibrium fluctuations of virial and potential energy Ingebrigtsen (), which are referred to as Rosklide-simple or just Rosklide systems BacherNature (). OCP and strongly coupled Yukawa systems belong to this class and, thus, RT scaling can be quite useful to describe thermodynamic properties of these systems.

Previously, a variant of the RT scaling has been successfully used to obtain practical expressions for the internal energy and pressure of 3D Yukawa fluids with across coupling regimes PractUP (). The expression for the thermal component of the excess energy suggested in this work is

(7) |

where is the coupling parameter at the fluid-solid phase transition. More recently, it has been also pointed out that for Yukawa fluids sufficiently close to freezing an even simpler scaling, , is more appropriate Solid ().

The first question to be addressed now is which of the scalings, OCP [Eq. (3)] or RT [Eq. (7)] is more appropriate in the OCP limit. Figure 1 shows the comparison of these scaling with the exact MC simulation data from Ref. Caillol99 (). To plot the RT scaling we have used for the OCP in three dimensions, as reported in Ref. Hamaguchi96 () (the true for the OCP can be slightly higher, around OCP2014 (); Dubin1999 (); we use here for the sake of consistency with the Yukawa case, see below). Figure 1 demonstrates that the RT scaling describes fairly well the numerical data, but the OCP scaling is somewhat more accurate, especially in the regime . This is absolutely not surprising, since the OCP scaling is nothing but the best simple fit to the numerical data shown in Fig. 1.

The second relevant question is which of the scalings is more relevant for Yukawa systems near the OCP limit (i.e., at weak screening). To answer this, we have used the numerical data for the excess energy of Yukawa fluids tabulated in Ref. CG (), subtracted the static contribution [Eq. (5)], and plotted the remaining thermal component as a function of in Fig. 2. To do that we used the fit for , proposed in Ref. Hamaguchi96 ():

(8) |

which is applicable for . The corresponding numerical data for are shown in Fig. 2 along with the OCP and RT scalings. We see that the OCP scaling is a better approximation for . The deviations between these two scalings are particularly observable for intermediate coupling , as indicated in Fig. 2. The observation that the OCP scaling provides better description for the thermal component of the excess energy of weakly screened () Yukawa systems allows us to improve the previous practical approach to the thermodynamics of Yukawa fluids in this regime. The modifications to the expressions derived in Ref. PractUP () are straightforward and we do not discuss this further in detail, except for a relevant comment below.

It is clear that the thermodynamic functions depending on the total excess energy (like e.g. pressure and compressibility modulus), will be only slightly affected by the differences between the OCP and RT scalings of the thermal energy, since the static energy component dominates for soft repulsive potentials. However, other thermodynamic functions can be more sensitive to these differences. For example, the reduced specific heat at constant volume is determined by the reduced thermal energy alone,

(9) |

In Figure 3 we plot the available numerical data for for the OCP, along with the analytical expressions based on the OCP and RT scalings. The scattering of the numerical data points does not allow us to clearly discriminate between the OCP and RT scalings. Nevertheless, based on the comparison shown in Figs. 1 and 2, it is not unreasonable to expect that the OCP scaling is more appropriate for very soft repulsive potentials, which implies, for instance, , instead of (according to the RT scaling) in this regime.

## Iii Yukawa fluids in two dimensions

As in the previous section, we first consider the OCP limit. In two-dimensions, two different systems are actually referred to as the OCP. The first is characterized by logarithmic interactions, , where is an arbitrary length scale Totsuji79 (); Caillol82 (); Leeuw82 (). The logarithmic potential emerges from the solution of the two-dimensional Poisson equation around the central point-like particle Totsuji79 (). In this case the system is conventionally characterized by the coupling parameter , which does not depend on the particle density. In the second system, the interaction is of conventional Coulomb-type, , but the interacting particles are confined to the 2D plane. In this second case, the coupling parameter is defined similarly to the OCP (or Yukawa) system in 3D, , with being the 2D analog of the Wigner-Seitz radius (defined in the Introduction).

The numerical results for the internal energy of both OCP systems are available in the literature Caillol82 (); Leeuw82 (); Totsuji78 (); Gann79 (). The first interesting point to address, is wether the scaling of the thermal component of the excess energy in two-dimensions is also quasi-universal for soft repulsive systems (like the two OCPs considered here), similarly to the 3D case. For the static component of the energy in the 2D case we use the triangular lattice sums (Madelung energies), which are for the logarithmic and for the Coulomb potential, respectively Caillol82 (); Gann79 (). (The ion disc model, an analog of the ISM in 3D, can be constructed for logarithmic interactions Caillol82 (), but we are not aware of any such construction for the Coulomb interaction in 2D). Then, subtracting the static component from the full excess energy (when necessary) we can obtain the thermal energy component. The resulting dependence of on the reduced coupling parameter is shown in Figure 4. To produce this figure we assumed for the logarithmic interaction OCP2014 (); Caillol82 (); Leeuw82 () and for the Coulomb interaction Hartmann05 (). The closeness of values for these two systems is likely a coincidence, since the physical meaning of the coupling parameters is quite different for the logarithmic and Coulomb interactions.

The numerical data for the thermal component of the excess energy for 2D OCP systems shown in Fig. 4 tend to collapse on a single quasi-universal curve, which is analogous to the OCP scaling in three dimensions. Some scattering of the data points is present, but no clear systematic trend is observed, indicating that this can simply reflect the level of accuracy of the simulation results. The dependence of on has a logarithmic character (in contrast to the power-law scaling in 3D) and we found convenient to approximate the data by the following functional dependence:

(10) |

with the coefficients and . Also shown in Fig. 4 are the numerical data for the thermal energy of the fluid with inverse power law (IPL) interaction, , with Allen83 (). To produce this plot we have used the notation , and took at the fluid-solid phase transition (in fact, this is the fluid-hexatic transition) Kapfer (). We see that the universal scaling only holds for sufficiently soft interactions, which is not the case for the IPL case (the thermal energy is systematically lower). However, the fit of the form (10) can still be reasonably applied, see the dashed curve in Fig. 4 (the fit yields and for the IPL).

Now, by analogy with approach discussed in the previous section, we apply the observed soft-repulsion scaling of to evaluate the excess energy of weakly screened Yukawa systems in 2D geometry. The static component of the excess energy is associated with the triangular lattice sum (Madelung energy), which for the Yukawa interaction can be approximated by Totsuji04 ()

(11) |

Note that this static energy accounts for particle-particle correlations and the energy of the neutralizing medium (this is why the energy is negative). The energy of the neutralizing medium (plasma) compensates the energy coming from particle-particle interactions when no correlations are present. It can be therefore calculated using the virial energy equation with the radial distribution function , that is

(12) |

where is the Yukawa potential. The term is sometimes referred to as the positive Hartree part, while the total energy is then called the correlational part Hartmann05 (). Obviously, the energy of the SCYS in two dimensions can be obtained as a sum of the Hartree and the correlational parts. (The total energy may also include the contribution from particle-sheath interactions, similarly to the 3D case. However, this term does not affect main thermodynamic quantities of interest and we drop it from further consideration).

The dependence for two dimensional Yukawa systems has been approximated in Ref. Hartmann05 () by the following fit

(13) |

This fit underestimates by the transitional coupling parameter of the OCP, but describes relatively well the melting points found from the bond angular correlation analysis (see Fig. 6 of Ref. Hartmann05 ()) up to .

The combination of Eqs. (10) - (13) allows us to calculate the energy of strongly coupled Yukawa fluids in 2D. We compare our results with the energy data tabulated in Ref. Totsuji04 (), which are in good agreement with the subsequent simulations in Ref. Hartmann05 (). The comparison is shown in Figure 5. The agreement is rather good, the maximum relative deviation from the numerical results does not exceed few percent and is much smaller than that for the majority of the data points.

MD | Theory | MD | Theory | MD | Theory | |
---|---|---|---|---|---|---|

0.5 | 128.93 | 144.90 | 152.89 | |||

1.0 | 43.59 | 48.90 | 51.55 | |||

1.5 | 21.43 | 23.97 | 25.24 | |||

2.0 | 12.99 | 14.48 | 15.22 |

The data tabulated in Ref. Totsuji04 () do not cover the regime of very strongly coupled fluids, near the fluid-solid phase transition. We, therefore, performed additional simulations for 2D Yukawa fluids near freezing. Simulations have been performed on graphics processing unit (NVIDIA GTX 960) using the hoomd-blue software HB1 (); HB2 (). We used Yukawa particles in a rectangular box with the periodic boundary conditions. The cutoff radius for the potential has been chosen . Simulations have been performed in the canonical ensemble () with the Langevin thermostat at a temperature corresponding to the desired target coupling parameter . We first equilibrated the system for one million time steps and measured and averaged the excess energy for another million of steps. The results of this simulation along with their comparison with the theoretical predictions is summarized in Table 1. The agreement between the theory and numerical experiment is excellent for . In this regime the deviations between the theoretical and MD results are all within the simulation uncertainty. For , the theory systematically overestimates the excess energy obtained in simulations, but the relative difference is below . This qualitative trend is not unexpected. As the potential softness decreases, the OCP scaling tends to overestimate the actual thermal component of the excess energy. An example, corresponding to IPL potential, is shown in Fig. 4.

Summarizing this Section, we have proposed simple and reliable practical method to estimate thermodynamics of weakly screened Yukawa fluids in two dimensions. The accuracy of this method is very good, at least for . In the next sections we provide an exemplary calculation of several thermodynamic functions.

## Iv Thermodynamic functions of two-dimensional Yukawa fluids

We use the conventional thermodynamic identities LL () to calculate several useful quantities for Yukawa fluids in two dimensions. In doing so it is convenient to operate with the Yukawa system phase state variables and . If the density and temperature of the neutralizing medium are kept fixed, we have and . This implies

The practical model expression for can then be employed to calculate various thermodynamic quantities.

The first example deals with the reduced isochoric heat capacity, which is related to the excess energy via . Since the static contribution is linear in , is completely determined by the thermal component of the excess energy, see Eq. (9). We have calculated the dependence of on for several values of () and plotted the resulting curves in Fig. 6.

The second quantity we consider is the reduced excess free energy of Yukawa fluids. It is related to the excess energy by a standard integration,

(14) |

An exact result would be obtained if the dependence is known in the entire range of coupling, starting from . This is not the case here, because we have proposed an approximation for the sufficiently strong coupling regime (say ). The conventional procedure can be applied to start the integration from and add the weak coupling correction . The latter can be calculated using the analytical results available in the literature (see e.g. Totsuji04 ()). However, this contribution is vanishingly small at and can be omitted. Since both and are linear in in the weak coupling limit, we will not introduce serious errors by starting the integration in Eq. (14) from . The obtained results for the reduced free energy are shown in Fig. 7. These compare very favorably with the data presented in Fig. 3 of Ref. Totsuji04 () (including the regime shown in the inset of this figure).

Finally, we calculate the compressibility factor (reduced pressure),

(15) |

and the isothermal compressibility modulus,

(16) |

These quantities can be of considerable interest in the context of wave phenomena in Yukawa fluids DAV (). The calculated and , as functions of , are plotted in Fig. 8 for several values of the screening parameter . Our present calculation for the compressibility factor compares very favorably with the MD simulation results presented in Fig. 4 of Ref. Totsuji04 () (again, including the regime shown in the inset of this figure). Thus, simple and accurate practical tool to estimate thermodynamic properties of 2D Yukawa fluids is now available.

## V Conclusion

The Rosenfeld-Tarazona (RT) scaling Rosenfeld1998 (); Rosenfeld2000 () of the thermal component of the excess energy of simple fluids with soft repulsive interactions is a useful tool to estimate their thermodynamic properties. Previously, the approach based on the RT scaling arguments was applied to three-dimensional Yukawa fluids and very good accuracy was documented PractUP (); Solid (). In this paper we have shown that there is still room for some slight improvements in the regime of weak screening. In this regime, the one-component-plasma-like scaling of the thermal energy is superior to the RT scaling and thus even higher accuracy in estimating thermodynamic properties can be reached. Moreover, we have documented the existence of similar scaling (quasi-universal behavior) of the thermal energy for two-dimensional fluids with soft repulsive interactions (the functional form of the scaling varies significantly between the two- and three-dimensional cases). Based on that the total excess energy of weakly screened two-dimensional Yukawa fluids has been evaluated and good agreement with both the existing and new results from direct molecular dynamics simulations has been observed. To further demonstrate the application of this approach, various important thermodynamic quantities of the two-dimensional Yukawa fluids have been evaluated in a wide range of coupling.

###### Acknowledgements.

At a late stage of this work SAK has been supported by the A*MIDEX grant (Nr. ANR-11-IDEX-0001-02) funded by the French Government “Investissements d’Avenir” program.### Footnotes

- Also at Joint Institute for High Temperatures, Russian Academy of Sciences, Moscow, Russia

### References

- E. J. Meijer and D. Frenkel, J. Chem. Phys. 94, 2269 (1991).
- C. F. Tejero, J. F. Lutsko, J. L. Colot, and M. Baus, Phys. Rev. A 46, 3373 (1992).
- R. T. Farouki and S. Hamaguchi, J. Chem. Phys. 101, 9885 (1994).
- S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, Phys. Rev. E 56, 4671 (1997).
- J. M. Caillol and D. Gilles, J. Stat. Phys. 100, 933 (2000).
- H. Totsuji, M. S. Liman, C. Totsuji, and K. Tsuruta, Phys. Rev. E 70, 016405 (2004).
- S. A. Khrapak and H. Thomas, Phys. Rev. E 91, 023108 (2015).
- S. A. Khrapak, N. P. Kryuchkov, S. O. Yurchenko, and H. M. Thomas, J. Chem. Phys. (2015, in press).
- Y. Rosenfeld and P. Tarazona, Mol. Phys. 95, 141 (1998).
- M. Baus and J. P. Hansen, Phys. Rep. 59, 1 (1980).
- S. Ichimaru, Rev. Mod. Phys. 54, 1017 (1982).
- E. H. Lieb and H. Narnhofer, J. Stat. Phys. 12, 291 (1975).
- N. D. Mermin, Phys. Rev. 171, 272 (1968).
- S. A. Khrapak and A. G. Khrapak, Phys. Plasmas 22, 044504 (2015).
- J. M. Caillol, J. Chem. Phys. 111, 6538 (1999).
- S. A. Khrapak and A. G. Khrapak, Phys. Plasmas 21, 104505 (2014).
- G. S. Stringfellow, H. E. DeWitt, and W. L. Slattery, Phys. Rev. A 41, 1105 (1990).
- D. H. E. Dubin and T. M. O’Neil, Rev. Mod. Phys. 71, 87 (1999).
- S. A. Khrapak, A. G. Khrapak, A. V. Ivlev, and H. M. Thomas, Phys. Plasmas 21, 123705 (2014).
- A. Ivlev, H. Löwen, G. Morfill, and C. P. Royall, Complex Plasmas and Colloidal Dispersions: Particle-resolved Studies of Classical Liquids and Solids (World Scientific, Singapore, 2012).
- S. Hamaguchi and R. T. Farouki, J. Chem. Phys. 101, 9876 (1994).
- Y. Rosenfeld, Phys. Rev. E 62, 7524 (2000).
- T. S. Ingebrigtsen, A. A. Veldhorst, T. B. Schrøder, and J. C. Dyre, J. Chem. Phys. 139, 171101 (2013).
- A. K. Batcher, T. B. Schrøder, and J. C. Dyre, Nature Communications 5, 5424 (2014).
- S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, J. Chem. Phys. 105, 7641 (1996).
- S. G. Brush, H. L. Sahlin, and E. Teller, J. Chem. Phys. 45, 2102 (1966).
- J. P. Hansen, Phys. Rev. A 8, 3096 (1973).
- H. Totsuji, Phys. Rev. A 19, 2433 (1979).
- J. M. Caillol, D. Levesque, J. J. Weis and J. P. Hansen, J. Stat. Phys. 28, 325 (1982).
- S. W. de Leeuw and J. W. Perram, Physica A 113, 546 (1982).
- H. Totsuji, Phys. Rev. A 17, 399 (1978).
- R. C. Gann, S. Chakravarty and G. V. Chester, Phys. Rev. B 20, 326 (1979).
- P. Hartmann, G. J. Kalman, Z. Donko, and K. Kutasi, Phys. Rev. E 72, 026409 (2005).
- M. P. Allen, D. Frenkel, W. Gignac, and J. P. McTague, J. Chem. Phys. 79, 4206 (1983).
- S. C. Kapfer and W. Krauth, Phys. Rev. Lett. 114, 035702 (2015).
- HOOMD-blue web page: http://codeblue.umich.edu/hoomd-blue
- J. A. Anderson, C. D. Lorenz, and A. Travesset, J. Comput. Phys. 227, 5342 (2008).
- L. D. Landau and E. M. Lifshitz, Statistical Physics (Elsevier, 2005).
- S. A. Khrapak and H. M. Thomas, Phys. Rev. E 91, 033110 (2015).