Thermodynamics and dynamics of two-dimensional systems with dipole-like repulsive interactions
Sergey Khrapak Nikita Kryuchkov Stanislav Yurchenko





Sergey.Khrapak@dlr.de
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


Bauman Moscow State Technical University, 2nd Baumanskaya street 5, 105005 Moscow, Russia


st.yurchenko@mail.ru
Bauman Moscow State Technical University, 2nd Baumanskaya street 5, 105005 Moscow, Russia





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. 


 Thermodynamics and dynamics of two-dimensional systems with dipole-like repulsive interactions
    Stanislav O. Yurchenko
    December 30, 2023
===============================================================================================





§ 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 <cit.>. 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 <cit.> and bubbles <cit.>  for synthesis of novel optical materials <cit.>, chemical sensors, and catalysis <cit.>, 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 <cit.>. As a result, the interaction between colloidal particles is similar to that between vertically oriented dipoles <cit.> and can be in the first approximation described by a pairwise repulsive inverse-power law (IPL) potential decaying as ∼ 1/r^3 with the interparticle separation r. Direct experimental measurements of colloidal interaction potential in such systems by the laser tweezers method <cit.> or using other approaches <cit.> 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 <cit.>. A similar shape of the interaction potential is observed in two-dimensional colloidal systems of paramagnetic particles exposed to an external magnetic field <cit.>. 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 attractive<cit.>  or repulsive <cit.>.

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 ∝ 1/r^3, at large
distances from the test charge <cit.>.
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 <cit.>, there is a certain parameter regime, where the IPL scaling ∝ 1/r^3 is relevant <cit.> (see, in particular, Fig. 3 in Ref. <cit.>).

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 ∝ 1/r^3 repulsion have been extensively studied (see, for instance, Refs. <cit.> 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 <cit.> 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 <cit.>. 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 <cit.> with main applications to complex (dusty) plasmas. Here we apply the same arguments to the 2D system with 1/r^3 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) <cit.>, also referred to as the quasi-localized charge approximation (QLCA) in the plasma-related context <cit.>, 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 <cit.>). Previously, QLCA has been applied with certain success to dipole-like systems in 2D <cit.>. 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  <ref> 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 <ref> main results obtained for the fluid phase are summarized, including accurate expressions for thermodynamic quantities and detailed analysis of collective modes. In Section <ref> 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 <ref> 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 <ref>.




§ METHODS





 §.§ 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

    ϕ(r) = ϵ(σ/r)^3,

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

    Γ=ϵ/T(σ/a)^3,

where T is the temperature (in energy units), a=(πρ)^-1/2 is the 2D Wigner-Seitz radius, and  ρ=N/V is the areal density of N particles occupying the 2D volume (i.e. surface) V. 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 Γ≫ 1.

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 <cit.>. 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 <cit.>. For softer interactions the liquid-hexatic transition is continuous, with correlations consistent with the Berezinsky-Kosterlitz-Thouless-Halperin-Nelson-Young (BKTHNY) scenario <cit.>. For the IPL family of potentials (∝ r^-n) the transition between the two regimes occurs at about n≃ 6 <cit.>. 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 <cit.> and colloidal experiments <cit.>.
However, the hexatic phase occupies a relatively narrow region on the phase diagram and its properties will not be addressed in this work.




 §.§ 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 <cit.>. These MD simulations have been done in the NVT ensemble at different temperatures using N = 4 × 10^4 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 1.5 T and 2 T for the fluid and crystal phases, respectively. The numerical time step of Δ t = 2.4 × 10^-4√(ma^2Γ/ϵ) was chosen. All simulation runs were performed for 6× 10^5 time steps, where the last 4× 10^5 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 25ρ^-1/2. The internal energy and pressure were corrected accordingly, which resulted in a relative error of about 2× 10^-5, sufficient for the range of problems considered here.



 §.§ Thermodynamic relations


The main thermodynamic quantities of interest in this work are the internal energy U, Helmholtz free energy F, and pressure P of the system. The following thermodynamic definitions are useful <cit.>

    U=-T^2(∂/∂ TF/T)_V = ∂(F/T)/∂( 1/T)|_V, 
    
    P=-(∂ F/∂ V)_T.

In addition, U and P can be calculated using the integral equations of state <cit.>

    U= N(T+ ρ/2∫d𝐫 ϕ(r)g( r)),
        PV = N(T - ρ/4∫d𝐫  rϕ'(r)g( r)),

where g(𝐫) 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: u=U/NT, f=F/NT, and p=PV/NT and divide the thermodynamic quantities into the kinetic (ideal gas) and potential (excess) contributions, so that u=1 + u_ ex (in 2D), f=f_ id+f_ ex, and p=1+p_ ex. Finally, it is useful to operate with the single coupling parameter Γ, instead of temperature and density. Since Γ∝ a^-3T^-1∝ρ^3/2T^-1, the transformation of standard thermodynamic relations to their dimensionless form is governed by

    ∂ X/∂ρ=3Γ/2ρ∂ X/∂Γ, 
       ∂ X/∂ T=-Γ/T∂ X/∂Γ,

where X 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:

    p_ ex=3/2u_ ex.

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 <ref>.



§ FLUIDS
 



 §.§ 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 ∝ 1/r^3 potential (this relation is meaningful for both crystalline and fluid phases). The corresponding lattice sum has been evaluated previously with a very high accuracy <cit.>. The proportionality constant between the static energy and the coupling parameter (Madelung constant) is

M≃0.798512.

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

    u_ fl = MΓ + u_ th.

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 g(r)≃ 1. 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 (MΓ≫ 1), 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 (u_ th∼ 1). This implies that even moderately accurate approximations for u_ th result in a very accurate estimation of the total excess energy u_ fl of strongly coupled fluids. The remaining step is therefore to identify an appropriate approximation for u_ th.

Based on the previous results for other soft interacting particle systems in 2D (such as one-component-plasma <cit.> and weakly screened Yukawa systems <cit.>), 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 <cit.>, but has a quite different functional form. The functional form suggested for 2D systems with soft pairwise repulsive interactions is

    u_ th =a ln (1+bΓ).

The validity of this functional form at sufficiently strong coupling is documented in Fig. <ref>, where numerical data from the present MD simulations along with those reported previously <cit.> are plotted. The best fit of the MD data obtained in this work yields a=0.27284 and b=2.2357. The fit is valid in the range 1≲Γ≲ 70. Combining Eqs. (<ref>) and (<ref>) we write for the excess energy of the strongly coupled fluid phase

    u_ fl=MΓ+a ln (1+bΓ).

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. (<ref>) is ≃ 0.4 at Γ=1, it decreases to ≃ 0.1 at Γ = 10, and further drops to ≃ 0.03 at Γ=50. 




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

    f_ fl = ∫_0^ΓdΓ' u_fl(κ, Γ')/Γ'.

However, one must pay some attention to the fact that expression (<ref>) is not applicable all the way down to Γ=0. 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 <cit.>

    f_1=ρ/2∫[1-e^-ϕ(r)/T]d r=∫_0^∞[1-e^-Γ/x^3]xdx,

which can be evaluated analytically. 




The excess energy in this regime is u_ ex=23f_1. However, careful comparison with the results from numerical simulations shows that Eq. (<ref>) is applicable only at very weak coupling, Γ≲ 0.05. On the other hand, the strong coupling scaling (<ref>) is justified only for Γ≳ 1. An approximation for the intermediate regime, basically based on an appropriate combination of Eqs. (<ref>) and (<ref>) is described in Appendix <ref>. Using this approximation [Eq. (<ref>)] below Γ=10 and Eq. (<ref>) at higher values of Γ, and performing integration in Eq. (<ref>) we finally obtained a simple and accurate analytical approximation for the fluid free energy in the strong coupling regime, Γ≳ 10:

    f_ fl=MΓ-a Li_2(-bΓ)+0.381,

where Li_2(x) 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 (<ref>) and (<ref>) 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 Γ≃ 69, with a very narrow coexistence gap (see Section <ref>). 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 <ref>.




 §.§ Collective modes


It is well known that fluids can exhibit different collective dynamics depending on the regime of coupling and correlations <cit.>.
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 <cit.>. 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 <cit.>. In recent years it was successively applied to strongly coupled one-component plasma in both 2D and 3D <cit.> as well as 2D and 3D Yukawa fluids, mostly in the context of complex (dusty) plasmas <cit.>. Application to the 2D IPL3 system was described in Refs. <cit.>. 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), g(r), and the pair interaction potential ϕ(r). The compact expressions are <cit.>

    ω_ L^2=ρ/m∫∂^2 ϕ(r)/∂ z^2 g(r) [1-cos(kz)]d r,


    ω_ T^2=ρ/m∫∂^2 ϕ(r)/∂ y^2 g(r) [1-cos(kz)]d r,

where ω is the frequency, k is the wave number, and z=rcosθ is the direction of the propagation of the longitudinal mode (the particles are confined to the zy plane). Here the subscripts “L” and “T” 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 <ref> for completeness. 



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

    lim_k→ 0ω_ L/T^2/k^2 = C_ L/T^2.

Similarly to the IPL system in the 3D case <cit.>, 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:

    C_ L^2=33/8v_ T^2 u_ ex,      C_ T^2=3/8v_ T^2 u_ ex,

where v_T=√(T/m) is the thermal velocity of the particles.
Here we used the relation that follows directly from energy or pressure equations (<ref>),

p_ex=3/2u_ex=3/4Ω_0^2a^2/v_T^2∫_0^∞g(x)dx/x^2,

where  Ω_0^2= 2πρϵσ^3 /ma^3 is the conventional 2D frequency scale.  Note an immediate consequence of Eq. (<ref>), C_ L/C_ T=√(11), for the IPL3 system in 2D.

A simplest model g(r), 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

    g(x)=θ(x-R),

where θ(x) is the Heaviside step function and the radius of the correlational hole R is of order unity in reduced units (the distances are expressed in units of a 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 <cit.>. It was also used to calculate the dispersion relation of Coulomb bilayers and superlattices at strong coupling <cit.>. Physically, the model form of Eq. (<ref>) seems sensible in the present context, because the main contribution to the long-wavelength dispersion corresponds to long
length-scales, where g(x)≃ 1. For soft enough interactions, this regime provides dominant contribution to the integrals in Eqs. (<ref>) and (<ref>). The excluded volume effect for x≤ R 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 <cit.> and one-component plasma with logarithmic interactions in 2D <cit.>.


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

    ω_ L^2=Ω_0^2/R^3{3/2-q^3R^3+J_1(qR)/2qR[6+2q^2R^2-2q^4R^4+π q^5 R^5 H_0(qR)
    ]
    -J_0(qR)/2[
    6-2q^2R^2-2q^4R^4+π q^4 R^4 H_1(qR)
    ]
    },

where q=ka is the reduced wave number, J_0(x) and J_1(x) are Bessel functions of the first kind, H_0(x) and H_1(x) are the Struve functions of order 0 and 1, respectively. The dispersion relation of the transverse mode is remarkably more simple,

    ω_ T^2=Ω_0^2/R^3{3/2-3J_1(qR)/qR}.


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

    R=Γ/u_ ex.

It is straightforward to verify that with this definition of R, the low-q series expansion of Eqs. (<ref>) and (<ref>) will reproduce the acoustic velocities given by (<ref>). Note also that in the strong coupling regime the excess energy is mainly associated with the static contribution, u_ ex≃ MΓ. In this regime the radius of the correlational hole is practically constant, R≃ 1/M≃ 1.2523.



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 <cit.>.
Specifically, we calculated

    𝒜_ L, T(𝐪,ω) ∝ q Re∫ dt <j_ L, T(𝐪,t)j_ L, T(-𝐪,0)> e^iω t,

where j_ L(𝐪,ω) and j_ T(𝐪,ω) are the projections of velocity current 𝐣(𝐪,t)∝∑_i 𝐯_i(t) exp[i 𝐪𝐫_i(t)] to longitudinal and transversal directions, respectively.
Here the summation is performed over all particles in the system,
𝐫_i(t) is radius-vector of the i-th particle and 𝐯_i(t) is its velocity. Since fluids are isotropic we can average
𝒜_ L, T(𝐪,ω) over all directions of the wave vector 𝐪 to get the dependence on q=|𝐪|.

Figure <ref> 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 Γ=56 (a) - (c) and Γ=28 (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 𝒜_ L, T(q,ω) by the Cauchy distribution,

    f(ω)∝1/(ω-ω_0)^2+α^2+1/(ω+ω_0)^2+α^2,

for each value of q. Examples of the obtained dependencies ω_0(q) and α(q) are plotted in Fig. <ref>. 

The long-wavelength portions of the dispersion relations ω_0(q) obtained from MD simulation are plotted in Figs. <ref>(c) and <ref>(f). Here they are compared with QCA dispersion relations. The solid curves correspond to the “full” QCA with the actual RDF g(r) obtained in MD simulations and substituted in Eqs. (<ref>) and (<ref>). The black dashed curves correspond to the “simplified” QCA given by analytical expressions (<ref>) and (<ref>). The agreement between the two versions of QCA and MD dispersions is satisfactory at sufficiently long wavelengths (q≲ 2 for the longitudinal and q≲ 3 for the transverse mode). This (low-q) regime corresponds to long length-scales, where both actual and model RDF are similar, g(r)≃ 1.  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 g(r). 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 g(r). 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 <cit.>. At the same time, anharmonicity is related to the short-range region of the interaction potential, which corresponds to large q 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 q-gap for the transverse waves propagation <cit.> 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. <ref>(b) and (e) that the width of the q-gap decreases with increasing correlations (i.e. increasing Γ), in agreement with previous studies on collective excitations in various kinds of liquids <cit.>.    

Regarding the longitudinal mode, the asymptotic behavior of the dispersion relation at q≫ 1 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-q asymptote ω^2 ∝ k^2 v_ T^2 is consistent with MD simulation results. The proportionality coefficient 9/4 has been previously suggested in Ref. <cit.>. Our preliminary analysis indicates that the coefficient 2 can be more appropriate. The small (∼ 10%) 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.   



§ CRYSTALS


The reduced excess energy of a 2D crystalline lattice in the harmonic approximation is u_ h=MΓ +1.
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 <ref>. The anharmonic corrections are fitted using the standard functional form <cit.>

    u_ anh=A_1/Γ+A_2/Γ^2+A_3/Γ^3.

The coefficients determined from the fitting are A_1=2.47672, A_2=-148.77, and A_3=13507.4 (curve in Fig. <ref>). The resulting excess internal energy of the solid phase is

    u_ s=MΓ+1+A_1/Γ+A_2/Γ^2+A_3/Γ^3.


The reduced excess Helmholtz free energy can then be evaluated in the following manner <cit.>. First, we divide it into anharmonic and harmonic contributions

    f_ s = ∫_∞^ΓdΓ' u_anh(κ, Γ')/Γ'+f_ h,

where f_ h 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 <cit.>. In 2D geometry the resulting expression for the reduced harmonic free energy is <cit.>

    f_ h= MΓ + lnΓ + 1+ 1/2N∑_ k, slnω_ s^2( k)/Ω_0^2,

where ω_ s( k) is the frequency of a phonon with wavenumber k and polarization s, and the sum on k 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. <ref>. The corresponding harmonic entropy constant has been evaluated as Σ = 0.09284 (a related approach to estimate Σ using the QCA dispersion relations is briefly discussed in Appendix <ref>). Thus, the reduced Helmholz free energy of the IPL3 crystalline solid is

    f_ s = MΓ + lnΓ+0.09284-A_1/Γ-A_2/2Γ^2-A_3/3Γ^3.

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 Γ_ m≃ 69 (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 Γ_ L≃ 69.2 and Γ_ S≃ 69.4. This is comparable to the results previously reported in the literature <cit.> and is particularly close to the results from the Brownian dynamics simulations <cit.>. Note that a very narrow coexistence gap obtained, ΔΓ/Γ_ S≃ 0.003 (here ΔΓ = Γ_ S-Γ_ L), 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 Γ_ m≃ 75. As a final remark, we note that we have not considered the existence of the hexatic phase. 



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 <cit.>. 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 <ref>.

In Fig. <ref> an example of the crystalline RDF g(r) is presented. Visual similarity between the obtained theoretical RDF and MD data is the same as in previous applications of the IM approach <cit.>. Using these highly accurate RDFs, pressure and excess energy can be readily obtained using Eqs. (<ref>). In its simplest harmonic form (β=0) the IM approach provides a relative error in the excess energy smaller than ≃ 10^-3 in the worst case near the melting point. Taking into account anharmonicity, with the anharmonic correction coefficient β=7.54 [see Eq. (<ref>)], obtained from MD simulations of the IPL3 crystal, reduces the relative error to ≃ 10^-4.

The longitudinal and transverse sound velocities of a perfect IPL3 2D lattice are given by Eqs. (<ref>) combined with u_ ex=MΓ (we remind that QCA reduces to the standard phonon theory of solids in the limit T=0). The final result is

    C_ L≃ 1.8149 (ϵ/m)^1/2(σ/a)^3/2, 
        C_ T≃ 0.5472 (ϵ/m)^1/2(σ/a)^3/2.







§ 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 <cit.>. The 2D analogues of these elastic moduli are

    K_∞=2ρ T-πρ^2/4∫_0^∞drr^2g(r)[ϕ'(r)-rϕ”(r)],

the high frequency limit of the bulk modulus <cit.>, and

    G_∞=ρ T+πρ^2/8∫_0^∞drr^2g(r)[3ϕ'(r)+rϕ”(r)],

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

    mρ C_ L^2= K_∞+G_∞-3ρ T,     
    mρ C_ T^2= G_∞ -ρ T.

Useful relations between C_ L and C_ T in both 3D and 2D have been recently discussed <cit.>. In particular, the high-frequency (instantaneous) sound velocity C_∞, directly related to the instantaneous bulk modulus, can be introduced

    C^2_∞=K_∞/mρ = 2v_ T^2+C_ L^2-C_ T^2=v_ T^2(2+15/4u_ ex).

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 <cit.>

    C_ s =1/m(∂ P/∂ρ)_S = v_ T√(γμ).






The sound velocities C_ L,  C_∞, and C_ s for several values of the coupling parameter Γ, corresponding to the strongly coupled fluid phase, are summarized in Table <ref>. It is observed that C_ L overestimates the adiabatic sound velocity C_ s. However, the difference is rather small, as should be expected for the soft interaction potential studied in this work <cit.>.
(For soft interactions at strong coupling one normally observe C_ L≫ v_ T and C_ L≫ C_ T, which implies C_ L∼ C_∞≃ C_ s <cit.>). The instantaneous sound velocity C_∞ is just slightly above the adiabatic sound velocity C_ s, the difference practically disappears with the increase in the coupling strength. The general inequality C_ s≤ C_∞ was established by Schofield <cit.>. 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 <cit.>. This tendency, however, breaks down in the case of extremely steep hard-sphere-like interactions, where C_ L, and C_∞ are all diverging (in 2D and 3D  <cit.>), whilst C_ s remains finite <cit.>.



§ 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.






§ 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 u_ ex. For example, the compressibility Z=PV/NT is

    Z=1+p_ ex=1+3/2u_ ex.

The inverse reduced isothermal compressibility modulus μ=(1/T)(∂ P/∂ρ)_T is

    μ = 1+3/2u_ ex+9Γ/4∂ u_ ex/∂Γ.

The reduced isochoric heat capacity c_ V=(1/N)(∂ U/∂ T)_V is

    c_ V=1+u_ ex-Γ∂ u_ ex/∂Γ.

The adiabatic index γ = c_ P/c_ V is (for the considered potential and dimensionality)

    γ = 1+(3c_ V-1)^2/4μ c_ V.

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



§ INTERNAL ENERGY AT INTERMEDIATE COUPLING


The dependence of the reduced excess energy u_ex on the coupling parameter Γ obtained from MD simulations at weak and moderate coupling is shown in Fig. <ref> along with the scalings at strong coupling (<ref>) and weak coupling (<ref>). 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

    u_ ex =2/3[1-ξ(Γ)] f_1+ ξ(Γ)[MΓ+Aln(1+B Γ^s)],

where

    ξ(Γ)=[1+e^-C(Γ-Γ_0)]^-1.

is the smooth step function.
Fitting MD data resulted in the following coefficients: A=0.4791, B=1.2198, s=0.6044, C=428.216, and Γ_0=2.25× 10^-2. The corresponding curve is also plotted in Fig. <ref> documenting excellent agreement with the MD results. 






§ THERMODYNAMICS QUANTITIES OF THE IPL3 MELT


In Table <ref> we have tabulated some reduced thermodynamic quantities of the IPL3 fluid near the boundary of the fluid-solid coexistence (IPL3 melt) at Γ = 69. The quantities displayed are internal thermal energy (u_ th), isochoric heat capacity (c_ V), adiabatic index (γ= c_ P/c_ V), Helmholtz free energy (f_ ex), excess internal energy (u_ ex), excess entropy (s_ ex=u_ ex-f_ ex), compressibility (Z=1+p_ ex), longitudinal elastic velocity (C_ L/v_ T), and transverse elastic velocity (C_ T/v_ T).






§ EXPLICIT EXPRESSIONS FOR THE DISPERSION RELATIONS


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

ϕ(r)=ϵ_0 f(r/a),

where ϵ_0 is the energy scale [the subscript is used here to demonstrate difference in energy scales compared to Eq. (<ref>); for the IPL3 potential we have ϵ_0=ϵ(σ/a)^3]. The integrals
in Eqs.  (<ref>) and (<ref>) can be simplified taking into account d r= rdr dθ (in 2D), kz=krcosθ and that the derivatives of the interaction potential are

∂^2 ϕ/∂z^2=ϕ”(r)z^2/r^2+1/rϕ'(r)(1-z^2/r^2),
    
and

∂^2 ϕ/∂y^2=Δϕ(r) - ∂^2 ϕ/∂z^2.

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

1/2π∫_0^2π [1-cos(krcosθ)]cos^2θd θ=1/2-J_1(kr)/kr+J_2(kr)

and 

1/2π∫_0^2π [1-cos(krcosθ)]d θ=1- J_0(kr),

where J_0(x), J_1(x) and J_2(x) denote the Bessel functions of the first kind, related via

J_0(x)+J_2(x)=2J_1(x)/x. 

Introducing the reduced distance x=r/a we obtain


    ω_ L^2= Ω_0^2∫_0^∞ g(x)dx{f'(x)[1/2-J_1(qx)/qx]+xf”(x)[1/2+J_1(qx)/qx-J_0(qx)]
    },

and

    ω_ T^2= Ω_0^2∫_0^∞ g(x)dx{f'(x)[1/2+J_1(qx)/qx-J_0(qx)]+xf”(x)[1/2-J_1(qx)/qx]}.


Here Ω_0^2=2πρϵ_0 /m= 2πρϵσ^3 /ma^3 is the nominal 2D frequency and  q=ka is the reduced wave number. For the IPL3 potential with f(x)=1/x^3 the dispersion relations become   

    ω_ L^2= 3Ω_0^2/2∫_0^∞g(x)dx/x^4[3-3J_0(qx)+5J_2(qx)],

and

    ω_ T^2= 3Ω_0^2/2∫_0^∞g(x)dx/x^4[3-3J_0(qx)-5J_2(qx)].

To within some minor difference in notation, Eqs. (<ref>) and (<ref>) coincide with Eqs. (16) and (17) from Ref. <cit.>. These expressions were previously used to generate the dispersion curves with the input of g(r) data obtained in MD computer simulations <cit.>. A simplification, which does not require the accurate knowledge of the RDF, is discussed in Sec. <ref>.



§ 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 4π^2ρ we can re-write the harmonic entropy constant  [the last two terms in Eq. (<ref>)] as

    Σ=1+1/4∫_0^2[lnω_ L^2(q)/Ω_0^2
    +lnω_ T^2(q)/Ω_0^2]qdq,

where ω_ L,T(q) 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. (<ref>) and (<ref>)  as was done in Ref. <cit.> for a 2D OCP with logarithmic interactions. In this way we have obtained for the present case of IPL3 in 2D

    Σ_0≃ 0.6862,

which significantly overestimates the actual harmonic entropy constant [subscript “0” in Eq. (<ref>) denotes zero approximation]. In order to improve the accuracy we have also used the analytical QCA expressions (<ref>) and (<ref>) in Eq. (<ref>). 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 <cit.>.
Taking R=1/M=1.25233 we have obtained in this approximation

    Σ_ QCA≃ 0.1336,

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.



§ INTERPOLATION METHOD FOR CALCULATING RDF IN 2D CRYSTALS


The anisotropic RDF g(𝐫) of a crystal is written in the form <cit.>

    g(𝐫) = 1/ρ∑_αp_α(𝐫-𝐫_α),

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

    p_α(𝐫) ∝exp[-ϕ(𝐫+𝐫_α)/k_ BT-b_α(𝐞_α·𝐫)-
    . 
                       .
    -(𝐞_α·𝐫)^2/2 a_α^2-
    𝐫^2-(𝐞_α·𝐫)^2/2 a_⊥α^2].

The normalization constant as well as the parameters a_,⊥α^2, b_α are defined by the conditions <cit.>

    ∫d𝐫 p_α(𝐫)=1,     ∫d𝐫 𝐫p_α(𝐫)=0, 
       ∫d𝐫 (𝐞_α·𝐫)^2 p_α(𝐫)=σ_α^2,
       ∫d𝐫 [𝐫^2-(𝐞_α·𝐫)^2] p_α(𝐫)= σ_⊥α^2,

where 𝐞_α=𝐫_α/r_α is the unit vector in the direction of 𝐫_α,
σ_,⊥^2 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 β <cit.>

    σ_,⊥α^2 = σ_,⊥α^2 [1+β Nσ_1^2/V],

where the tildes denote the mean-squared displacement (MSD) calculated using BvK phonon spectra (see Ref. <cit.>),
σ_1^2 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 <cit.>.




