Simple estimation of thermodynamic properties of Yukawa systems

# Simple estimation of thermodynamic properties of Yukawa systems

S. A. Khrapak, A. G. Khrapak, A. V. Ivlev, and G. E. Morfill Max-Planck-Institut für extraterrestrische Physik, D-85741 Garching, Germany
Joint Institute for High Temperatures RAS, 125412 Moscow, Russia
July 12, 2019
###### Abstract

A simple analytical approach to estimate thermodynamic properties of model Yukawa systems is presented. The approach extends the traditional Debye-Hückel theory into the regime of moderate coupling and is able to qualitatively reproduce thermodynamics of Yukawa systems up to the fluid-solid phase transition. The simplistic equation of state (pressure equation) is derived and applied to the hydrodynamic description of the longitudinal waves in Yukawa fluids. The relevance of this study to the topic of complex (dusty) plasmas is discussed.

###### pacs:
52.27.Lw; 52.25.Kn

## I Introduction

Yukawa systems are many-particle systems characterized by the pair interaction potential of the form

 U(r)=ϵ(σ/r)exp(−r/σ), (1)

where and are the energy and length scales, and is the distance between two particles. This potential is often used to describe interactions in systems of charged particles immersed in a neutralizing medium. Two well known examples are colloidal dispersions and complex (dusty) plasmas IvlevBook (); FortovBook (). A remarkable property, explaining the significance of Yukawa systems in soft condensed matter research, is that by varying it is possible to explore the extremely broad range of interaction steepness: From extremely soft Coulomb interactions () to very hard, almost hard sphere interactions ().

Various aspects of Yukawa systems have been the subject of study in the last several decades. This includes fluid-solid and solid-solid phase transitions and the emerging phase diagram, equilibrium transport properties, wave modes, confined systems and finite clusters, mixtures, various non-equilibrium phase transitions, etc. The attention has been focused on both two and three dimensional systems. It would be almost impossible to give credits to all related original works here, the reader is referred to books IvlevBook (); FortovBook () and some review papers Lowen (); FortovPR (); MorfillRMP (); Bonitz (); Chaudhuri () instead.

Thermodynamic properties of Yukawa systems have also been extensively investigated. Molecular dynamics simulations HamaguchiJCP1994 (); Farouki (); HamaguchiJCP_1996 (); Hamaguchi () as well as integral equation theory (in the hypernetted chain approximation) Kalman2000 () have been used to calculate accurately the system energy. Since differentiations and integrations are required to obtain other thermodynamic quantities, their accurate determination remains a demanding computational task. Rather high accuracy is required in some cases. For instance, when locating fluid-solid and solid-solid phase transitions the free energies of respective phases have to be known with extreme accuracy since the smallest change in the free energy of either phase can result in a significant deviation from the actual coexistence line Hamaguchi (); DeWitt2001 (). In some other cases no such accuracy is required and it would be valuable to have simple analytical expressions instead, which allow to estimate main thermodynamic properties of the system.

The purpose of this paper is to discuss such an approach for estimating thermodynamic properties of Yukawa systems. Although not extremely accurate, it is very simple and allows to catch the essential qualitative properties of these systems.

In the following we consider an idealized model consisting of point-like charged particles in the neutralizing surrounding medium. The mobile medium is responsible for screening so that the resulting pair interaction potential between the particles has the Yukawa form (1). The main emphasize is on complex (dusty) plasmas and the relevance of the present idealized model to these systems will be discussed towards the end of the paper. In particular, we will point out that the model itself does not account for some important properties of complex plasmas. This implies that approximate analytical schemes which potentially can be extended to account for these properties are not irrelevant, although highly accurate data for an idealized model do exist.

The paper is organized as follows. In Section II we specify the model. In Section III we briefly remind the standard Debye-Hückel approximation for weakly coupled Yukawa systems. The improvement of this model, which is the main subject of this paper, is described in Section IV. A limiting case of Yukawa systems – the one-component plasma limit is briefly discussed in Section V. We then proceed with the derivation of an approximate equation of state for Yukawa systems in Section VI. Its application to the analysis of waves in strongly coupled (fluid) Yukawa systems is described in Section VII. Section VIII presents discussion and conclusion.

## Ii Model

We consider the two-component system consisting of microparticles of charge and density and neutralizing medium, characterized by the charge and density (subscript denotes “medium”). In equilibrium the system is quasineutral, so that

 Qn0−enm0=0, (2)

where the subscript denotes unperturbed quantities. (If the neutralizing medium is comprised of several species and some of them are oppositely charged, this should be taken into account: e.g., for an electron-ion medium). It is conventional to characterize such a system by two dimensionless parameters:

 Γ=Q2aT    and    κ=akm, (3)

where is the Wigner-Seitz radius, is the system temperature (in energy units), and is the inverse screening length (Debye radius) associated with the total density of the neutralizing medium. In our case , for the electron-ion medium . In principle, the particle species and surrounding medium can be characterized by different temperatures, but this is not important for the present consideration. The coupling parameter is roughly the ratio of the Coulomb interaction energy, evaluated at the mean interparticle separation, to the kinetic energy. The screening parameter is the ratio of the interparticle separation to the screening length.

Another inverse screening length-scale is associated with the particle component, . The quantity characterizes linear screening when both the particle and surrounding medium are responsible for it. Note that and, therefore, the relation between and takes the form .

The main quantities we will be dealing with in the following are the internal energy , Helmholtz free energy , and pressure , associated with the particle component. In reduced units these are

 u=U/NT,   f=F/NT,   p=PV/NT, (4)

where is the number of particles in the volume (so that ).

## Iii Debye-Hückel approximation

The Debye-Hückel (DH) approximation corresponds to the limit of extremely weak coupling, . The electric field around a (test) particle is screened due to rearrangement of neutralizing medium and the particles themselves. Linearizing Boltzmann distributions for the both components and substituting this into the Poisson equation yields

 ϕ(r)=(Q/r)exp(−kΣr) (5)

for the electrical potential distribution around the test particle. The reduced excess energy is

 uex=12QT[ϕ(r)−Qr]r→0=−12Q2kΣT. (6)

This can be easily rewritten in terms of and as

 uex(κ,Γ)=−12Γκ√1+3Γ/κ2. (7)

In the limit (where DH approximation is reliable) we have . In the one-component-plasma (OCP) limit, screening comes only from the particle component () and we recover the familiar result  Hansen ().

Other thermodynamic functions are easily obtained from the internal energy. For instance, the excess free energy in the fluid phase can be obtained via the integration

 fex=∫Γ0dΓ′uex(κ,Γ′)/Γ′. (8)

In the DH approximation this integration is straightforward and yields

 fex(κ,Γ)=−κ39[(1+3Γκ2)3/2−1]. (9)

The reduced pressure is

 p=n(∂f∂n)T. (10)

It is more convenient to rewrite this derivative in terms of and . In doing so we fix the density of the surrounding medium and observe that

 ∂Γ∂n=13Γn   and   ∂κ∂n=−13κn, (11)

since and . The pressure then becomes

 p=Γ3∂f∂Γ−κ3∂f∂κ. (12)

Note that this definition is different from that used in Ref. Farouki (), where the densities of neutralizing species (electrons and ions) were taken in constant proportion to (which resulted in the scaling ). In the next section we will explain why the present choice seems more appropriate.

Applied to the excess free energy of the DH model this relation yields the excess pressure

 pex=−12Γκ(1+3Γκ2)1/2+κ39[(1+3Γκ2)3/2−1]. (13)

In the OCP limit this reduces to the conventional expression .

The applicability of the DH approximation requires coupling to be small. To get an idea about its accuracy, let us compare the values of (i.e. at ) calculated with the help of equation (9) with the “exact” numbers obtained from molecular dynamics (MD) simulations in Refs. Hamaguchi (); Farouki (). This comparison is shown in Table  1. It is evident that even in the regime , the DH approximation is not characterized by high accuracy. Other approaches are needed and in the next Section we discuss one of the possible improvements.

## Iv Debye-Hückel plus hole approximation

The Debye-Hückel plus hole (DHH) approximation allows to reduce inaccuracy of the Debye-Hückel theory with respect to evaluating thermodynamics properties of moderately and strongly coupled OCP. The term “Debye-Hückel plus hole” is conventionally associated with the work by Nordholm Nordholm (), although similar arguments were used earlier Iosilevskiy (). The main idea behind the DHH approximation is that the exponential particle density must be truncated close to a test particle so as not to become negative upon linearization. Below we apply it to the model Yukawa system described in Section II.

The main equations of this approximation are as follows. The electrical potential around a test particle is given by the Poisson equation

 Δϕ=−4π(Qn−enm). (14)

The neutralizing medium (each species when multi-component) follows the Boltzmann distribution, which can be linearized. Other particles are absent in the sphere (hole) of radius around a test one. Outside the sphere, their density also follows the Boltzmann distribution which can be linearized. This can be written as

 n={0,r≤hn0(1−Qϕ/T),r>h. (15)

The quasineutrality condition (2) also holds.

For the potential inside the hole we then have yielding a general solution of the form

 ϕin(r)=(A1/r)exp(−kmr)+(A2/r)exp(kmr)+A3, (16)

where . Outside the sphere the potential satisfies , which gives the following solution vanishing at :

 ϕout(r)=(B/r)exp(−kΣr). (17)

The two solutions, Eqs. (16) and (17), should be matched at the hole boundary, which yields and . The two additional conditions are (implying that the particle density vanishes at the hole boundary) and (implying that tends to as ). This constitutes the full set of equations necessary to determine the hole radius as a function of and .

Introducing the reduced hole radius we get after some algebra the following transcendent equation for :

 xκ2[(1+√1+3Γ/κ2)ex+(1−√1+3Γ/κ2)e−x]+3Γ[(x−1)ex+(x+1)e−x]−2Γκ3=0. (18)

The reduced excess energy is which yields

 (19)

The first two terms on the right-hand side of Eq. (19) correspond to the particle-particle correlations in the DHH approximation, the third term represents the excess free energy of the surrounding medium, and the last term is the free energy of the sheath around each particle (see also Eq. (10) of Ref. Hamaguchi ()). Equation (19) can be therefore rewritten as

 uex=upp+um+ush.

Figure 1 shows the results of calculating the excess energy from Eqs. (18), (19) and comparison with the numerical results obtained in Ref. Hamaguchi () in the regime . Here the reduced excess energy is plotted versus the reduced coupling parameter , where is the coupling parameter at which fluid-solid phase transition occurs. The values of for a number of are tabulated in Table X of Ref. Hamaguchi (); various analytical fits for the dependence are also available Hamaguchi (); Vaulina (); KM (); JCP (). Figure 1 demonstrates that the DHH approximation is rather accurate up to . In this regime typical deviations of DHH from MD simulations do not exceed few percent. For stronger coupling, DHH systematically overestimates the (negative) excess energy. As fluid-solid phase transition is approached, the difference between DHH and MD simulations amounts to at , and reduces to at .

Figure 2 shows the comparison between excess energies calculated using DHH approximation and obtained in MD simulations in the regime . The qualitative picture remains the same as in the previously considered case. The DHH approximation provides good accuracy up to . Here the difference between DHH and MD results does not normally exceed . For stronger coupling DHH is again systematically overestimating . On approaching the fluid-solid transition the inaccuracy of DHH is for . As increases, the contribution from the sheath becomes dominant. In this regime and DHH becomes virtually more and more accurate. Already for one can hardly observe any difference between DHH and MD simulations in Fig. 2.

The excess free energy can be calculated from Eqs. (19) and (8). To get an idea how DHH can improve the conventional DH approximation in the weakly coupled regime, the values of have been calculated and listed in Table 1. Considerable improvement is evident. DHH approach underestimates the free energies from numerical simulations by approximately at and at .

Overall, we observe that DHH approximation reduces to the DH theory in the limit of weak coupling (), but remains relatively accurate up to , where DH theory is already grossly wrong. For even stronger coupling the accuracy is merely qualitative. This is not surprising, since the approach cannot catch the essential structural properties of the fluid state and thus cannot properly describe particle-particle correlations. Consequently, the DHH approximation is completely useless when thermodynamic quantities have to be known with sufficient accuracy (e.g. in the context of fluid-solid phase transition). In such cases direct numerical simulations are required Farouki (); Hamaguchi (). On the other hand, the appealing simplicity of the DHH approximation suggests to apply it when no such accuracy is necessary. One particular example will be given in Section VII.

The energy associated with particle-particle correlations can be also calculated via the energy equation HansenBook ()

 upp=2πn0T∫∞0r2V(r)g(r)dr, (20)

where is the radial distribution function and is the pair interaction energy. In the limit of weak correlations the pair distribution function tends to unity everywhere, . In this case we easily get so that and cancel each other exactly. The only remaining contribution to the excess free energy is due to sheaths . This will not give any contribution to the pressure, as evident from Eq. (12) and thus if particle-particle correlations are absent. This result, which has to be expected, critically depends on the model relation between and . In particular, non-zero excess pressure would be obtained if the model of Ref. Farouki () was used.

The DHH approximation is equivalent to the following form of [cf. Eq. (15)]

 g(r)={0,r≤h1−Qϕout(r)/T,r>h. (21)

It is straightforward to verify that integration in Eq. (20) then yields the first two terms () in the right-hand side of Eq. (19).

## V The OCP limit

The one-component-plasma is an idealized system of point charges immersed in a neutralizing uniform background of opposite charges. It corresponds to the limit of the model under consideration. Various aspects of the OCP systems have been extensively studied in the literature. Among them are thermodynamic properties and, especially, the equation of state. For the dependence of the excess energy on in the fluid phase it is conventional to use an expression of the form , which was obtained using the variational hard sphere approach (yielding the exponent DeWitt (). Later, it has been observed that the exponent yields somewhat better agreement with simulations Stringfellow (). The resulting fit for the excess energy in the fluid phase as proposed in Ref. HamaguchiJCP_1996 () is

 uex(Γ)=−0.899Γ+0.565Γ1/3−0.207−0.031Γ−1/3. (22)

Figure 1 demonstrates that it can be safely used in the regime , i.e for (we remind that in the OCP limit). On the other hand the linear DH approaxch is accurate only up to Hansen () and starts to overestimate significantly the actual energy at (see Fig. 1).

Concerning the DHH approximation, the hole radius is directly obtained from Eq. (18) expanding terms in series around . This yields

 kph=[1+(3Γ)3/2]1/3−1. (23)

The energy is then obtained from Eq. (19), where two terms survive in the considered limit, . Using Eq. (23) this becomes

 uex=−14{[1+(3Γ)3/2]2/3−1}. (24)

Equations (23) and (24) coincide with those in Refs. Nordholm (); Iosilevskiy (). In the limit of very small , Eq. (24) reduces to the DH result, but it remains adequate at much higher than the DH approach do. Figure 1 demonstrates that the DHH approximation provides reasonable agreement up to , where the fit (22) starts to work. In the strongly coupled regime , the DHH approximation yields the correct scaling , but the coefficient of proportionality is too low ( instead of ). Note that in this strongly coupled regime the simplest ion-sphere model, yielding  Baus (), reproduces the leading term of the fit (22) with impressive accuracy.

To conclude this section we point out that the hole radius defined by Eq. (23) represents the distance of the minimum separation between the OCP particles in the DHH approximation. The particle-particle interaction can thus be viewed as the strong short-range hard-sphere repulsion at plus weak long-range Debye-Hückel repulsion at . It is reasonable to apply linear plasma response formalism to describe momentum transfer in distant collisions between OCP particles. In doing so the inverse hole radius should be used as an estimate of the maximum wave vector entering into the kinetic definition of the Coulomb logarithm. The resulting effective Coulomb logarithm reduces to the classical expression in the regime of weak coupling, but remains meaningful for strong coupling, too. It has been recently shown that such an approach would describe reasonably the relaxation rate and the self-diffusion coefficient of OCP over the entire region of coupling, up to the fluid-solid transition CoulLog ().

## Vi Towards an equation of state

In this section we derive an equation for the excess pressure of the particle component in the DHH approximation. As has been discussed earlier, the excess free energy of the sheath does not contribute to the excess pressure. The excess pressure arising from the particle-particle correlations can be conveniently evaluated from the virial pressure equation HansenBook (). Since the terms corresponding to particle-particle correlations and neutralizing medium cancel each other exactly in the limit of no correlations [], the resulting expression for the excess pressure becomes

 pex=−2πn03T∫∞0r3V′(r)[g(r)−1]dr. (25)

Combining with the expression (21) for in the DHH approximation we get after some algebra

 pex=12Γκ2[−3+e−x(3+3x+x2)−−e−x1+√1+3Γ/κ2(x+x2+x1+√1+3Γ/κ2)]. (26)

Figure 3 shows the calculated as a function of for several values of (same values as used in Figs. 1 and 2). Excess pressure is always negative and decreases as the fluid-solid transition is approached. For the pressure curves are relatively close to each other, as has been already pointed out in Ref. Farouki (). This suggests to use a more accurate fit based on Eq. (22) in this regime (not shown in the Figure). With further increase in , the excess pressure tends to be less and less negative at a given value of . (Note however that comparison in terms of would apparently be more appropriate here). Thus, the OCP fit (22) becomes completely irrelevant for . We will discuss this issue further in the context of the particle density waves in Yukawa systems (dusty plasmas).

## Vii Dust Acoustic Waves

The minimalistic model for the dust acoustic waves (DAW) DAW () in dusty plasmas yields the following dispersion relation

 ω2ω2p=q2q2+κ2+q23Γγμp, (27)

where is the wave frequency, is the frequency-scale associated with the particle component (dust plasma frequency), is the reduced wavenumber, is the adiabatic index, and is the inverse reduced isothermal compressibility Kaw1998 (). The inverse compressibility is related to the excess pressure via

 μp=1+pex+Γ3∂pex∂Γ−κ3∂pex∂κ, (28)

The dispersion relation (27) can be easily derived using the Boltzmann response of the neutralizing medium along with the simplest hydrodynamic description of the particle component. Substituting these into the Poisson equation and linearizing will immediately yield Eq. (27). In complex (dusty) plasmas the neutralizing medium normally consists of positively charged ions and negatively charged electrons. Each component can be characterized by its own temperature, but this is not essential for the present consideration, since only the actual value of is affected. Note that in the limit Eq. (27) coincides with the phenomenological hydrodynamic dispersion relation of the OCP model (see e.g. Eq. (4.51) from Ref. Baus (); in this case, the leading terms yield ). Note also that in the original derivation of Ref. DAW () the particle component pressure was neglected at all, corresponding to (i. e. the assumption of cold particles was used).

Another dispersion relation suggested in the literature for DAWs in the strongly coupled regime reads Kaw1998 (); Kaw2001 ()

 ω2ω2p=q2q2+κ2+q23Γ(3+415uex), (29)

This type of dispersion relation originates from the sum-rule analysis of the OCP in the long-wavelength limit (e.g. Eq. (4.52) from Ref. Baus (); see also Eq. (21) from Ref. Abramo ()).

In the context of complex plasmas, both dispersion relations (27) and (29) neglect a number of properties specific to these systems, including e.g. collisions between different components (of particular importance are particle-neutral collisions), particle charge variations, external and internal forces acting on the particles (except the electrical one), etc. They are however appropriate for comparison with idealized computer experiments designed to study the effect of strong coupling on wave dispersion in Yukawa systems. We take the results of Ref. MD_DAW (), where wave dispersion relations in the fluid phase of Yukawa systems were obtained using molecular dynamics simulations. Simulations were performed for several state points characterized by certain values of and parameters. These state points are shown in Fig. 4, representing the sketch of the phase diagram of Yukawa systems. All the investigated state points correspond to rather strong coupling – they are located just below the melting curve.

Comparison between the numerical results and theory is shown in Fig. 5. Symbols and vertical bars represent simulation results for the longitudinal waves and their uncertainties. The solid curves correspond to the dispersion relation (27) with the compressibility evaluated using the DHH approximation (we also assume for such strong coupling). The dashed curves correspond to the dispersion relation (29) with the excess energy evaluated from the OCP fit of Eq. (22).

Figure 5 demonstrates that the theoretical OCP dispersion relation (29) fits nicely the numerical data in the regime . This is clearly the consequence of the fact that the presence of weak screening does not lead to significant deviations from the equation of state of the OCP, as has been pointed out in Ref. Farouki (). However, as screening becomes more pronounced (), the dispersion relation (29) becomes progressively less accurate. At it is completely off the simulation data. Note that this defect cannot be cured by taking into account screening and using more accurate results (either DHH or exact numerical) for the excess energy . The actual excess energy decreases (become more negative) with (see e.g. Table III from Ref. Hamaguchi ()), which implies that the situation would become even worse if exact values of are used instead of the OCP results. On the other hand, Eq. (27) is less accurate in the weakly screened regime (), in particular in the short-wavelength domain, but demonstrates reasonable agreement with the numerical results at stronger screening. Taking into account that the accuracy of the DHH approximation also increases in the strongly screened regime we can conclude that the simplest hydrodynamic description (27) combined with the DHH approximation for the inverse compressibility provides a reasonable compromise between the accuracy and simplicity/speed of calculations in the regime .

Another theoretical approach for the waves in strongly coupled Yukawa fluids (dusty plasmas) is based on the quasilocalized charge approximation Rosenberg (); Kalman (). The theory has been shown to agree very well with the numerical simulation data Kalman (). Since this approach also requires evaluation of the system internal energy, simple approximations similar to that considered in the present paper can again be of certain value.

We conclude this section with the following general observation. The excess pressure is a negative decreasing function of . Thus, the pressure term results in negative contribution to the dispersion relation, provided exceeds some critical value. In the OCP limit this transition is known as the onset of negative dispersion and recent numerical simulations locate it at  Mithen (). [Equations (27) and (29) yield and , respectively]. This implies that in the strongly coupled regime the group velocity becomes negative () at large as actually seen in all cases shown in Fig. 5. This feature is peculiar to the longitudinal modes in solids, which indicates that there is no qualitative difference between the dispersion properties of (strongly coupled) liquid and crystalline Yukawa systems (dusty plasmas). In the long-wavelength limit the longitudinal waves exhibit acoustic behavior () and their phase velocity is somewhat decreased due to the effect of strong coupling.

## Viii Discussion and Conclusion

We have discussed a simple analytical approach to estimate the thermodynamic properties of idealized Yukawa systems. The model considered consists of point-like charges embedded in a neutralizing medium, which is responsible for the exponential screening and the Yukawa pair interaction potential between the particles. Accurate numerical results exist for this model and these were used for comparison with our approximation. Although the obtained analytical results do not yield very high accuracy (in particular, in the regime of strong coupling and weak screening), they provide convenient formulas to describe the essential qualitative properties of Yukawa systems.

We note, however, that the idealized model does not account for some important properties of real systems. Some of these properties, which are relevant to complex plasmas are as follows: (i) Particles are not point-like, the typical ratio of the particle size to the plasma screening length can vary in a relatively wide range; (ii) There is a wide region around the particle where the ion-particle interaction is very strong, which results in non-linear screening; (iii) Particle charge is not fixed, but depends on complex plasmas parameters (e.g. on the particle density); (iv) The average density of ions and electrons is not fixed, but is related to the particle density and charge via the quasineutrality condition. Most of these properties are also to a large extent relevant to colloidal dispersions.

Clearly these properties can considerably affect the thermodynamics. From this perspective, accurate results for an idealized model can be considered as reference data for more advanced models. Extension of simple analytical approximations, similar to that discussed in this paper, would be a reasonable strategy to study the relative importance of the above mentioned properties. We leave this for future work.

###### Acknowledgements.
The authors would like to thank Satoshi Hamaguchi for providing numerical data on wave dispersion relations in Yukawa fluids. We appreciate funding from the Russian Foundation for Basic Research, Project No. 13-02-01099, and from the European Research Council under the European Union’s Seventh Framework Programme, Grant Agreement No. 267499.

## References

• (1) 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).
• (2) Complex and dusty plasmas: From Laboratory to Space, edited by V. E. Fortov and G. E. Morfill (CRC Press, Boca Raton, 2010).
• (3) H. Löwen, Phys. Rep. 237, 249 (1994).
• (4) V. E. Fortov, A. V. Ivlev, S. A. Khrapak, A. G. Khrapak, and G. E. Morfill, Phys. Rep. 421, 1 (2005).
• (5) G. E. Morfill and A. V. Ivlev, Rev. Mod. Phys. 81, 1353 (2009).
• (6) M. Bonitz, C. Henning, and D. Block, Rep. Progr. Phys. 73, 066501 (2010).
• (7) M. Chaudhuri, A. V. Ivlev, S. A. Khrapak, H. M. Thomas, G. E. Morfill, Soft Matter 7, 1287 (2011).
• (8) S. Hamaguchi and R. T. Farouki, J. Chem. Phys. 101, 9876 (1994).
• (9) R. T. Farouki and S. Hamaguchi, J. Chem. Phys. 101, 9885 (1994).
• (10) S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, J. Chem. Phys. 105, 7641 (1996).
• (11) S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, Phys. Rev. E 56, 4671 (1997).
• (12) G. J. Kalman, M. Rosenberg, and H. DeWitt, J. Phys. IV France 10, 403 (2000).
• (13) H. DeWitt, W. Slattery, D. Baiko, and D. Yakovlev, Contrib. Plasma Phys. 41, 251 (2001).
• (14) J. P. Hansen, Phys. Rev. A 8, 3096 (1973).
• (15) S. Nordholm, Chem. Phys. Lett. 105, 302 (1984).
• (16) V. K. Gryaznov and I. L. Iosilevskiy, Numerical methods in fluid mechanics 4, 166 (1973); For english translation see e-print arXiv:0903.4913 (2009).
• (17) O. S. Vaulina and S. A. Khrapak, JETP 92, 228 (2001); O. Vaulina, S. Khrapak, and G. Morfill, Phys. Rev. E 66, 016404 (2002).
• (18) S. A. Khrapak and G. E. Morfill, Phys. Rev. Lett. 103, 255003 (2009).
• (19) S. A. Khrapak, M. Chaudhuri, and G. E. Morfill, J. Chem. Phys. 134, 241101 (2011).
• (20) J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic Press, 2006).
• (21) H. E. DeWitt and Y. Rosenfeld, Phys. Lett. 75A, 79 (1979).
• (22) G. S. Stringfellow, H. E. DeWitt, and W. L. Slattery, Phys. Rev. A 41, 1105 (1990).
• (23) M. Baus and J. P. Hansen, Phys. Rep. 59, 1 (1980).
• (24) S. A. Khrapak, Phys. Plasmas 20, 054501 (2013).
• (25) N. N. Rao, P. K. Shukla, and M. Y. Yu, Planet Space Sci. 38, 543 (1990).
• (26) P. K. Kaw and A. Sen, Phys. Plasmas 5, 3552 (1998).
• (27) P. K. Kaw, Phys. Plasmas 8, 1870 (2001).
• (28) M. C. Abramo and M. P. Tosi, Il Nuovo Cimento 21, 363 (1974).
• (29) H. Ohta and S. Hamaguchi, Phys. Rev. Lett. 84, 6026 (2000); S. Hamaguchi and H. Ohta, Phys. Scripta T89, 127 (2001).
• (30) M. Rosenberg and G. Kalman, Phys. Rev. E 56, 7166 (1997).
• (31) G. Kalman, M. Rosenberg and H. E. DeWitt, Phys. Rev. Lett. 84, 6030 (2000).
• (32) J. P. Mithen, J. Daligault and G. Gregory, AIP Conf. Proc. 1421, 68 (2012).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters