Sound Velocities of Lennard-Jones Systems Near the Liquid-Solid Phase Transition.

Longitudinal and transverse sound velocities of Lennard-Jones systems are calculated at the liquid-solid coexistence using the additivity principle. The results are shown to agree well with the "exact" values obtained from their relations to excess energy and pressure. Some consequences, in particular in the context of the Lindemann's melting rule and Stokes-Einstein relation between the self-diffusion and viscosity coefficients, are discussed. Comparison with available experimental data on the sound velocities of solid argon at melting conditions is provided.

1. Introduction

The high frequency (instantaneous) elastic moduli and the corresponding sound velocities are important characteristics of condensed fluid and solid phases, which affect and regulate wave propagation [1], the instantaneous Poisson’s ratio [2], the coefficient in the Stokes–Einstein relation [3,4], the Lindemann melting rule [5,6], as well as some other simple melting rules [7], the relaxation time in the shoving model [8,9], just to give few examples.

Several useful empirical observations exist. For example, it is well known 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 [10,11,12]. A close value (≃9.5) has been reported from experiments with solid argon at the melting temperature [13]. Values close to ≃9 have also been reported for solid hydrogen and deuterium along the melting curve [14]. Rosenfeld pointed out that this “quasi-universal” property is also shared by the hard-sphere (HS) model [12]. More recently, it has been demonstrated that this property is also exhibited by the purely repulsive soft inverse-power-law (IPL) model in a wide range of IPL exponents [15]. An important related question is how the presence of long-range attraction (in addition to short-range repulsion) can change this picture.

It has been recently demonstrated that the Lindemann’s criterion of melting can be re-formulated for two-dimensional (2D) classical solids using statistical mechanics arguments [6]. With this formulation the expressions for the melting temperature are equivalent in three dimensions (3D) and 2D. An important consequence of this formulation is that the ratio of the transverse sound velocity to the thermal velocity is predicted to have a quasi-universal value along the melting curve (which can be different in 3D and 2D). This condition has been verified on soft repulsive interactions (in particular, using the Yukawa or exponentially screened Coulomb potential) in both 2D and 3D, where it works reasonably well [6,16]. A natural question arises whether this condition remains valid and useful also for potentials with a long-ranged attraction.

Motivated by these and related questions, we have investigated in detail how the longitudinal and transverse sound velocities behave at the liquid–solid phase transition of 3D Lennard-Jones (LJ) systems. The high-frequency (instantaneous) velocities have been calculated using standard expressions along the liquidus and solidus of LJ system. To facilitate the calculations, we elaborate on the principle of additivity of the melting and freezing curves put forward by Rosenfeld [17,18]. We are then able to quantify the behavior of sound velocities at the liquid–solid phase transition and make comparison with purely repulsive soft sphere systems. Towards the end of the paper, we will also present a comparison with available experimental data on sound velocities of solid argon at melting conditions.

2. Approach

2.1. Formulation

For an arbitrary spherically symmetric pairwise repulsive potential , the longitudinal and transverse velocities can be expressed (in 3D) as follows [19,20,21] and Here, and are the longitudinal and transverse elastic sound velocities, m is the particle mass, the sums run over all neighbours of a given particle, and primes denote derivatives of the interaction potential with respect to the distance r.

The representation of sound velocities used above is based on the relations between the sound velocities and elastic moduli, and , where is the particle number density. In solids M and G are the conventional longitudinal and shear elastic moduli and the summation is over ideal lattice sites at zero temperature. For simplicity we assume that sound velocities are isotropic. This essentially corresponds to some effective sound velocities, averaged over the directions of propagation. The bulk elastic modulus is . Elastic modes softening at finite temperatures is not accounted for in this formulation. In fluids, M and G correspond to the so-called infinite frequency (instantaneous) longitudinal and shear moduli [22] (often denoted as and ) and the summation should be performed using the actual liquid structure. This summation is usually replaced by integration involving the radial distribution function : [19,22,23]. For our present purposes summation just keeps the notation somewhat more compact. In the liquid state additional kinetic terms should appear in Equations (1) and (2), but these are numerically small near the liquid–solid coexistence and are therefore omitted for simplicity.

In the plasma-related context, expressions (1) and (2) with integration instead of summations are familiar as the quasi-localized charge approximation [24,25,26,27,28].

Important thermodynamic properties of a system of interacting particles are the internal energy and pressure. For pairwise interactions the excess (over the ideal gas) contributions to the energy, , and pressure, , can be expressed via summations similar to those used above and Reduced units have been used, , , where the temperature T is measured in energy units. Trivial manipulation with Equations (1) and (2) allows us to obtain the relation between the sound velocities and the excess pressure where denotes the thermal velocity. Equations (5) is known as Cauchy relation, it is valid independently of whether kinetic terms are included or not (simply because they cancel out if retained) [22]. For the Lennard-Jones potential we can also express the sound velocities in terms of and [22]. The corresponding expressions are and The Cauchy relation is obviously satisfied.

2.2. Inverse-Power-Law Model

As an important reference case let us first consider the soft-sphere model near the fluid–solid phase transition [15]. The IPL potential is defined as where and are the energy and length scales, and n is the IPL exponent. In this case the sound velocities are directly related to the excess pressure of the system. The corresponding relations are directly obtained from Equations (1), (2), and (4) and read Note that for the transverse sound velocity does not become negative, because in this regime the neutralizing background should be included, which makes the excess pressure negative [29]. We will not consider this regime.

The sound velocities have been evaluated along the fluid–solid coexistence boundaries using the coexistence properties tabulated in Ref. [30]. The inverse IPL exponent, referred to as the softness parameter varies in the range . The solid near melting is in the bcc phase for sufficiently soft interactions with and forms an fcc solid otherwise. Since the thermodynamic properties of the bcc crystal at melting are very nearly the same as those of the fcc crystal [30], it suffices to perform calculations for the fcc-fluid phase transition. The results are shown in Figure 1. The following main trends are observed: (i) Very weak dependence of both longitudinal and transverse sound velocities on the softness parameter in the considered range of softness; (ii) Numerical values are comparable to those of the hard-sphere fluids at the freezing packing fraction [2]; (iii) The difference between the sound velocities in the fluid and solid phases is very tiny and can normally be neglected; (iv) the longitudinal sound velocity exhibits a minimum at , while the transverse sound velocity decreases continuously as s increases; (v) the ratio of sound velocities, , decreases monotonously from ≃0.53 to as s increases in the range considered.

Figure 1

Reduced longitudinal and transverse sound velocities at the fluid–solid coexistence of the soft sphere model versus the softness parameter s. Sound velocities are expressed in units of the thermal velocity . Upper curves are for the longitudinal mode, lower curves are for the transverse mode. Solid (dashed) curves correspond to the solid (fluid) boundary of the fluid–solid coexistence.

Note that we cannot trace the transition to the HS limit using Equations (9) and (10). They predict divergence of sound velocities as (or ), which contradicts finite values in the HS limit [31,32]. The origin behind the unphysical divergence of the conventional expressions for the instantaneous elastic moduli when approaching the HS limit has been identified and discussed [2,33]. The conventional expressions should not be applied for .

3. Results

3.1. Additivity of Melting Curves

The Lennard-Jones potential is where and are again the energy and length scales (or LJ units), respectively. The density, temperature, pressure, and energy expressed in LJ units are , , , and . Relation to the reduced excess units is straightforward: and .

It has been long known, from the results of Monte Carlo (MC) simulations, that details of the interaction potential have relatively little effect on the structure of fluids near the melting temperature, in particular when extreme cases of HS and Coulomb interactions are excluded from consideration [34]. The same concerns some simple melting characteristics such as the Lindemann ratio, reduced free volume, amplitude of the first maximum of the structure factor. These observations allowed Rosenfeld to formulate his principle of additivity of the melting curves [17,18]. This principle states that if the stable pairwise interaction potential represents a linear combination of stable repulsive potentials, then the temperature and pressure along the liquidus and solidus can also be expressed as a linear combinations of temperatures and pressures corresponding to individual repulsive potentials at freezing and melting [17]. Further support to the Rosenfeld’s point of view is provided by the concept of isomorphs [35], which are the curves along which structure and dynamics in properly reduced units are invariant to a good approximation [36,37]. Many simple systems, including the LJ case exhibit isomorphism. Melting and freezing curves appear as approximate (although not exact) isomorphs [38,39,40], and this can simplify considerably calculation of system properties at melting and freezing.

The fact that the melting and freezing curves are approximate isomorphs indicates that structures are nearly invariants when properly scaled units are used. For instance, if the distance is measured in units of characteristic interparticle separation a, , then the sums are independent of density to a good approximation. This implies For an ideal zero-temperature crystal, would correspond to the lattice sum of the IPL-n potential. Using this property the expressions for the excess reduced energy, pressure, and sound velocities can be written in the compact form as where is the required quantity and and are the corresponding numerical coefficients. These coefficients are summarized in Table 1. Equation (13) can be referred to as additivity of excess pressure, energy, and sound velocities along melting and freezing curves.

Table 1

Coefficients and from Equation (13).

X=	u_ex	p_ex	cl2/vT2	ct2/vT2
C_12	2	8	2965	725
C_6	2	4	765	125

A scaling, analogous to that of Equation (13), was discussed in Refs. [40,41]. From numerical data it was observed that vary noticeably with density close to the triple point for both phases. Away from the triple point, they become almost constant. We have chosen to take the sums as constant and to evaluate them from the known coexistence data for IPL and potentials [30,42]. The results are summarized in Table 2. This procedure can be straightforwardly generalized to the case of a general LJ () potential.

Table 2

Solid-liquid coexistence data [42] and numerical values of the sums for and IPL potentials.

n	P_∗	ρ∗sol	ρ∗liq	Σ^sol	Σ^liq
n=12	23.74	1.211	1.167	4.325	5.214
n=6	105.0	2.358	2.330	7.829	8.106

Figure 2 and Figure 3 demonstrate the accuracy of the additivity principle (13) when applied to calculate the reduced pressure and excess energy of the LJ systems at the fluid–solid coexistence. MC data are those tabulated in Ref. [43]. Note that is constant at the coexistence, while is not, because density is used in the normalization. We have used solid densities and to plot the curve in Figure 2. Very similar results would be obtained with liquid densities and . Overall, the agreement between MC results and theory is rather convincing.

Figure 2

Reduced coexistence pressure, , of the Lennard-Jones system versus the reduced temperature . Symbols correspond to MC results [43]. The curve is calculated using the additivity principle (13).

Figure 3

Reduced excess energy per particle in LJ units, , of the Lennard-Jones system versus the reduced temperature . Symbols correspond to the MC results [43]. The curves are calculated using the additivity principle. Solid symbols and curve correspond to the solid side of the liquid–solid coexistence (solidus); open symbols and the dashed curve correspond to its liquid side (liquidus).

3.2. Sound Velocities of the LJ System

The sound velocities of the LJ system at the liquid–solid coexistence are plotted in Figure 4. Curves are calculated using the additivity principle (13). Symbols correspond to relations (6) and (7) using the MC data for excess energy and pressure [43]. The agreement between these two methods is good everywhere, except in the vicinity of the triple point. This is not so surprising, because in this region the two large terms associated with IPL-12 repulsion and IPL-6 attraction are almost comparable in magnitude, so that even a small relative inaccuracy in each term can result in a much greater relative inaccuracy of their difference.

Figure 4

Reduced longitudinal and transverse sound velocities of the LJ model versus the reduced temperature . Upper symbols and curves are for the longitudinal mode, lower symbols and curves are for the transverse mode. Symbols are the results of calculation using relations (6) and (7) using MC data from Ref. [43]. Solid and open symbols correspond to the boundaries of the solid and liquid phases, respectively. Solid (dashed) curves correspond to the solid (liquid) coexistence boundary and are plotted using the additivity principle (13).

The reduced sound velocities are to a good approximation constant there with and . The ratio of sound velocities is approximately constant with . The difference between the sound velocities at the solid and liquid coexistence boundaries is vanishingly small away from the triple point, similarly to the IPL case.

Now we can elaborate on some consequences of the obtained results. The presence of long range attraction seems not to affect (or, more precisely, affects rather weakly) the sound velocities in the vicinity of the fluid–solid phase transition. For comparison, the sound velocities of IPL-12 system at melting is and . For the IPL-6 system we have and .

The Lindemann’s criterion of melting can be expressed as [5,6] where is the Debye frequency, a is the characteristic interparticle separation (the Wigner–Seitz radius is used here), and C is expected to be a quasi-universal constant. In this form the Lindemann expression for the melting temperature applies to both 3D and 2D solids (the constants C can be different for the 2D and 3D cases) [6]. The Debye frequency in 3D can be expressed via the longitudinal and transverse sound velocity as For soft repulsive interactions the strong inequality usually holds and this can be used to further simplify the melting condition [6]. In the considered case , and we use the actual values and to get . The resulting value is somewhat below , previously reported for a one-component Coulomb plasma (which corresponds to an extremely soft interaction potential with ) [6]. On the other hand, constancy of either or , combined with the scaling (13) immediately yield freezing and melting equations of the form which is a very robust result reproduced in a number of various theories and approximations [17,18,38,39,40,41,44,45].

Zwanzig’s result for the relation between the self-diffusion and viscosity coefficients of liquids [3] can be expressed in the form of the Stokes–Einstein (SE) relation with a coefficient that depends on the ratio of the transverse to longitudinal sound velocities [4], where D is the coefficient of self-diffusion, , is the viscosity, and is the SE coefficient. For the SE coefficient becomes , which is in remarkable agreement with recent extensive MD simulation results demonstrating that not too far from the fluid–solid coexistence [46].

3.3. Comparison with Experiment

Measurements of the longitudinal and transverse ultrasonic wave velocities in compressed, solidified argon for pressures up to 6 kbar (600 MPa) corresponding to melting temperatures in the range 123–206 K have been reported in Ref. [13]. Experimental results, in the form of the temperature dependence of the reduced sound velocities and , are plotted in Figure 5. We observe that the reduced sound velocities are constant to a very good accuracy in the temperature range investigated. The ratio is close to 0.5 in agreement with theoretical expectations. The numerical values and are however by about lower than the theory predicts (we have used K and cm as LJ parameters for argon [47]). In this context, we should remind that the theory of sound velocities considered here is idealized in a sense that it takes into account the temperature effect on the system structural properties, but does not take into account the effect of thermal fluctuations. The latter are known to somewhat reduce the values of elastic constants (elastic moduli softening) [48] and become important as the melting transition is approached. For many metallic solids the reduction in the shear modulus (compared to zero-temperature conditions) amounts to 20%–30% [49,50]. Similarly, a reduction was observed in the one-component plasma model [51]. Thus, the trends observed for the LJ system are not unique. Additionally, we should not ignore the fact that the pairwise Lennard-Jones potential function itself may not be the best representation of a real interaction potential in solid argon.

Figure 5

Reduced longitudinal and transverse sound velocities in compressed solidified argon at the melting temperature in the range 123–206 K. Symbols correspond to the experimental results tabulated in Ref. [13].

4. Discussion and Conclusions

The longitudinal and transverse sound velocities of the Lenard-Jones system have been evaluated at the liquid–solid coexistence. Two methods have been employed, one uses the additivity principle and the other uses relations between sound velocities and excess energy and pressure. The first method is simple but approximate, while the second is exact, but requires the knowledge of the excess energy and pressure. The agreement between the two methods is rather good, the deviations are only observable near the triple point. This is not surprising, because in this region the repulsive and attractive contributions are comparable in magnitude, so that even a small relative inaccuracy in each of these terms can result in a much greater relative inaccuracy of their difference. Nevertheless, even near the triple point the results based on the additivity principle demonstrate acceptable accuracy.

The calculated ratios of sound velocities to the thermal velocity are practically constant along the melting and freezing curves. The difference between the sound velocities in the solid and liquid phases is insignificant. The numerical values and are comparable with those of repulsive soft sphere (IPL) and HS models. Thus, long-range attraction seems not to affect sound velocities considerably.

The latter conclusion should not be understood too literally. The structure of the LJ system in the vicinity of the fluid–solid phase transition is merely determined by the short-range repulsive branch of the interaction potential at distances about the mean inter-particle separation [52]. In this range the potential can be approximated by the IPL shape with an effective IPL exponent (generally, ) [52]. The value of is considerably affected by the attractive term of the LJ potential. However, since and are only weakly dependent on (see, e.g., Figure 1), the effect of attractive force is small in the considered context. The opposite situation is not impossible in other circumstances. For example, the presence of a long-range attractive (dipole-like) interaction can considerably suppress the sound velocity in complex plasma fluids (fluids composed of macroscopic charged particles immersed in a plasma environment) [53,54].

The obtained results have been analysed in the context of the Lindemann’s melting rule and Stokes–Einstein relation. The constancy of and is consistent with the functional form for the dependence of melting and freezing temperatures on density, which emerge in various theories and approximations. The ratio of the transverse and longitudinal sound velocities allows to evaluate the numerical coefficient in the Stokes–Einstein relation, and the result agrees remarkably well with that from recent MD simulations.

A comparison with available experimental data on the sound velocities of solid argon at melting conditions has been provided. It is demonstrated that the ratio in experiment is close to 0.5 in agreement with theoretical expectations. The ratios and are, however, about lower than the theory predicts. This can be a consequence of elastic moduli softening on approaching the melting temperature, which is not taken into account in the theory. Additionally, the LJ potential may not be the best model of real interactions in solid argon.

As a final remark, we should point out that the results presented here can be easily generalized to the case of the () Lennard-Jones potential.