Collective modes in simple melts: Transition from soft spheres to the hard sphere limit

We study collective modes in a classical system of particles with repulsive inverse-power-law (IPL) interactions in the fluid phase, near the fluid-solid coexistence (IPL melts). The IPL exponent is varied from n = 10 to n = 100 to mimic the transition from moderately soft to hard-sphere-like interactions. We compare the longitudinal dispersion relations obtained using molecular dynamic (MD) simulations with those calculated using the quasi-crystalline approximation (QCA) and find that this simple theoretical approach becomes grossly inaccurate for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\,\gtrsim 20$$\end{document}n≳20. Similarly, conventional expressions for high-frequency (instantaneous) elastic moduli, predicting their divergence as n increases, are meaningless in this regime. Relations of the longitudinal and transverse elastic velocities of the QCA model to the adiabatic sound velocity, measured in MD simulations, are discussed for the regime where QCA is applicable. Two potentially useful freezing indicators for classical particle systems with steep repulsive interactions are discussed.

The problem of describing collective motion in liquids is one of the long standing major topics in condensed matter physics [1][2][3][4] . One of the theoretical approaches to collective motion in neutral liquids is, in certain sense, analogous to that applied to elastic waves in disordered solids. Comparable expressions for the dispersion relations of these elastic waves can be obtained from the analysis of forth and second frequency moments of the dynamical structure factor S(k, ω) 3,5 , (that is from the sum rules) or from the variational calculation of the frequency spectra of phonons in classical fluids, performed by Zwanzig 6 . Particularly enlightening derivation is due to Hubbard and Beeby 7 who considered their simple theory as either a generalization of the random phase approximation or, alternatively, as a generalization of the phonon theory of solids. Taking this latter viewpoint, the approach was also referred to as the quasi-crystalline approximation (QCA) 8 and we keep this notation for the class of considered theories throughout the paper. In many cases the QCA-based description is in fair agreement with the dispersion curves deduced experimentally for various simple liquids, at least in the vicinity of the first frequency maximum [7][8][9][10][11] .
In the context of plasma physics, an analogue of the QCA, known as the quasi-localized charge approximation (QLCA), has been initially proposed as a formalism to describe collective mode dispersion in strongly coupled charged Coulomb liquids 12 . In the last couple of decades the QCA/QLCA approach has been numerously applied to collective wave spectra in one-component-plasma (OCP) in both two-dimensions (2D) and three-dimensions (3D) [12][13][14][15] , magnetized 3D OCP 16,17 , 2D and 3D Yukawa systems in the context of complex (dusty) plasmas [18][19][20][21][22][23] , 2D dipole systems 24 , etc. In general, the QCA/QLCA model has been documented to describe rather good collective oscillation spectra for soft repulsive interactions between charged particles in the regime of strong coupling (when the potential energy of interaction is much greater than the particle kinetic energy). Interestingly, being developed for the strongly coupled liquid state, QCA reduces to the conventional phonon-dispersion relation in the special case of cold crystalline solid and can describe real part of the dispersion of the electron (Langmuir) and ion (ion acoustic) waves in classical weakly-coupled (gaseous) plasmas.
Thus, QCA and (using present notation) its variations have been applied to various real and model systems in a wide interdisciplinary context. It is the appealing simplicity combined with reasonable accuracy, which made QCA a commonly used approximation in studies related to condensed matter, soft matter, and plasma physics.
Within the QCA the dispersion relations are explicitly given in terms of the radial distribution function and the pair interaction potential by means of very compact expressions [see Eqs. (2) and (3) below]. More involved theories are available, but they require more information about systems under investigation. For example, the mode coupling theory (MCT) 25 can also successively describe collective excitations of simple classical liquids [26][27][28] . However, compared to QCA, it requires additionally information about three-particle correlations and the calculation procedure is much more complex 27 . Although MCT demonstrates some improvement over QCA, in particular in the vicinity of the first (roton) minimum of the longitudinal dispersion curve 28 , QCA can still be considered as a very useful zero approximation to the dispersion relation.
Some deficiencies of the QCA approach (in the original Hubbard and Beeby formulation) are well recognized. For example, it does not include direct thermal effects associated with the free movement of the particles in a liquid (although this is not expected to be an important issue at sufficiently strong coupling). Damping effects are also excluded in the original theory (in complex plasmas, for instance, the collisions between charged particles and neutral atoms or molecules represent an important damping mechanism, which has to be accounted for). Possible modifications to QCA in order to remove these deficiencies have been discussed in the literature 18,22,29 , and demonstrated to be helpful in certain situations. Overall, it is well documented that the QCA model is a reliable and accurate tool to describe collective modes for strongly interacting particles, provided the interaction potential is soft enough (e.g., in plasma-related systems).
An important issue is, therefore, related to the applicability limits of QCA from the side of steep (hard-sphere-like) interactions. Recently, it has been observed that the sound velocity in a system of particles interacting via the inverse-power-law (IPL) repulsive potential ∝ − V r r ( ) n , calculated using the QCA approach, diverges as ∝ n as the potential exponent n increases 30 . This tendency contradicts the finite value of the sound velocity in a system of hard spheres 31 . Such contradiction is a strong indication that QCA, being a very good approach for soft interactions, loses its accuracy and eventually applicability as the softness of the interaction potential decreases. The purpose of the present paper is to provide an approximate location of the applicability limits of QCA-based approaches.
The most direct way is chosen in this study. We consider the IPL family of potentials, n where ε is the energy scale and σ is the length scale, and perform extensive molecular dynamics (MD) simulations of the static and dynamical properties of the IPL fluids by increasing the exponent n from n = 10 to n = 100. The phase state of the studied IPL systems corresponds to a dense fluid, just at the boundary of the fluid-solid coexistence region (hence we call it melt throughout the paper). One of the output of the simulations is the calculated current fluctuation spectra, which contains information about the dispersion of collective modes. We then compare the dispersion curves from MD simulations with those calculated with the help of QCA approach and identify the condition when significant deviations emerge. As n increases, the steepness of the potential (1) increases and in the limit → ∞ n it approaches the interaction of hard spheres (HS) of diameter σ. Thus, the chosen strategy seems appropriate for the main purpose of the present study. Note that collective excitations in the IPL fluid with n = 12 have been recently investigated 32 from a somewhat different perspective.

Results
The equilibrium radial distribution function g(r) is required as an input to calculate the dispersion relations of collective modes within the QCA approach (see Methods for details). Thus, we briefly address some important structural properties first. Figure 1 shows the radial distribution functions, g(r), obtained from numerical simulations. The distances are expressed in units of the Wigner-Seitz radius, πρ = − a (4 /3) 1/3 , throughout the paper. We observe that the long-range part is very little affected by the variation of the potential steepness. On the other hand, the shape of the first coordination shell [region around the first maximum of g(r)] is rather sensitive to the potential steepness. In particular, we see that the value of g(r) at the first maximum increases considerably with n.
The familiar Raveché-Mountain-Streett freezing rule 33 postulates that the ratio of the first non-zero minimum to the first maximum of g(r), that is = R g g / min max , remains constant at crystallization. This rule was, however, put forward to describe freezing of the Lennard-Jones fluid and it is not very useful for the IPL family of potential studied here. We see in the inset in Fig. 1 that the ratio R drops by almost a factor of two when n increases from 10 to 100. Similar observation was previously reported 34 . We also see, however, that the value of g(r) at the first minimum remains practically constant, . ± .  g 0 55 0 02 min . Previously, the value of g min has been used to characterize the fluid-solid (freezing) phase transition in the Lennard-Jones system 35,36 , and the glass transition in metallic alloys 37 . To which extent this simple freezing indicator applies to other simple model fluids should be investigated in future work.
The dispersion relations of the longitudinal modes in IPL melts are shown in Fig. 2. The color background corresponds to the spectral decomposition of the longitudinal current fluctuations. The maximum magnitude (red color) marks the location of dispersion curves as measured in the MD numerical experiment. The black dashed curves correspond to the calculations based on the QCA model with the radial distribution functions taken from MD simulations. Reasonable agreement between QCA results and the current fluctuations analysis is observed in the vicinity of the first frequency maximum for sufficiently soft interactions with = n 10 (a), .  n 14 3 (b), and .  n 16 6 (c). For steeper potentials the deviations are very well observable for = n 20 (d) and  n 33 (e). Finally, for a very steep interaction with n = 50 (f) the QCA curve is completely off MD simulation results.
In Fig. 3 we compare the sound velocities obtained using different approaches. The black diamond, corresponding to = ∞ n , is the adiabatic sound velocity of the hard sphere system at the fluid side of the fluid-solid coexistence 31 . The green symbols (c s ) correspond to the linear fit of the MD data on current fluctuation spectra in . All the sound velocities are expressed in units of the thermal velocity v T . It is observed that c s and ∞ c practically coincide for  n 25, while c L overestimates these sound velocities by approximately a factor of 3/ 5 in the considered regime. These observations are discussed in detail in the next Section. The blue symbols (left axis) correspond to the familiar Raveché-Mountain-Streett ratio g g / min max at freezing 33 . The red symbols correspond to the value of the RDF at the first non-zero minimum, g min . We observe that the later quantity is contained in a very narrow range at freezing of the IPL systems.

Discussion
One of the main observations is that the QCA model, designed to describe elastic modes in simple fluids, has limited applicability with respect to the interactions softness. It is only applicable for sufficiently soft interactions.
In terms of the IPL family studied here, the potential exponent n should satisfy  n 20. Even when the QCA reasonably describes the dispersion near the first frequency maximum, the predicted elastic sound velocity C L is generally higher than the adiabatic sound velocity C s (see Fig. 3). The magnitude of C s can still be estimated using the QCA by means of the instantaneous sound velocity , related to the instantaneous bulk modulus (see section Methods for details). For sufficiently large exponent n, the ratio of elastic to instantaneous sound velocities is approximately 3/ 5. When the potential becomes softer this ratio decreases and reach unity for n = 3. The case n = 3 is special for the IPL model since for ≤ n 3 a uniform neutralizing background should be applied 38 . In addition, the dispersion relations for longitudinal waves are non-acoustic in the long-wavelength limit for n = 1 and n =2 12, 39 . A general remark regarding the regime of soft interactions is appropriate here. In soft interacting dense fluids, not too far from the fluid-solid transition, the longitudinal sound velocity is normally much higher than the transverse sound velocity,  C C L T . This implies that ∞  C C L and QCA provides a rather good estimate of the actual thermodynamic sound velocity. The very fact that the QCA elastic and thermodynamic sound velocities are rather close to each other has been numerously documented for weakly screened Yukawa (screened Coulomb) systems at strong coupling [40][41][42] .
It is observed in Fig. 3 that there is a very wide range of potential softness, at least from n = 10 to = ∞ n , where the longitudinal sound velocity near freezing varies weakly when approaching its HS asymptotic value ( . 12 5). This is unlikely to be a specific property of the IPL system. For example, it is well recognized that the ratio of the sound to thermal velocity of many liquid metals and metalloids has about the same value 10 at the melting temperature 43,44 . This observation has been already discussed in the context of the HS system 31 and, more recently, in the context of soft-sphere interactions 30 . We just remark here that the approximate constancy of the sound velocity (in units of thermal velocity) near the fluid-solid phase transition can perhaps be of some use as an approximate dynamical freezing indicator for systems with steep repulsive interactions.
In an attempt to understand the failure of the QCA consideration in the limit of steep interactions we evaluated some additional dynamical properties of the system. The results are summarized in Figs 4-6. Figure 4 demonstrates how the velocity autocorrelation functions (VAF) vary upon increasing the potential steepness. The shape of VAFs is changed considerably as n increases. For ∼ n 10 a clear "caging" behaviour of particles trapped and oscillating in slowly fluctuating potential wells, reminiscent particle dynamics in OCP 45 , is observed. For ∼ n 100 these oscillations are heavily suppressed and the shape of the VAF tends to that of the HS fluids near the freezing point 46 . The amplitude of the first minimum of the VAF shifts considerably upwards with the increase of n (see inset). Thus, a pronounced oscillatory caging behaviour seems important for the application of the QCA, in qualitative agreement with the physical picture behind the development of the QLCA model 12 .
Related to this is the observation by Brazhkin et al. 47,48 , that the caging dynamics changes significantly upon increasing the potential steepness of the IPL potential. Namely, they calculated the potential-to-kinetic energy ratio of the IPL fluids near the melting curve and found that it drops below unity at ∼ n 30. We repeated this calculation using the data summarized by Agrawal and Kofke 34 and present the obtained results in the inset of Fig. 6. Based on the observations, the particle motion inside the cages can be characterized as harmonic in the soft interaction regime and purely collisional absolutely non-harmonic in the steep interaction regime 47 . The transition between these two regimes at ∼ n 30 is close to the point where QCA failure occurs. The dynamical crossover is likely to be the reason. In Fig. 5 individual particle trajectories in a small cubic box, placed near the center of the computation cell, are shown for two cases, n = 10 and n = 50. The trajectories are color coded according to the particles energy. Trajectories shown in (a) appear relatively smooth and continuous. In contrast, we observe in (b) strong breaks in particles trajectories for steeper potential with n = 50, which result from short hard-core-like collision events. This is further illustrated in Fig. 6, which shows probability distribution functions for accelerations (in arbitrary units). There is a significant tail corresponding to higher accelerations in the case n = 50, compared to the case n = 10. This implies, quite expectedly, that as the interaction steepens, the number of short time collisions with higher accelerations increases and hence the characteristic collision time decreases. In this context, it should be pointed out that the necessity of a careful consideration of the time-scale of two-particle collisions was discussed already by Schofield 49 . In particular, he argued that if there is a sharp hard core, then the duration of a collision τ D tends to zero. If high frequency elastic solid-like excitations exist, they can only appear above a background of width τ − D 1 49 . One more consideration can be based on the observed behaviour of the instantaneous bulk modulus. Figure 3 demonstrates that the actual adiabatic bulk modulus, as inferred from the sound velocities measured in MD simulations, remains finite and smoothly approaches its HS limiting value. In contrast, the high frequency bulk modulus diverges as ∝ ∞ K n when n increases. It is instructive to consider the main steps of the derivation of . The time is normalized in such a way, that all VAF curves intersect zero at the same point, Ω is roughly the Einstein frequency 3 . The inset shows the dependence of the minimum value of VAF on the IPL exponent n, revealing the saturation-like behavior at large n. Figure 5. Particle trajectories within a cubic box with the side length  a 5 , chosen near the center of the simulation cell, demonstrating the particle motion in IPL melts for n = 10 (a) and n = 50 (b). The trajectories are color-coded according to particle energy (transition from red to blue corresponds to the decrease in energy). The time interval during which the trajectories are followed corresponds roughly to one hundred of inverse Einstein frequencies.
the analytical expression for ∞ K 49,50 . A particularly simple route to derive its excess (associated with potential interactions between the particles) part has been discussed previously 9,29 . The starting point is the virial expression for the excess pressure, which should be differentiated with respect to the particle density. In doing so a term containing an explicit derivative ρ ∂ ∂ g r ( )/ appears. The latter clearly depends on the thermodynamic process. In the infinite frequency (instantaneous) limit one assumes that no relaxation is allowed and the density change occurs without any structural rearrangement. Mathematically, this implies that g r ( ) does not change if r is scaled by ρ −1/3 . This allows to express ρ ∂ ∂ g r ( )/ in terms of ∂ ∂ g r r ( )/ and then, integrating by parts, to recover the conventional results 49,50 . It is the assumption of no structural rearrangement [independence of ρ g r ( ) 1/3 of ρ] that causes problems. While well justified for sufficiently soft interactions, it is clearly not applicable to HS like interaction, because an intrinsic length scale -the hard sphere diameter -emerges. Thus, it is legitimate to say that the divergence of the high frequency elastic moduli in the limit of very steep interactions is unphysical. Rather, the conventional expressions must not be applied in this limit.
Of course, all these considerations are merely qualitative. Only direct comparison between the dispersion relations from MD simulations and the QCA model allowed us to quantify the limits of the applicability of elastic approaches to collective modes in fluids as  n 20 for the IPL model. On the other hand, none of these considerations is directly related to the special shape of the IPL potential studied here. Therefore, one should expect that elastic approaches such as QCA also fail for other systems of strongly interacting particles, provided the interaction potential becomes too steep.
The main results of our investigation can be formulated as follows: (i) The QCA description of collective modes fails for sufficiently steep interactions; For the IPL model considered here this failure occurs at ∼ n 20. (ii) For softer potentials QCA can provide a reasonable description of the high-frequency portion of the longitudinal dispersion relation (in particular, in the vicinity of the first maximum); However, the QCA elastic sound velocity overestimates the actual sound velocity. (iii) The instantaneous sound velocity, evaluated using the instantaneous bulk modulus, is an appropriate measure of the actual (adiabatic) sound velocity observed in MD simulations; This instantaneous sound velocity ∞ C ( ) is related to the elastic longitudinal C ( L ) and transverse C ( T ) sound velocities of the QCA approach via 2 . (iv) In the limit of soft interactions (not considered here, but previously extensively studied using Coulomb and Yukawa model systems), the strong inequality  C C L T generally holds, which implies that the longitudinal QCA elastic sound velocity, instantaneous sound velocity, and the thermodynamic sound velocity are all close to each other.
(v) The conventional expressions for high frequency (instantaneous) elastic moduli should not be employed in the hard-sphere interaction limit, where they diverge.
The main finding can also be briefly reformulated as follows. The QCA model was formulated originally to describe collective motion in simple fluids without any explicit limitations regarding the shape of the interaction potential. It turns out, however, that this model works very well for soft interactions (in particularly in the plasma-related context), but fails when the interaction becomes sufficiently steep. We located this failure here using the IPL family of potentials.
In addition to this, two important relevant observation have been reported. It was found that the value of radial distribution function g(r) at the first minimum remains practically constant for all investigated IPL fluids  Figure 6. Probability distribution functions (PDF) of the particle accelerations in IPL melts for two different exponents: n = 10 (blue solid curve) and n = 50 (red dashed curve). The inset shows the ratio of the potential to kinetic energies stored in the system of interacting particles. The ratio drops below unity at  n 30, which is close to the estimated failure of the QCA. Note that for n > 30 the IPL melts are strongly correlated, but not strongly coupled.
with > > n 100 10 near the fluid-solid coexistence. This quantity is shown to exhibit much smaller relative variation than the ratio of the first minimum to the first maximum of g(r) appearing in the Raveché-Mountain-Streett freezing criterion. For the same range of steepness, the ratio of the sound to thermal velocities also remains practically constant, slowly increasing towards the hard-sphere limiting value as n increases. These two observations can potentially be used as approximate structural and dynamical freezing indicators for simple systems with steep repulsive interactions.

QCA dispersion relations.
Within the QCA approach the dispersion relations of the longitudinal and transverse modes are written as 7,8 , where ρ is the particle number density, m is the particle mass, ω is the frequency, k is the wave number, and θ = z rcos is the direction of the propagation of the longitudinal wave. The frequency spectra are thus determined by the first and second derivatives of the interaction potential V r ( ) and the equilibrium radial distribution function of the fluid g(r); the explicit expression for ω L (after intergration over angles) can be found in the literature 39,51 .
As displayed in Eqs (2) and (3), the dispersion relations within QCA are completely determined by the potential interaction between the particles. In related approaches (e.g. based on sum rules or generalized high-frequency elastic moduli) kinetic terms can naturally appear additively to the potential terms ( k v 3 T 2 2 for the longitudinal mode and k v T 2 2 for the transverse mode, where = v T m / T is the thermal velocity) 3,6 . Shofield argued that these kinetic terms should not be included in the dispersion relation 49 . We can avoid further related discussion here by pointing out that the kinetic contributions are normally numerically small at liquid densities 49 . Since in the present paper we consider very dense liquids, just in the vicinity of the solidification, it would be more than sufficient to keep only the potential terms in the dispersion relations (2) and (3).
In the limit of long wavelengths the dispersion relations of the IPL fluid exhibit an acoustic dispersion and the longitudinal and transverse sound velocities can be introduced, Further, for the IPL systems these sound velocities can be easily related to the reduced excess pressure of the fluid, ρ = − p P T / 1 ex , as follows: which is a general (independent of the interaction potential) property of the QCA model in both two-and three dimensions 40 . In addition, the longitudinal and transverse sound velocities are related to the high frequency shear modulus ∞ G and the high frequency bulk modulus ∞ K 6, 50 , This is very similar to the conventional relations for elastic waves in an isotropic medium 52,53 . The "instantaneous" sound velocity can be introduced using the high frequency (instantaneous) bulk modulus 49 This instantaneous sound velocity ∞ C appears to be rather close to the conventional thermodynamic sound velocity for soft repulsive potentials, as has been documented for instance for the case of strongly coupled Yukawa fluids 40 . But this is not the case for steeper interactions as demonstrated in this paper.

Numerical simulations.
Simulations were performed on graphics processing unit (NVIDIA Tesla K80) using the HOOMD-blue software 54,55 . We used = N 55296 particles in a cubic box with periodic boundary conditions. Simulations were performed in the canonical ensemble (NVT) using the Langevin thermostat at a temperature ε = = T 1. A cut-off radius for the potential r max was set such that = × − V r ( ) 5 10 max 8 . The numerical time step was set to be approximately − 10 2 of the inverse Einstein frequency for each exponent n investigated. The system was first equilibrated for fifteen million time steps, and then,the particle positions and velocities were saved for 80 000 time steps.
SCIeNtIFIC REPORTS | 7: 7985 | DOI:10.1038/s41598-017-08429-5 The IPL exponent varied from n = 10 to n = 100. For each chosen exponent n the system state corresponded to the fluid side of the fluid-solid (fcc lattice) coexistence region. This was achieved by choosing the reduced density ρσ 3 according to the data tabulated by Agrawal and Kofke 34 .
For each simulation runs the following diagnostic tools were implemented. We calculated the radial distribution function g(r) and the normalized velocity autocorrelation function (VAF) 2 where t v( ) is the velocity of a particle at time t and 〈⋅⋅⋅〉 denotes an average over all particles. In addition, the longitudinal component of the particle current i N i t k r L 1 ( ) i was calculated. The Fourier transform in time was performed to obtain the current fluctuation spectrum and, hence, the longitudinal dispersion, similarly to our previous approach 23 .