Thermodynamics and dynamics of two-dimensional systems with dipole-like repulsive interactions

Thermodynamics and dynamics of two-dimensional systems with dipole-like repulsive interactions

Sergey A. Khrapak Aix Marseille University, CNRS, PIIM, 13397 Marseille, France Institut für Materialphysik im Weltraum, Deutsches Zentrum für Luft- und Raumfahrt (DLR), 82234 Weßling, Germany Joint Institute for High Temperatures, Russian Academy of Sciences, 125412 Moscow, Russia    Nikita P. Kryuchkov Bauman Moscow State Technical University, 2nd Baumanskaya street 5, 105005 Moscow, Russia    Stanislav O. Yurchenko Bauman Moscow State Technical University, 2nd Baumanskaya street 5, 105005 Moscow, Russia
July 14, 2019

Thermodynamics and dynamics of a classical two-dimensional system with dipole-like isotropic repulsive interactions are studied systematically using extensive molecular dynamics (MD) simulations supplemented by appropriate theoretical approximations. This interaction potential, which decays as an inverse cube of the interparticle distance, belongs to the class of very soft long-ranged interactions. As a result, the investigated system exhibits certain universal properties that are also shared by other related soft-interacting particle systems (like, for instance, the one-component plasma and weakly screened Coulomb systems). These universalities are explored in this article to construct a simple and reliable description of the system thermodynamics. In particular, Helmholtz free energies of the fluid and solid phases are derived, from which the location of the fluid-solid coexistence is determined. The quasi-crystalline approximation is applied to the description of collective modes in dipole fluids. Its simplification, previously validated on strongly coupled plasma fluids, is used to derive explicit analytic dispersion relations for the longitudinal and transverse wave modes, which compare satisfactory with those obtained from direct MD simulations in the long-wavelength regime. Sound velocities of the dipole fluids and solids are derived and analyzed.

I Introduction

Two- and quasi-two-dimensional (2D) interacting particle systems attract great scientific interest, since they play an important role in a broad range of phenomena operating at fluid and solid surfaces and various interfaces Kosterlitz and Thouless (1978); Kosterlitz (2017). Several relevant examples include atomic monolayers and thin films on a substrate, two-dimensional electron gas on the surface of liquid helium, vortices in thin-film semiconductors, metallic and magnetic layer compounds, smectic liquid crystals, colloidal particles at flat interfaces, and complex (dusty) plasma systems in ground-based conditions.

In the context of colloidal systems, apart from technological applications of colloidally-stabilized emulsions Clegg (2008); Stocco et al. (2009); Garbin et al. (2015) and bubbles Poulichet and Garbin (2015); Yurchenko et al. (2016a) for synthesis of novel optical materials van der Meer et al. (2016); Elsner et al. (2009); Gray et al. (2015); Garbin et al. (2015), chemical sensors, and catalysis Lee et al. (2004), many biologically important processes occur at interfaces. Importantly, there is a way to control these processes by attaching colloidal particles to soft matter interfaces.

When colloidal particles are trapped in oil/water or gas/water interfaces, electrical dipoles are usually associated with each interfacial particle Pieranski (1980). As a result, the interaction between colloidal particles is similar to that between vertically oriented dipoles Goulding and Hansen (1998); Oettel (2007); Park and Furst (2011); Kelleher et al. (2015, 2017) and can be in the first approximation described by a pairwise repulsive inverse-power law (IPL) potential decaying as with the interparticle separation . Direct experimental measurements of colloidal interaction potential in such systems by the laser tweezers method Park and Furst (2011); Aveyard et al. (2002); Kelleher et al. (2015) or using other approaches Bonales et al. (2011) generally confirm this assumption, although it is also clear that actual interactions can be very complicated, particularly in the regime where the separation is comparable to the particle size Girotto et al. (2016). A similar shape of the interaction potential is observed in two-dimensional colloidal systems of paramagnetic particles exposed to an external magnetic field Zahn et al. (1999); Zahn and Maret (2000); Gasser et al. (2010). 2D colloidal suspensions in external electric fields represent another important class of dipole-like interacting systems. An external electric field polarizes the particles and ion clouds in the solvent around them, inducing a (tunable) dipole-dipole interaction between the particles. Depending on the orientation of the external electric field with respect to the plane of particle confinement, the dipolar interaction potential can be either attractiveJuárez and Bevan (2009); Carstensen et al. (2015); Yakovlev et al. (2017); Ovcharov et al. (2017) or repulsive Dutcher et al. (2013); Tanaka et al. (2014); Ma et al. (2014); Woehl et al. (2015); Saini et al. (2016); Lotito and Zambelli (2017); Gong and Wu (2017).

In the context of plasma physics, it has been long known that the effective potential of a point test charge immersed in a flowing collisionless plasma is not screened exponentially, but falls off as , at large distances from the test charge Montgomery et al. (1968); Cooper (1969). This can be relevant to complex (dusty) plasmas, a collection of small solid particles in the neutralizing plasma medium. In a typical laboratory dusty plasma experiment the highly negatively charged identical micron-size particles form a horizontal (quasi 2D) layer above the bottom negatively biased electrode of a radio-frequency gas discharge, where the electric force directed upwards is able to balance the gravity force acting on the particles. A strong electric field required to balance the gravity produces significant ion flow, which makes electric potential distribution around the particles highly anisotropic. Although the actual interactions between the particles in these conditions are quite sophisticated and are governed by a competition between screening and plasma-wake mediated effects Vladimirov and Nambu (1995); Ivlev et al. (2012); Fortov et al. (2004, 2005); Chaudhuri et al. (2010); Hutchinson (2011); Ludwig et al. (2012), there is a certain parameter regime, where the IPL scaling is relevant Kompaneets et al. (2009, 2016a, 2016b) (see, in particular, Fig. 3 in Ref. Kompaneets et al. (2016b)).

Thus, dipole-like interactions occur in various two-dimensional physical systems such as ions and colloidal particles trapped at various interfaces, colloidal particles in external electric fields, paramagnetic particles exposed to external magnetic fields, electrical charges placed in a flowing collisionless plasma, etc. Not surprisingly, structural and dynamical properties, thermodynamics, phase transitions, collective motion and related phenomena in classical systems with repulsion have been extensively studied (see, for instance, Refs. von Grünberg et al. (2004); Keim et al. (2004); Lin et al. (2006); Löwen (1996); van Teeffelen et al. (2006, 2008); Golden et al. (2010, 2008a, 2008b); Dillmann et al. (2012); Haghgooie and Doyle (2005); Jaiswal et al. (2013) and references therein). The main purpose of this work is to put strong emphasis on the fact that the considered dipolar interaction belongs to the class of very soft long-ranged interactions, the limit opposite to the celebrated hard sphere interaction in three-dimensions (3D) and hard disc interaction in 2D. Based on this, a simple description of thermodynamic and dynamic properties is possible, using methods validated recently on other classical soft interacting particle systems, mainly in the plasma-related context.

Systems of soft interacting particles exhibit certain universal properties and there exist useful approximations, that are particularly suitable for this regime. In particular, the Rosenfeld-Tarrazona (RT) scaling Rosenfeld and Tarazona (1998); Rosenfeld (2000) of the thermal component of the internal excess energy on approaching the freezing transition allowed previously to construct a very simple practical approach to the thermodynamics of weakly screened Yukawa systems in 3D Khrapak and Thomas (2015a); Khrapak et al. (2015a); Khrapak (2016a). An analog of the RT scaling also exists in the 2D case (although of a quite different functional form) and this has been recently used to construct a simple thermodynamic description of one-component plasmas and weakly screened Yukawa systems in 2D Khrapak et al. (2015b); Semenov et al. (2015); Khrapak and Khrapak (2016); Kryuchkov et al. (2017) with main applications to complex (dusty) plasmas. Here we apply the same arguments to the 2D system with dipolar interactions to put forward simple and accurate expressions for the thermodynamic properties of the liquid state, which are (by construction) in excellent agreement with the MD simulation results. Combined with the accurate calculation of the thermodynamic functions of the crystalline solid (using MD simulations and the shortest graphs method, proposed recently by some of the present authors) we are also able to approximately locate the fluid-solid phase transition, as well as the narrow coexistence region.

We also elaborate on the properties of collective modes in 2D dipolar systems. Recent investigations demonstrated that the quasi-crystalline approximation (QCA) Hubbard and Beeby (1969); Takeno and Gôda (1971), also referred to as the quasi-localized charge approximation (QLCA) in the plasma-related context Golden and Kalman (2000), is a good approximation to describe elastic collective modes in dense fluids for the regime of soft interactions (though it fails in the limit of very steep hard-sphere/hard-disc interactions Khrapak et al. (2017)). Previously, QLCA has been applied with certain success to dipole-like systems in 2D Golden et al. (2009, 2010). Here we go somewhat further and derive simple analytic expressions, describing well the dispersion relations of the longitudinal and transverse elastic modes at sufficiently long wavelengths. The accuracy of these dispersion relations is demonstrated by comparing with the dispersions obtained from MD simulations. We demonstrate how these results can be useful in estimating the free energy of the crystalline solid. We also evaluate the high-frequency elastic moduli of the considered system and discuss relations to sound velocities operating in the strongly coupled fluid regime.

The rest of the article is organized as follows. In Section  II we describe in detail the system under investigation, provide necessary details about the performed MD simulations, and summarize main thermodynamic relations used in this work. In section III main results obtained for the fluid phase are summarized, including accurate expressions for thermodynamic quantities and detailed analysis of collective modes. In Section IV topics related to the thermodynamics of the crystalline phase are addressed. This includes thermodynamic functions, location of the fluid-solid phase transition, and sound velocities of an idealized lattice. Section V is focused on elastic moduli and their relations to the sound velocities in a dense fluid phase. This is followed by our conclusion in Section VI.

Ii Methods

ii.1 System description

We investigate a classical system of point-like particles in the 2D geometry, which are interacting via the pairwise repulsive inverse-third-power (IPL3) potential of the form


where and are the energy and length scales of the interaction. Phase behavior is conveniently described by the dimensionless interaction (coupling) parameter ,


where is the temperature (in energy units), is the 2D Wigner-Seitz radius, and is the areal density of particles occupying the 2D volume (i.e. surface) . The coupling parameter is roughly the ratio of the potential energy of interaction between two neighboring particles to their kinetic energy. The system is conventionally referred to as strongly coupled when the potential energy dominates, that is when .

At very low the system properties are similar to those of an ideal gas in 2D. When coupling increases the system forms a strongly coupled fluid phase, which can crystallize upon further increase in . The nature of the fluid-solid phase transition in 2D systems depends considerably on the potential softness Kapfer and Krauth (2015). For sufficiently steep repulsive interactions the hard-disk melting scenario holds: a first-order liquid-hexatic and a continuous hexatic-solid transition can be identified Bernard and Krauth (2011); Engel et al. (2013); Thorneywork et al. (2017). For softer interactions the liquid-hexatic transition is continuous, with correlations consistent with the Berezinsky-Kosterlitz-Thouless-Halperin-Nelson-Young (BKTHNY) scenario Kapfer and Krauth (2015); Kosterlitz (2017). For the IPL family of potentials () the transition between the two regimes occurs at about  Kapfer and Krauth (2015). The IPL3 system studied here thus belongs to the soft interaction class and the BKTHNY melting scenario should apply. This indeed has been observed both in numerical simulations Lin et al. (2006) and colloidal experiments Zahn et al. (1999); Zahn and Maret (2000); Deutschländer et al. (2014); Kelleher et al. (2015). However, the hexatic phase occupies a relatively narrow region on the phase diagram and its properties will not be addressed in this work.

ii.2 Computational details

To obtain the accurate thermodynamic properties of IPL3 fluids and crystals in 2D, as well as necessary information about the properties of collective modes, extensive MD simulations have been performed using LAMMPS package LAM (). These MD simulations have been done in the ensemble at different temperatures using particles in a simulation box with periodic boundary conditions and the Nose-Hoover thermostat. The initial systems configuration was chosen as an ideal hexagonal lattice and velocities were set according to the Maxwell distribution with the temperature equal to and for the fluid and crystal phases, respectively. The numerical time step of was chosen. All simulation runs were performed for time steps, where the last steps were used for energy and pressure calculation based on standard functions implemented in the LAMMPS package (in the case of fluids, to guarantee that equilibrium was reached, we performed three simulations with different initial conditions for each examined state point). The cutoff radius of the potential was set equal to . The internal energy and pressure were corrected accordingly, which resulted in a relative error of about , sufficient for the range of problems considered here.

ii.3 Thermodynamic relations

The main thermodynamic quantities of interest in this work are the internal energy , Helmholtz free energy , and pressure of the system. The following thermodynamic definitions are useful Landau and Lifshitz (2013a)


In addition, and can be calculated using the integral equations of state Hansen and MacDonald (2006); Frenkel and Smit (2001)


where denotes the radial distribution function, which is isotropic in gas and fluid phases and anisotropic in the crystalline phase.

We will use conventional reduced units: , , and and divide the thermodynamic quantities into the kinetic (ideal gas) and potential (excess) contributions, so that (in 2D), , and . Finally, it is useful to operate with the single coupling parameter , instead of temperature and density. Since , the transformation of standard thermodynamic relations to their dimensionless form is governed by


where is a thermodynamic function of interest. In addition, a simple relation between the reduced excess pressure and energy for the IPL3 interaction in 2D holds:


Other thermodynamic quantities can be readily evaluated when the excess internal energy is known. We summarize the main relations employed in this work in Appendix A.

Iii Fluids

iii.1 Thermodynamics of the fluid phase

The excess energy and pressure of the 2D IPL3 fluid have been determined using MD simulations. Based on these results, combined with our previous experience with thermodynamics of soft interacting particle systems in 2D geometry (mostly in the plasma-related context), simple and reliable analytical approximations are derived.

In the strongly coupled regime it is helpful to divide the thermodynamic properties, such as energy and pressure, into static and thermal contributions. The static contribution corresponds to the value of internal energy when the particles are frozen in a regular configuration and the thermal corrections arise due to deviations of the particles from these fixed positions (due to thermal motion). Here we relate the static energy to the lattice sum of the triangular lattice (Madelung energy) formed by particles interacting via the potential (this relation is meaningful for both crystalline and fluid phases). The corresponding lattice sum has been evaluated previously with a very high accuracy Topping (1927); van der Hoff and Benson (1953). The proportionality constant between the static energy and the coupling parameter (Madelung constant) is

Thus, the excess internal energy in the fluid phase can be expressed as


The usefulness of this approximation stems from the fact that the ratio of the thermal-to-static contribution is small for strongly coupled fluids with soft long-ranged interactions. In this case the static part is dominated by the cumulative contribution from large interparticle separations. This part is not very sensitive to the actual short-range order in the system since for large separations . It also does not change much across the fluid-solid phase transition, and thus the Madelung energy is an appropriate characteristic for both solid and fluid phases. Quantitatively, the static contribution is much larger than the kinetic energy (), by the definition of strong coupling. On the other hand, the thermal contribution comes from the particle thermal motion and its magnitude is of the order of the average particle kinetic energy (). This implies that even moderately accurate approximations for result in a very accurate estimation of the total excess energy of strongly coupled fluids. The remaining step is therefore to identify an appropriate approximation for .

Based on the previous results for other soft interacting particle systems in 2D (such as one-component-plasma Khrapak and Khrapak (2014, 2016) and weakly screened Yukawa systems Khrapak et al. (2015b); Kryuchkov et al. (2017)), the thermal component of the excess energy is expected to exhibit a certain scaling with on approaching the fluid-solid transition. This scaling is to some extent analogous to the RT scaling of the thermal component of excess energy in 3D Rosenfeld and Tarazona (1998); Rosenfeld (2000), but has a quite different functional form. The functional form suggested for 2D systems with soft pairwise repulsive interactions is


The validity of this functional form at sufficiently strong coupling is documented in Fig. 1, where numerical data from the present MD simulations along with those reported previously Golden et al. (2008a, b) are plotted. The best fit of the MD data obtained in this work yields and . The fit is valid in the range . Combining Eqs. (8) and (9) we write for the excess energy of the strongly coupled fluid phase


It is easy to ascertain that the first term is indeed dominant at strong coupling. For example, the ratio of the second to the first terms in Eq. (10) is at , it decreases to at , and further drops to at .

Figure 1: Thermal component of the reduced excess energy, , of a strongly coupled IPL3 system in 2D versus the coupling parameter . Circles correspond to the results of MD simulations performed in this work. Triangles are the results by Golden et al. Golden et al. (2008a, b). The curve is the analytical fit of Eq. (9).

Having a fit for the reduced excess energy, we can evaluate the reduced excess Helmholz free energy from the relation


However, one must pay some attention to the fact that expression (9) is not applicable all the way down to . Although the actual contribution to the free energy from the weak coupling regime is of minor importance at strong coupling, we have accounted it in the following manner. At very weak coupling, the virial expansion allows us to estimate the free energy with a reasonable accuracy. The first order correction to the ideal free energy is Landau and Lifshitz (2013a)


which can be evaluated analytically. The excess energy in this regime is . However, careful comparison with the results from numerical simulations shows that Eq. (12) is applicable only at very weak coupling, . On the other hand, the strong coupling scaling (10) is justified only for . An approximation for the intermediate regime, basically based on an appropriate combination of Eqs. (10) and (12) is described in Appendix B. Using this approximation [Eq. (39)] below and Eq. (10) at higher values of , and performing integration in Eq. (11) we finally obtained a simple and accurate analytical approximation for the fluid free energy in the strong coupling regime, :


where is a polylogarithm function. It is the last constant term, which is responsible for the contribution from the weak-coupling regime. Clearly, the first two terms are dominant at strong coupling.

Equations (10) and (13) represent our main results regarding thermodynamics of the IPL3 fluids in 2D. We will estimate below that the fluid-solid phase transition should occur at , with a very narrow coexistence gap (see Section IV). Some thermodynamic quantities of the IPL3 melt (fluid just at the boundary of the fluid-solid coexistence) obtained using the approach discussed here are summarized in Appendix C.

iii.2 Collective modes

It is well known that fluids can exhibit different collective dynamics depending on the regime of coupling and correlations Hansen and MacDonald (2006); Trachenko and Brazhkin (2016); Fomin et al. (2016); Yang et al. (2017). In the regime of weak correlations the dynamics is close to that in the ideal gas, and there exists only the longitudinal collective mode. On the other hand, in dense liquids not too far from the melting line, where interparticle correlations are strong, the transverse mode (one mode in the 2D case and two modes in the 3D case) can also be excited in addition to the longitudinal mode. It is this latter regime that will be mostly considered below.

A powerful theoretical approach to describe collective motion in classical systems of strongly interacting particles with soft pairwise interactions is the quasi-crystalline approximation Hubbard and Beeby (1969); Takeno and Gôda (1971). This approach can be considered as either a generalization of the random phase approximation or, alternatively, as a generalization of the phonon theory of solids (the latter explains why it is often referred to as QCA). In the context of plasma physics an analog of the QCA is known as the quasi-localized charge approximation, QLCA. It was initially proposed as a formalism to describe collective mode dispersion in strongly coupled charged Coulomb liquids Golden and Kalman (2000). In recent years it was successively applied to strongly coupled one-component plasma in both 2D and 3D Golden and Kalman (2000) as well as 2D and 3D Yukawa fluids, mostly in the context of complex (dusty) plasmas Rosenberg and Kalman (1997); Kalman et al. (2000, 2004); Hartmann et al. (2007); Donko et al. (2008); Hou et al. (2009); Upadhyaya et al. (2010, 2011). Application to the 2D IPL3 system was described in Refs. Golden et al. (2010, 2009). Here we discuss a procedure to derive simple explicit expressions for the longitudinal and transverse dispersion relations. We demonstrate that these dispersion relations are reasonably accurate at long wavelengths using the comparison with the results from MD simulations. Then, we also briefly discuss how the obtained results can be used to estimate the free energy of the IPL3 solid.

Within the QCA approach, the dispersion relations of elastic waves at strong coupling are directly expressed in terms of the radial distribution function (RDF), , and the pair interaction potential . The compact expressions are Hubbard and Beeby (1969); Takeno and Gôda (1971); Zwanzig (1967)


where is the frequency, is the wave number, and is the direction of the propagation of the longitudinal mode (the particles are confined to the plane). Here the subscripts “” and “” correspond to the longitudinal and transverse modes, respectively. The explicit expressions in the 2D case along with the expressions for the special case of the IPL3 system are provided in Appendix D for completeness.

In the long-wavelength () regime the dispersion relations of the IPL3 fluid exhibit acoustic dispersion and the longitudinal and transverse sound velocities can be introduced,


Similarly to the IPL system in the 3D case Khrapak et al. (2017), these sound velocities can be easily related to the reduced excess energy (or pressure). For the considered case of IPL3 in 2D we immediately get for the QCA elastic sound velocities:


where is the thermal velocity of the particles. Here we used the relation that follows directly from energy or pressure equations (5),

where is the conventional 2D frequency scale. Note an immediate consequence of Eq. (17), , for the IPL3 system in 2D.

A simplest model , which takes into account the existence of a correlational hole (which prevent strongly repulsive particles from closely approaching each other) and is unity at longer separation (where correlations are small), turns out to be quite useful for soft repulsive interactions. Mathematically, this simplest model RDF reads


where is the Heaviside step function and the radius of the correlational hole is of order unity in reduced units (the distances are expressed in units of here). A similar RDF was employed previously to analyze the main tendencies in the behavior of specific heat of liquids and dense gases at low temperatures Stishov (1980). It was also used to calculate the dispersion relation of Coulomb bilayers and superlattices at strong coupling Golden and Kalman (1993). Physically, the model form of Eq. (18) seems sensible in the present context, because the main contribution to the long-wavelength dispersion corresponds to long length-scales, where . For soft enough interactions, this regime provides dominant contribution to the integrals in Eqs. (43) and (44). The excluded volume effect for allows us to properly account for strong coupling. In addition, an appealing advantage of this simple RDF is that when substituted in the QCA (QLCA) expressions, it allows the analytical integration for certain interaction potentials. Particularly simple and elegant expressions have been recently derived for Yukawa systems and one-component plasma in 3D Khrapak et al. (2016a); Khrapak (2017); Khrapak and Khrapak (ress) and one-component plasma with logarithmic interactions in 2D Khrapak et al. (2016b).

For the considered IPL3 system in 2D the resulting expressions are not so elegant, although still tractable. We get for the longitudinal mode


where is the reduced wave number, and are Bessel functions of the first kind, and are the Struve functions of order and , respectively. The dispersion relation of the transverse mode is remarkably more simple,


The remaining step is the determination of the appropriate correlational hole radius . Previously, it was proposed to determine from the condition that the model form (18) delivers good accuracy when substituted into the energy and/or pressure equations Khrapak et al. (2016a). This has been demonstrated to work well for both Yukawa and OCP systems in 3D and 2D OCP with logarithmic interactions Khrapak et al. (2016a, b). Following the same procedure we obtain


It is straightforward to verify that with this definition of , the low- series expansion of Eqs. (19) and (20) will reproduce the acoustic velocities given by (17). Note also that in the strong coupling regime the excess energy is mainly associated with the static contribution, . In this regime the radius of the correlational hole is practically constant, .

Figure 2: Longitudinal and transverse wave dispersion relations in color-coded format as obtained from Eq. (22) for (a) - (c) and (d) - (f). Top panel [(a), (d)] corresponds to the longitudinal mode, bottom panel [(b), (e)] to the transverse mode. Circles, and symbols correspond to the frequencies and values, respectively, obtained by applying a fitting function (23) for every fixed value of dimensionless wave vector . White dashed lines correspond to the acoustic asymptotes and , where and are the longitudinal and transverse sound velocities. For detailed discussion about the sound velocities in the IPL3 fluid, see Sec. V. The dotted red lines correspond to the short-wavelength kinetic asymptote , with (see the text). The insets (c) and (f) show the long-wavelength portions of the dispersion relations. Here the red (blue) circles correspond to the longitudinal (transverse) dispersion relations as obtained from MD simulations. The red and blue solid curves display calculations using QCA approach of Eqs. (43) and (44) with the RDFs obtained from MD simulations. The dashed black curves correspond to simple analytical expressions of Eqs. (19) and (20).

In order to verify the quality of this simple analytical approximation, the dispersion relations of the IPL3 fluid have been obtained from MD simulations. We used the standard approach to compute phonon spectra in fluids, which is based on the longitudinal and transverse current correlation functions Golden et al. (2010); Ohta and Hamaguchi (2000); Fomin et al. (2016). Specifically, we calculated


where and are the projections of velocity current to longitudinal and transversal directions, respectively. Here the summation is performed over all particles in the system, is radius-vector of the -th particle and is its velocity. Since fluids are isotropic we can average over all directions of the wave vector to get the dependence on .

Figure 2 shows, in color-coded format, the dispersion relations of the longitudinal and transverse waves obtained in this way for the two fluid state points characterized by (a) - (c) and (d)-(f). In contrast to crystals, color coding of current fluctuation spectra for fluids can merely be used to illustrate qualitative properties. To get more quantitative information we fitted by the Cauchy distribution,


for each value of . Examples of the obtained dependencies and are plotted in Fig. 2.

The long-wavelength portions of the dispersion relations obtained from MD simulation are plotted in Figs. 2(c) and 2(f). Here they are compared with QCA dispersion relations. The solid curves correspond to the “full” QCA with the actual RDF obtained in MD simulations and substituted in Eqs. (14) and (15). The black dashed curves correspond to the “simplified” QCA given by analytical expressions (19) and (20). The agreement between the two versions of QCA and MD dispersions is satisfactory at sufficiently long wavelengths ( for the longitudinal and for the transverse mode). This (low-) regime corresponds to long length-scales, where both actual and model RDF are similar, . For shorter wavelengths the two versions of QCA behave differently. This regime corresponds to short distances and the correct account of short-range correlations existing in liquids is necessary. Not surprisingly, the full QCA with the actual RDF is more consistent with MD-generated dispersion relations.

Nevertheless, clear disagreement between the QCA and MD spectra is still observed at short wavelengths even with the use of accurate . The main reason for this is that QCA does not take into account effects of anharmonicity, which are causative, in particular, for damping of collective excitations. Indeed, the particles in liquid are considered within the framework of QCA as “frozen” near their equilibrium positions, whose statistics is determined by the actual RDF . Then, the excitation spectra are calculated in the harmonic approximation using perturbation theory for small displacements of particles around equilibrium positions. Account of particles’ jumps (important for the physics of liquid state) cannot be done within the framework of perturbation theory Trachenko and Brazhkin (2016). At the same time, anharmonicity is related to the short-range region of the interaction potential, which corresponds to large in the reciprocal space and results in the observed growing discrepancy between the QCA and MD spectra. It should be also pointed out that the disappearance of the transverse mode at long wavelengths and the existence of a -gap for the transverse waves propagation Yang et al. (2017) cannot be properly described within the conventional QCA approach, because damping effects (associated with anharmonic interactions) are not included. One can see from examples presented in Figs. 2(b) and (e) that the width of the -gap decreases with increasing correlations (i.e. increasing ), in agreement with previous studies on collective excitations in various kinds of liquids Trachenko and Brazhkin (2016); Yang et al. (2017); Schmidt et al. (1997); Bryk et al. (2017).

Regarding the longitudinal mode, the asymptotic behavior of the dispersion relation at is not described by QCA, because the kinetic effects are missing in this approximation. In this regime the characteristic scales of particle motion are much less than the average separation between the particles. The expected high- asymptote is consistent with MD simulation results. The proportionality coefficient has been previously suggested in Ref. Golden et al. (2010). Our preliminary analysis indicates that the coefficient can be more appropriate. The small () relative difference between these coefficients does not allow to discriminate between these two values using the obtained MD data. We, therefore, leave this issue for future work.

Iv Crystals

The reduced excess energy of a 2D crystalline lattice in the harmonic approximation is . We need to add a small anharmonic correction to get the total excess energy of a crystalline phase. This anharmonic correction has been evaluated using MD simulations, and the results are plotted in Figure 3. The anharmonic corrections are fitted using the standard functional form Dubin and O’Neil (1999)


The coefficients determined from the fitting are , , and (curve in Fig. 3). The resulting excess internal energy of the solid phase is


The reduced excess Helmholtz free energy can then be evaluated in the following manner Hamaguchi et al. (1997). First, we divide it into anharmonic and harmonic contributions


where is the reduced excess free energy in the harmonic approximation. It is calculated by adding the lattice and vibrational free energies and subtracting the free energy of the perfect gas Landau and Lifshitz (2013a); Alastuey and Jancovici (1981). In 2D geometry the resulting expression for the reduced harmonic free energy is Alastuey and Jancovici (1981)


where is the frequency of a phonon with wavenumber and polarization , and the sum on is over the first Brillouin zone in the reciprocal lattice. The sum of the last two terms is sometimes referred to as the harmonic entropy constant , which is determined by the phonon spectrum of the crystalline lattice. The latter has been evaluated for the IPL3 triangular lattice using the standard technique, the resulting phonon dispersion curves are shown in the inset of Fig. 4. The corresponding harmonic entropy constant has been evaluated as (a related approach to estimate using the QCA dispersion relations is briefly discussed in Appendix E). Thus, the reduced Helmholz free energy of the IPL3 crystalline solid is


This is our main result concerning the thermodynamics of the solid phase.

We can now estimate the location of the fluid-solid phase transition in the 2D IPL3 system by equating Helmholtz free energies of the corresponding phases. This yields (here the subscript “m” refers to melting). In a more detailed consideration we equate fluid and solid temperatures, pressures and chemical potentials to evaluate the location and width of the phase coexistence region. The standard procedure then yields and . This is comparable to the results previously reported in the literature Zahn et al. (1999); Zahn and Maret (2000); Löwen (1996); van Teeffelen et al. (2006); Jaiswal et al. (2013) and is particularly close to the results from the Brownian dynamics simulations Haghgooie and Doyle (2005). Note that a very narrow coexistence gap obtained, (here ), should be related to the very soft character of the interaction potential. The anharmonic terms are not very important for the location of the phase transition: With neglecting anharmonic corrections, the fluid and solid free energy intersection point moves to . As a final remark, we note that we have not considered the existence of the hexatic phase.

Figure 3: Anharmonic corrections to the reduced excess energy of the IPL3 crystal in 2D versus the inverse coupling parameter . Symbols represent the results from our MD simulation, the solid curve is a fit of Eq. (24).

Thermodynamic properties of the IPL3 crystals can also be evaluated based on purely theoretical approach using an interpolation method (IM) for pair correlations in classical crystals proposed recently by some of the present authors Yurchenko (2014); Yurchenko et al. (2015, 2016b). This approach allows us to compute RDF in the crystalline state based on the Born-von Karman (BvK) phonon spectrum and taking into account anharmonic corrections to the first correlation peak. The technical details of this approach are summarized in Appendix F.

In Fig. 4 an example of the crystalline RDF is presented. Visual similarity between the obtained theoretical RDF and MD data is the same as in previous applications of the IM approach Kryuchkov et al. (2017); Yurchenko et al. (2015, 2016b). Using these highly accurate RDFs, pressure and excess energy can be readily obtained using Eqs. (5). In its simplest harmonic form () the IM approach provides a relative error in the excess energy smaller than in the worst case near the melting point. Taking into account anharmonicity, with the anharmonic correction coefficient [see Eq. (51)], obtained from MD simulations of the IPL3 crystal, reduces the relative error to .

The longitudinal and transverse sound velocities of a perfect IPL3 2D lattice are given by Eqs. (17) combined with (we remind that QCA reduces to the standard phonon theory of solids in the limit ). The final result is

Figure 4: Example of the radial distribution function for the IPL3 crystalline lattice at . The symbols correspond to the MD simulation results, the solid curve is obtained using the IM approach. The inset demonstrates BvK phonon dispersion curves for the IPL3 triangular lattice.

V High frequency elastic moduli and sound velocities in the liquid state

The high frequency (instantaneous) elastic moduli for simple 3D fluids were derived by Zwanzig and Mountain Zwanzig and Mountain (1965). The 2D analogues of these elastic moduli are


the high frequency limit of the bulk modulus Khrapak (2016b), and


the high frequency limit of the shear modulus. The relations between the QCA elastic sound velocities derived in Section III.2 and the elastic moduli are:


Useful relations between and in both 3D and 2D have been recently discussed Khrapak (2016c). In particular, the high-frequency (instantaneous) sound velocity , directly related to the instantaneous bulk modulus, can be introduced


The main purpose of this Section is to demonstrate that this instantaneous sound velocity is extremely close to the conventional adiabatic sound velocity appearing in hydrodynamic description of fluids Landau and Lifshitz (2013b)

10 20 30 40 50 60 70
6.04 8.38 10.18 11.70 13.04 14.25 15.37
5.93 8.11 9.81 11.24 12.51 13.66 14.72
5.92 8.10 9.80 11.24 12.51 13.66 14.72
Table 1: Reduced sound velocities (in units of thermal velocity) of the IPL3 fluid in 2D evaluated for different values of the coupling parameter , corresponding to the strongly coupled fluid phase.

The sound velocities , , and for several values of the coupling parameter , corresponding to the strongly coupled fluid phase, are summarized in Table 1. It is observed that overestimates the adiabatic sound velocity . However, the difference is rather small, as should be expected for the soft interaction potential studied in this work Khrapak (2016a); Khrapak et al. (2017); Khrapak and Thomas (2015b). (For soft interactions at strong coupling one normally observe and , which implies  Khrapak et al. (2017)). The instantaneous sound velocity is just slightly above the adiabatic sound velocity , the difference practically disappears with the increase in the coupling strength. The general inequality was established by Schofield Schofield (1966). We see that from the side of soft long-ranged interactions, this inequality is very close to equality, the observation previously reported for soft IPL systems in 3D Khrapak et al. (2017). This tendency, however, breaks down in the case of extremely steep hard-sphere-like interactions, where , and are all diverging (in 2D and 3D  Khrapak et al. (2017); Khrapak (2016d)), whilst remains finite Rosenfeld (1999).

Vi Conclusion

We studied thermodynamics of two-dimensional IPL3 classical systems across coupling regimes, from the weakly non-ideal gas to the strongly coupled fluid and crystalline phases. Careful analysis of the extensive MD simulation results allowed us to put forward simple and physically suitable expressions for the thermodynamic properties (e.g. excess energy) of the investigated system. In particular, Helmholtz free energies of the fluid and solid phases have been derived and the location of the fluid-solid coexistence has been determined. The obtained results are comparable to those previously reported in the literature. A very narrow fluid-solid coexistence gap observed is likely related to the very soft nature of the interaction potential.

The QCA/QLCA approach has been applied to the description of collective modes of the IPL3 fluids. The use of a simplistic RDF has been suggested, based on previous results related to strongly coupled plasma fluids. This has allowed us to derive explicit analytic dispersion relations for the longitudinal and transverse modes, which have been checked against the results of direct MD simulations. Reasonable agreement in the long-wavelength regime has been observed. We also briefly pointed out that the obtained simple fluid dispersion relations can be of some use in estimating the harmonic entropy constant of the solid phase.

The expressions for various sound velocities have been examined. These include conventional longitudinal and transverse elastic sound velocities of the idealized IPL3 crystalline lattice, their analogues (based on a QCA/QLCA approximation) in the strongly coupled fluid state, as well as conventional adiabatic sound velocity of the IPL3 fluid. Additionally, expressions for the 2D high frequency (instantaneous) elastic moduli have been introduced and related to the sound velocities. One useful observation is that the instantaneous sound velocity (related to the instantaneous bulk modulus) is extremely close to the adiabatic sound velocity. This observation is very likely a general property of soft long-ranged potentials, related neither to the exact shape of the interaction potential, nor to the dimensionality.

Finally, we would like to point once more that the interaction potential studied in this work represents just one particular example of very soft long-ranged interactions. It is therefore important that the approaches used here can be directly (or with minor modifications) transferred and applied to other related soft-interacting particle systems.

This work has been supported by the A*MIDEX project (Nr. ANR-11-IDEX-0001-02) funded by the French Government “Investissements d’Avenir” program managed by the French National Research Agency (ANR). Numerical simulations and analysis have been supported by the Russian Science Foundation (RSF), Project No. 17-19-01691. Studies of fluid-solid transitions in 2D dusty plasmas and related model systems have been supported by RSF project No. 14-50-00124. We thank I. Semenov for careful reading of the manuscript.

Appendix A Thermodynamic relations used in this work

Here we express some of the reduced thermodynamic quantities of interest in terms of the reduced excess internal energy . For example, the compressibility is


The inverse reduced isothermal compressibility modulus is


The reduced isochoric heat capacity is


The adiabatic index is (for the considered potential and dimensionality)


The quantities and are used to calculate the conventional adiabatic sound velocity in IPL3 fluids.

Appendix B Internal energy at intermediate coupling

The dependence of the reduced excess energy on the coupling parameter obtained from MD simulations at weak and moderate coupling is shown in Fig. 5 along with the scalings at strong coupling (10) and weak coupling (12). Based largely on the combination of these two scalings, we constructed a practical approximation (interpolation) suitable for the intermediate coupling regime. The proposed interpolation is




is the smooth step function. Fitting MD data resulted in the following coefficients: , , , , and . The corresponding curve is also plotted in Fig. 5 documenting excellent agreement with the MD results.

Figure 5: The reduced excess energy, , of a weakly to moderately coupled IPL3 system in 2D versus the coupling parameter . Circles correspond to the results of MD simulations performed in this work, the green line corresponds to the weak coupling asymptote, Eq. (12), the orange line corresponds to the strong coupling asymptote of Eq. (10), and the red curve is the fit of Eq. (39) appropriate for the moderately coupled regime.

Appendix C Thermodynamics quantities of the IPL3 melt

In Table 2 we have tabulated some reduced thermodynamic quantities of the IPL3 fluid near the boundary of the fluid-solid coexistence (IPL3 melt) at . The quantities displayed are internal thermal energy (), isochoric heat capacity (), adiabatic index (), Helmholtz free energy (), excess internal energy (), excess entropy (), compressibility (), longitudinal elastic velocity (), and transverse elastic velocity ().

Table 2: Selected thermodynamic quantities (see the text for nomenclature) at the fluid boundary of the fluid-solid coexistence, at .

Appendix D Explicit expressions for the dispersion relations

Assume that a pairwise interaction potential can be written in the general form

where is the energy scale [the subscript is used here to demonstrate difference in energy scales compared to Eq. (1); for the IPL3 potential we have ]. The integrals in Eqs.  (14) and (15) can be simplified taking into account (in 2D), and that the derivatives of the interaction potential are


Integration over the angle is then performed with the help of the identities


where , and denote the Bessel functions of the first kind, related via

Introducing the reduced distance we obtain




Here is the nominal 2D frequency and is the reduced wave number. For the IPL3 potential with the dispersion relations become




To within some minor difference in notation, Eqs. (43) and (44) coincide with Eqs. (16) and (17) from Ref. Golden et al. (2010). These expressions were previously used to generate the dispersion curves with the input of data obtained in MD computer simulations Golden et al. (2010, 2009). A simplification, which does not require the accurate knowledge of the RDF, is discussed in Sec. III.2.

Appendix E Harmonic entropy constant from QCA dispersion relations

Taking into account that we have two (longitudinal and transverse) modes in a 2D lattice and approximating the first Brillouin zone by a disk with the area we can re-write the harmonic entropy constant [the last two terms in Eq. (27)] as


where correspond to the angularly averaged longitudinal and transverse phonon dispersion curves. As a simplest rough estimate one can approximate the phonon spectrum by its isotropic acoustic asymptote, Eqs. (16) and (29) as was done in Ref. Alastuey and Jancovici (1981) for a 2D OCP with logarithmic interactions. In this way we have obtained for the present case of IPL3 in 2D


which significantly overestimates the actual harmonic entropy constant [subscript “” in Eq. (46) denotes zero approximation]. In order to improve the accuracy we have also used the analytical QCA expressions (19) and (20) in Eq. (45). This approach is based on the observation that angularly averaged lattice dispersion relations show remarkable similarity to isotropic QCA dispersion relation, in particular within the first Brillouin zone Sullivan et al. (2006); Hartmann et al. (2007). Taking we have obtained in this approximation


which is considerably closer to the actual harmonic entropy constant. Thus, the fluid QCA dispersion relations in the strong coupling limit can be of some use in quickly estimating the free energy of the solid phase. Note that the harmonic entropy constant contributes only a small fraction of the total free energy for soft long-ranged interactions.

Appendix F Interpolation method for calculating RDF in 2D crystals

The anisotropic RDF of a crystal is written in the form Yurchenko et al. (2016b)


where the summation is over all the nodes , and each individual peak has the shape


The normalization constant as well as the parameters are defined by the conditions Yurchenko et al. (2016b)


where is the unit vector in the direction of , is the mean squared displacement for longitudinal and transverse directions, respectively.

The effect of the temperature dependence of phonon spectra can be taken into account by introduction of the anharmonic correction coefficient  Kryuchkov et al. (2017)


where the tildes denote the mean-squared displacement (MSD) calculated using BvK phonon spectra (see Ref. Yurchenko et al. (2015)), is the total MSD for the nearest neighbors.

Contrary to 3D crystals, in 2D cases the mean squared displacements diverge logarithmically. The resulting correlation peaks become less localized, so the overlap of the neighboring peaks is generally stronger for 2D crystals. Nevertheless, it turns out that the IM approach can be applied also in this case, in essentially the same way as for 3D crystals Yurchenko et al. (2016b).


Comments 0
Request Comment
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
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description