Thermodynamics of freezing and melting

Abstract

Although the freezing of liquids and melting of crystals are fundamental for many areas of the sciences, even simple properties like the temperature–pressure relation along the melting line cannot be predicted today. Here we present a theory in which properties of the coexisting crystal and liquid phases at a single thermodynamic state point provide the basis for calculating the pressure, density and entropy of fusion as functions of temperature along the melting line, as well as the variation along this line of the reduced crystalline vibrational mean-square displacement (the Lindemann ratio), and the liquid’s diffusion constant and viscosity. The framework developed, which applies for the sizable class of systems characterized by hidden scale invariance, is validated by computer simulations of the standard 12-6 Lennard-Jones system.

Introduction

Melting is the prototypical first-order phase transition^1,2,3. Its qualitative description has been textbook knowledge for a century, but it has proven difficult to give quantitatively accurate predictions. This is the case not only for the kinetics of freezing and melting, which are exciting and highly active areas of research^4,5,6,7,8; there is not even a theory for calculating, for example, the entropy of fusion as a function of temperature along the melting line in the thermodynamic phase diagram.

The everyday observation that matter sticks together but is at the same time almost impossible to compress⁹ is modelled, for example, in the system proposed by Lennard-Jones (LJ) in 1924 (ref. 10). Here, particles interact via a pair potential that as a function of distance r is a difference of two inverse power-law terms: υ_LJ(r)=4ɛ((r/σ)⁻¹²−(r/σ)⁻⁶). The first term reflects the fact that the repulsive ‘Pauli’ forces are harsh and short-ranged, the negative term models the softer, longer ranged attractive van der Waals forces. The 1970s led to the development of highly successful thermodynamic perturbation and integral-equation theories for simple liquids^{11,12,13,14,15,16}. Their main ingredient is the assumption that the structure of a dense, monatomic fluid closely resembles that of a collection of hard spheres^14,16,17,18. Confirming this, the structure of melts of, for example, metallic elements near freezing is close to that of the hard-sphere system^15,16,18,19. The term ‘structure’ generally refers to the entire collection of spatial equal-time density correlation functions, but our focus below is on the pair correlation function (in the form of its Fourier transform, the structure factor) as the most important structural characteristic.

Since the hard-sphere system has only a single nontrivial thermodynamic state parameter, the packing fraction, the phase diagram is basically one-dimensional, which implies that the system has a unique freezing/melting transition. On the basis of this, for simple systems one expects invariance along the freezing and melting lines of structure and dynamics in proper units, as well as of thermodynamic variables like the relative density change upon melting and the melting entropy²⁰. Empirical freezing and melting rules, which follow from the hard-sphere melting picture and are fairly well obeyed for most simple systems, include the fact that the ratio between the crystalline root-mean-square atomic displacement and the nearest-neighbor distance—known as the Lindemann ratio—is constant and about 0.1 along the melting line; this is the famous Lindemann melting criterion from 1910 (refs 20, 21, 22, 23, 24, 25). In the hard-sphere model the Lindemann ratio is universal at melting because, as mentioned, there is just a single melting point. Thus, for systems well described by the hard-sphere model the Lindemann ratio is predicted to be invariant along the melting line. Other empirical rules, which are predicted by the hard-sphere picture and reasonably well obeyed by many systems, include the facts that in properly reduced units the liquid’s self-diffusion constant and viscosity are invariant along the freezing line^26,27, the Hansen–Verlet rule^17,28 that the amplitude of the first peak of the liquid static structure factor is about 2.85 at freezing, or Richard’s melting rule³ that the entropy of fusion ΔS_fus is about 1.1k_B (which in a more modern and accurate version is the fact that the constant-volume entropy difference across the density–temperature coexistence region is close to 0.8k_B (refs 23, 29)).

The below study shows how the thermodynamics of freezing and melting for a large class of systems may be predicted to a good approximation from computer simulations carried out at a single coexistence state point. In particular, the theory developed quantifies the deviations from the above mentioned hard-sphere predicted melting-line invariants^{16,22,30,31,32}. The theory is validated by computer simulations of the standard 12-6 LJ system.

Results

General theory

It is well-known that adding a mean-field attractive term to the hard-sphere model broadens the coexistence region, which on the other hand, narrows if the repulsive part is softened^{13,16,33,34,35,36}. Such terms are therefore expected to modify the hard-sphere predicted invariances along the freezing and melting lines. As an illustration, Fig. 1a shows that in reduced units there is approximate identity of structure along the LJ freezing line, but the structure is not entirely invariant as seen in the inset where the dashed line marks the predicted maximum based on simulations at T₀=2.0ɛ/k_B, if the structure were invariant.

Figure 1: Structure of the LJ liquid.

(a) Liquid structure factor along the freezing line³⁷ showing results from T=0.7ɛ/k_B, which is close to the triple point, to T=3.4ɛ/k_B. The hard-sphere model predicts that the height of the first peak is invariant along the freezing line as indicated by the blue dashed line in the inset. Small, but systematic deviations are observed. (b) Liquid structure factor along the isomorph crossing the freezing line at temperature T₀=2.0ɛ/k_B (henceforth used as the liquid reference isomorph), demonstrating structural invariance to a much higher degree. This is the basis for the theory proposed in the present paper.

In order to develop a quantitative theory of freezing and melting, we take as starting point the ‘hidden scale invariance’ property of systems³⁸ characterized by a potential-energy function U(R), where R=(r₁, r₂,…,r_N) is the collective coordinate of the system’s N particles, which to a good approximation obeys the scaling condition.³⁹

Here, λ is a scaling factor and it is understood that the sample container undergoes the same scaling as the configuration; thus λ>1 corresponds to a density decrease and λ<1 to a density increase. This form of scale invariance is exact only for systems with Euler-homogeneous interactions (plus a constant)¹³. It is a good approximation, however, for the condensed phases of many systems in which this property is not obvious from inspection of the analytical expression for U(R), thus the term ‘hidden scale invariance’^39,40,41,42. Equation (1), which is formally equivalent to the conformal-invariance condition U(R_a)=U(R_b)⇔U(λ R_a)=U(λ R_b), implies invariance of structure and dynamics along the configurational adiabats in the phase diagram³⁹. These lines are referred to as isomorphs⁴². It was very recently shown by Maimbourg and Kurchan⁴³ that in high dimensions all pair-potential systems obey hidden scale invariance in their condensed phase. Experimentally, hidden scale invariance has been demonstrated directly and indirectly for molecular van der Waals bonded liquids and polymers^44,45,46. Further evidence for the existence of isomorphs comes from computer simulations of single-component systems^40,42 as well as, for example, of glass-forming mixtures⁴⁷, nanoflows⁴⁸, molecular models³⁸ and molecular dynamics (MD) simulations of the dynamics of most metallic elements based on quantum-mechanical, density-functional-theory potentials⁴⁹. Isomorphs have also been demonstrated in simulations of out-of-equilibrium situations like zero-temperature shear flows of glasses or nonlinear steady-state liquid flows (see, for example, ref. 38 and its references). It is important to emphasize, however, that not all condensed matter exhibits hidden scale invariance; for instance, water is a notable exception⁴¹. The general picture is that most metals and organic van der Waals bonded systems obey equation (1) to a good approximation in the condensed-phase part of their thermodynamic phase diagram, whereas systems with strong directional bonding generally do not³⁸. The former systems are simpler than the latter because their phase diagrams are effectively one-dimensional in regard to structure and dynamics, reminiscent of the hard-sphere system. Systems with hidden scale invariance are sometimes referred to as Roskilde (R) simple^{35,50,51,52,53,54,55,56,57,58,59,60,61,62} to distinguish them from simple systems traditionally defined as pair-potential systems¹⁶. The theory presented below makes use of R simple systems’ almost one-dimensional phase diagrams³⁸ and gives corrections to the hard-sphere picture of melting and freezing calculated by the first-order Taylor expansions. Figure 2 illustrates the idea.

Figure 2: Illustration of the main idea of the theory.

The freezing and melting lines are both close to isomorphs along which basically everything is known because the reduced-unit structure and dynamics are invariant to a very good approximation. Properties along the freezing and melting lines are estimated via first-order Taylor expansions by moving from an isomorph to the freezing or melting line; the two reference isomorphs (a liquid and a solid one) are determined from computer simulations at T₀=2.0ɛ/k_B. Details are given in the ‘Methods’ section.

Along an isomorph the structure is invariant in the reduced-unit system defined⁴² by the length unit ρ^−1/3 (ρ≡N/V is the number density and V is the system volume), the energy unit k_BT (T is the temperature) and the time unit (m is the particle mass). Figure 1b shows the LJ liquid’s static structure factor S(q) along an isomorph close to the freezing line (used below as the liquid-state reference isomorph) plotted for a range of temperatures. A comparison with Fig. 1a confirms the recent finding of Heyes and Branka³² that the freezing line is not an exact isomorph, although it is close to one.

The melting pressure as a function of temperature, p_m(T), can be predicted from information obtained at a single coexistence reference state point. The details about how this works are given in the ‘Methods’ section. The argument may be summarized as follows. Recalling that the entropy as a function of density and temperature is a sum of an ideal-gas term and an ‘excess’ term S_ex (ref. 16), isomorphs are the phase-diagram lines of constant excess entropy for any system obeying equation (1)^39,42. A computer simulation at the liquid/solid reference state point generates a series of configurations R₁⁰,…, R_n⁰. Scaling each of these uniformly to density ρ one obtains configurations representative for the state point with density ρ and temperature T on the isomorph through the reference state point³⁹ in which T is identified from the configurational temperature expression⁶³ k_BT=〈(∇U)²〉/〈∇²U〉. The average potential energy U and virial W at the state point (ρ, T) are likewise found by averaging over the scaled configurations. The key assumption here is that the canonical probabilities of the scaled configurations are identical to those of the original configurations, which follows from equation (1)³⁹ (thus no new MD simulations are required). As shown in the ‘Methods’ section, in conjunction with the excess isochoric specific heat C_V^ex calculated from the potential-energy fluctuations of the scaled configurations (C_V^ex=〈(ΔU)²〉/k_BT²) and the so-called density-scaling exponent also calculated from the fluctuations (γ=〈ΔUΔW〉/〈(ΔU)²〉), one has enough information to determine the thermodynamics of freezing and melting, as well as the variation along the melting line of isomorph-invariant properties like the Lindemann ratio and the reduced-unit viscosity.

The LJ system

For LJ type systems, the general procedure described above may be implemented analytically by making use of the fact that because the structure is isomorph invariant, it is possible to calculate the variation of the average potential energy and other relevant quantities analytically along an isomorph. This is done as follows. In reduced coordinates the pair correlation function is isomorph invariant . Consequently, for pairs of LJ particles at distance r the thermal average 〈r⁻ⁿ〉 scales with density as ρ^n/3 along an isomorph. Thus 〈r⁻ⁿ〉 ∝ ρ^n/3 with a proportionality constant that only depends on S_ex, implying that the average potential energy U is of the form⁶⁴ in which is the density relative to the reference state-point density and A₆(S_ex)<0 derives from the attractive term of the LJ potential. Since T=(∂U/∂S_ex)_ρ, one has . It follows that if the five quantities S_ex, A₁₂(S_ex), A₆(S_ex), and are known, the excess Helmholtz free energy, U–TS_ex, is known along the reference isomorph. The required quantities are easily determined from reference state-point simulations (see the ‘Methods’ section)—for instance the reference state-point’s potential energy and virial give two linear equations for determining A₁₂(S_ex) and A₆(S_ex). Once the excess Helmholtz free energy is known along the reference isomorph, the Gibbs free energy is found by adding the ideal-gas Helmholtz free energy and the pV term (pV=Nk_BT+W in which the virial is given⁴² by ).

Comparing theory to simulation results for the LJ system

Following the above procedure, we generated two reference isomorphs for the LJ system starting from the coexistence state point with temperature T₀=2.0ɛ/k_B, a liquid-phase isomorph and a crystal-phase isomorph. Gibbs free energy of the liquid phase at coexistence, G_l(T), is found by utilizing the fact that the freezing line is close to an isomorph. Since (∂G/∂p)_T =V, a good approximation to G_l at coexistence is

Here, p_m(T) is the coexistence pressure to be determined; G_l^I(T) is the Gibbs-free energy, V_l^I(T) the volume and p_l^I(T) the pressure along the liquid-state reference isomorph. These quantities are all known functions of the (relative) density on the isomorph henceforth denoted by , which for temperature T is found by solving .

An analogous expression applies for the crystal’s Gibbs free energy, of course, again involving only parameters determined from reference state-point simulations. The coexistence pressure is determined by equating the liquid and solid phases’ Gibbs free energies. As shown in the ‘Methods’ section (equation (21)), this results in p_m(T)(V_l^I(T)−V_s^I(T))=C₁(T)+C₂(T)−C₃(T) in which C₁(T) is the difference between U_s^I(T)−(T/T₀)U_s^I(T₀) and the analogous term for the liquid reference isomorph (here U_s^I(T) is the crystal’s potential energy along the reference isomorph), C₂(T) is the difference between and the analogous liquid term and C₃(T) is the difference between (T/T₀)W_s^I(T₀) and the analogous liquid term.

Figure 3a,b compare the theoretically predicted p_m(T) to the coexistence pressure computed numerically by means of the interface-pinning method³⁷. The density of the crystalline and liquid phases may also be computed by means of a first-order Taylor expansion working from the reference isomorph (see the ‘Methods’ section). Figure 3c compares the predicted (ρ,T) phase diagram based on equation (26) to that obtained by the interface-pinning MD simulations. Finally, Fig. 3d shows the predicted and simulated fusion entropy ΔS_fus and enthalpy ΔH_fus, the latter quantity being of course measured in experiments as the latent heat. In all cases there is good agreement between theoretical prediction and simulations.

Figure 3: Theoretical predictions (full red curves) and results of MD simulations (black dots) for the LJ system.

The theoretical predictions are based on simulations at the coexistence reference state point indicated by an arrow in each figure (T₀=2.0ɛ/k_B), the MD simulations employed the interface-pinning method³⁷, see the ‘Methods’ section. No fitting was done in these figures—the only input to the theory is properties of the coexisting liquid and crystal at the reference temperature. (a) Temperature–pressure phase diagram. The green dashed line marks the expectation based on a linear extrapolation of the Clausius–Clapeyron relation dp_m/dT=ΔS_fus/ΔV_fus from the reference state point, that is, assuming that the entropy of fusion and the volume change are both constant. (b) The same data plotted with a linear pressure axis. (c) The freezing and melting lines in the density–temperature phase diagram; the coloured area is the coexistence region. (d) Fusion entropy (main panel) and enthalpy (inset).

Having in mind the fact that the pressure at the triple point is very low for the LJ system, we estimate the triple point temperature to T_tp=0.688(2)ɛ/k_B from the theory by assuming zero pressure. This is within the statistical uncertainty of the triple point temperature computed with the interface-pinning method. For comparison, a linear extrapolation of the Clausius–Clapeyron relation from the reference temperature (the green dashed lines in Fig. 3a,b) predicts a triple point temperature of 0.909(2)ɛ/k_B.

Since the melting line is not an isomorph, the Lindemann ratio is not invariant along it. The theory estimates the deviation from a constant Lindemann ratio by a first-order Taylor expansion from the reference isomorph (see Fig. 2 and the ‘Methods’ section). Figure 4a demonstrates good, though not perfect agreement between theory and numerical computations of the Lindemann ratio. The liquids’ self-diffusion constant plays an important role for the crystal growth rate as expressed, for example, in the Wilson–Frenkel law^65,66. This motivated us to use the theory also for calculating the liquid’s diffusion constant variation along the freezing line (Fig. 4b). Another important component for crystal growth is the thermodynamic driving force on the crystal–liquid interface, which is the Gibbs free energy difference between the two phases, ΔG≅(T_m−T)ΔS_fus (ΔS_fus is shown on Fig. 3d). Finally, Fig. 4c shows the viscosity along the freezing line. In all cases the blue dashed lines mark the prediction if the dynamics were invariant in reduced units, that is, if the freezing/melting lines were isomorphs.

Figure 4: Predictions (red curves) versus results of computer simulations (black dots) for three properties along the freezing/melting lines of the LJ system.

The blue dashed lines show the predictions if perfect invariance of structure and dynamics in reduced units applies along the freezing/melting lines, the arrows indicate the reference state point upon which the predictions are based. (a) Lindemann ratio along the melting line. (b) Self-diffusion constant along the freezing line. (c) Viscosity along the freezing line.

Discussion

The theory presented above predicts the thermodynamics of freezing and melting from a single coexistence state point. The theory also enables one to calculate the deviations from the invariance of several quantities along the melting line predicted by the hard-sphere melting picture^{22,16,30,31,32}. The theory is analytic for LJ type systems, that is, systems involving a pair potential that is a difference of two inverse power laws, but the framework developed applies to any system with hidden scale invariance, including molecular systems. The theory works well for the LJ system, with the largest deviations found close to the triple point where the structure is less invariant along the reference isomorph (Fig. 1b).

Having established a firm foundation for the thermodynamics of freezing and melting for R simple systems, it is our hope that it will soon be possible to address the exciting questions of how nucleation and growth proceed, processes that are not well understood even for simple systems beyond the hard-sphere system⁶⁷. It seems likely that variations of the nucleation and growth mechanisms along the melting line can be analyzed in the same way as above, that is, by utilizing the fact that the freezing and melting lines are close to isomorphs along which the dynamics is invariant to a quite good approximation.

It is not clear to which degree this approach to melting can be generalized to quantum systems for which an outstanding question is the possible existence of a zero-temperature quantum fluid of metallic hydrogen. The quantum nature of the proton modifies classical melting, for example by increasing the value of the Lindemann ratio⁶⁸. It would be interesting to investigate whether melting of quantum crystals may be understood in the above framework, but this awaits the development of an isomorph theory for quantum systems. In ongoing work we are addressing another open question, namely whether the above can be generalized to deal with more realistic systems, for instance metals for which density-functional-theory computer simulations nowadays give realistic representations of the physics and have demonstrated hidden scale invariance for most metals⁴⁹.

Methods

Computer simulations

We studied a LJ system of N=5,000 particles with pair interactions truncated and shifted at 6σ. Coexistence pressures, p_m, are computed with the interface-pinning method³⁷ in which coexistence state points are determined by computing the thermodynamic driving force on a solid-liquid interface. Table 1 lists the energy U₀ and virial W₀ at the reference temperature T₀=2ɛ/k_B for both the liquid and crystal states at coexistence. The A₁₂ and A₆ coefficients (for the liquid and the crystal separately) are computed from reference state-point data using equation (8) below. The derivatives of the A coefficients with respect to excess entropy, and , are computed from reference state-point data using equation (11) with the γ₀'s listed in Table 1. Melting pressures (Fig. 3a,b) are computed from reference state-point data using equation (21) in which the potential energies along the two reference isomorphs are expressed in equation (6). The densities along the liquid and crystal reference isomorphs are found as functions of temperature by inversion of equation (9) (upper equation). The second derivatives of the A coefficients, and , are given by equation (15) where the reference state point excess heat capacity and are listed in Table 1. The freezing and melting densities (Fig. 3c) are computed from the pressures by combining equations (22) and (25). The entropy of fusion ΔS_fus (Fig. 3d) is computed by combining equations (27, 28, 29, 30). The value of the Lindemann ratio L of the crystal at the reference temperature, L₀ and its temperature derivative along an isochore, (∂L/∂T)_ρ, are listed in Table 1. By letting X=L in equations (32) and (38), we arrive at the prediction shown in Fig. 4a. Similarly, the predictions of the self-diffusion constant D (Fig. 4b) and viscosity η (Fig. 4c) are found by letting and , respectively. D is determined from the long-time limit of the mean-square displacement; η is computed using the SLLOD algorithm as detailed in ref. 27 except that in the present work we increased the number of particles to 4,096 and used the above-mentioned larger cutoff.

Table 1 Quantities characterizing the two reference state points in coexistence.

We proceed to describe the theory in detail. The reference state point is selected at coexistence, that is, with known temperature T₀ and pressure p₀. There are two different reference densities, a solid and a liquid one, below denoted, respectively, by ρ_s,0 and by ρ_l,0. In the density–temperature phase diagram there are two reference isomorphs. The arguments developed in the next two sections refer to either one of these.

Isomorph characteristics of arbitrary R simple systems

As mentioned, the temperature–pressure reference state point defines two reference density–temperature state points, a liquid and a solid one. Let us focus on one of these with density ρ₀ and temperature T₀ (thus dropping in this and the next subsection subscripts s and l). From an NVT MD equilibrium simulation (for example, with a Nosé–Hoover thermostat) n configurations R₁⁰,R₂⁰,…,R_n⁰ are sampled. In order to map out the reference isomorph parametrized by density, one first identifies the temperature T such that (ρ,T) is on the isomorph through the reference state point (ρ₀,T₀). This is done as follows. If the configurations scaled uniformly to density ρ are denoted by R₁,R₂,…,R_n in which R_i=(ρ₀/ρ)^1/3R_i⁰, the temperature T is determined from the standard configurational temperature expression (in which the averages are over the n sampled configurations)

This determines the function where we define the relative density along the isomorph by with superscript I indicating ‘isomorph’ (thus T(1)=T₀). By averaging the potential energy U(R) and the virial W(R)≡(−1/3)R˙∇U(R) over the scaled configurations one identifies the functions and . is found from the scaled configurations’ potential energy via C_V^ex=〈(ΔU)²〉/k_BT² in which . The density-scaling exponent may be found either via the statistical-mechanical identity^42,69 γ=〈ΔUΔW〉/〈(ΔU)²〉 or simply by taking the derivative of an analytical approximation to the the function .

As shown in the below subsection ‘The melting-line pressure’, one now has enough information to calculate the pressure along the melting line, p_m(T). To calculate the liquid and solid densities along the melting line (see subsection ‘The freezing- and melting-line densities’ below) one needs to know the below three partial derivatives. Denoting the derivative of the virial along the isomorph with respect to by and recalling that and (refs 42, 69), the three required quantities are given by

The entropy of fusion ΔS_fus is calculated by use of equations (27, 28, 29, 30) below. The three quantities needed here are given by

Isomorph characteristics of generalized LJ pair potentials

The above quantities may be calculated analytically for generalized LJ pair potentials, that is, for systems of particles interacting via pair potential(s) given as a sum or difference of two inverse power laws, r^−m and r⁻ⁿ. The derivation given below applies for any exponents m>n>0 and for general multi-component systems; its subsequent application to freezing and melting deals with single-component systems only.

Invariance of the structure along an isomorph implies that the thermodynamic average potential energy at a given state point, U, may be written (in this section the superscript I is dropped on the reference isomorph density) in which the two A coefficients are functions only of the excess entropy S_ex. For simplicity of notation we shall not indicate the S_ex dependence. The first and second order derivatives of A_m with respect to S_ex are marked by and and likewise for A_n.

The identity for the virial implies

At the reference state point , so for determining A_m and A_n from reference state-point data we have the following two equations:

This implies

From the identity and the definition of the density-scaling exponent, , we get

For determining and from reference state-point data one has

This implies

In order to arrive at equations for and , we first note that , that is, . This implies that . If we define a thermodynamic quantity B by

one has

The two equations for determining and from reference state-point data are thus

This implies

In summary, we have shown that for each of the two reference isomorphs the six numbers A_m, A_n, , , and may be found from reference state-point simulations determining: (1) the potential energy U₀, (2) the virial W₀, (3) the temperature T₀, (4) the excess isochoric specific heat , (5) the density-scaling exponent γ₀ and (6) the derivative of C_V^ex along the isomorph via the quantity B₀ defined in equation (12). The first three quantities are determined directly. The next two quantities are determined from fluctuations at the reference state point: and γ₀=〈ΔWΔU〉/〈(ΔU)²〉. Finally, the quantity B₀ is most accurately found from simulations along the reference isomorph carried out close to the reference state point, although in principle B₀ can be calculated from fluctuations at the reference state point (those needed are of third order and consequently of considerable numerical uncertainty). We calculated B₀ numerically by directly applying equation (12); alternatively, following the methods used in ref. 70 one may rewrite B as B=(γT/C_V^ex)[1+(∂ ln γ/∂ ln T)_ρ] and evaluate B₀ from the (rather weak) constant-density temperature variation of γ at the reference state point.

The melting-line pressure

In the temperature–pressure phase diagram the freezing and melting lines are identical. This section shows how to calculate the pressure on this line as a function of temperature, p_m(T), which is determined by equating the liquid and solid phase’s Gibbs free energies. Recalling that V=(∂G/∂p)_T we estimate these from the Gibbs free energies along the isomorphs, G_l^I(T) and G_s^I(T), as follows (below F_l^I(T) is the Helmholtz free energy along the liquid reference isomorph and likewise for the solid)

The coexistence condition G_l(T,p_m)=G_s(T,p_m) leads to

If F_id is the ideal-gas Helmholtz free energy, the Helmholtz free energy along the liquid isomorph is given by

An analogous expression applies for the solid isomorph’s Helmholtz free energy, , of course. The two constants and are not known, but one needs only their difference. This is determined from the equilibrium condition at the reference state point, G_l(T₀, p₀)=G_s(T₀, p₀) as expressed in equation (17), leading, since pV=Nk_BT+W and F_id(T, ρ_l)−F_id(T, ρ_s)=Nk_BT ln(ρ_l/ρ_s), to

The coexistence condition, equation (17), thus becomes (dropping the explicit temperature dependencies)

or, in terms of the relative density along the respective isomorphs ,

In the case of an arbitrary potential there is no analytical expression for the (average) potential energy as a function of density. Here, the density (of each phase) is the control parameter and T is identified from equation (3), resulting by numerical inversion in two functions and . In the case of generalized LJ pair potentials, for a given temperature T the functions and are found by solving equation (9) (in general numerically, but analytically for the 12-6 LJ system), using the A′ coefficients of equation (11). The potential energy along the isomorphs is given by equation (6).

The freezing- and melting-line densities

We work from the respective reference isomorphs knowing as functions of temperature the coexistence pressure, and the pressure along the reference isomorphs. From this information one calculates the solid and liquid densities by moving on an isotherm from the reference isomorph to the freezing/melting line (Fig. 2). In both cases we define the isothermal difference DW≡W(T)−W^I(T). Here and thoughout the paper D refers to isothermal differences between the reference isomorph and the freezing/melting line.

At any given temperature T the density of the liquid/solid at coexistence is calculated from

If is known, we can determine from equation (22).

The following standard identity is used

In the case of an arbitrary potential, the three terms on the right hand side are calculated from equation (4). For the generalized LJ case, these terms are expressed in terms of the A coefficients by making use of equations (6) and (9), resulting in

We thus have in the generalized LJ case

In both the arbitrary potential case and that of generalized LJ systems, the equation for the density ρ(T)=N/V(T) is found from equation (22) solved numerically in the form

The entropy of fusion

In this section we calculate the constant-pressure entropy of fusion ΔS_fus. One way to do this is to use the Clausius–Clapeyron equation dp_m/dT=ΔS_fus/ΔV_fus in which we now know all quantities except ΔS_fus. An alternative method similar to the above proceeds as follows. Across the melting line one has ΔG_fus=0, that is, ΔE_fus−TΔS_fus+p_m(T)ΔV_fus=0 (E is the total energy). Since the kinetic energy is the same for liquid and solid at the given temperature T, one has ΔE_fus=ΔU_fus and thus

This equation is used for evaluating ΔS_fus from interface-pinning simulations. It is also used for predicting ΔS_fus(T) by proceeding as follows. We have predictions for p_m=p_m(T) and for ΔV_fus=V_l(T)−V_s(T). The missing term is ΔU_fus=ΔU_fus(T), which is estimated via

The partial derivatives refer to the respective reference isomorph as in the last section, and these are evaluated like those of W. Thus,

In the case of an arbitrary potential, the three terms on the right hand side are calculated from equation (5). For the generalized LJ case, these terms may be expressed in terms of the A coefficients of the reference isomorph as follows

We now have all information required for calculating the entropy of fusion.

Melting-line variation of isomorph invariants

We finally turn to the problem of evaluating how much an isomorph invariant X—in casu the reduced vibrational crystalline mean-square displacement, the reduced liquid-state diffusion constant, and the reduced liquid-state viscosity—varies along the freezing/melting line. The starting point is that

On the one hand

On the other hand we have the standard fluctuation formula

Combining these equations at the reference state point leads to (where subscript 0 denotes an equilibrium average at the reference state point)

Consider next an arbitrary temperature T on the freezing/melting line. If DS_ex is the difference between crystal (respectively) liquid excess entropy at melting and that of the corresponding reference isomorph at the same temperature and Dρ likewise is the difference between crystal (respectively) liquid density at melting and that of the corresponding reference isomorph, we estimate X via

Equation (4) implies

Thus we have

This implies

in which the partial derivative is evaluated at the reference state point. If X is a thermodynamic quantity, one may use the fluctuation expression, equation (34), to rewrite this as follows

This expression may be used in the case of an arbitrary potential, as well as for generalized LJ systems for which analytical expressions are available.

Data availability

The data presented in this study are available from the corresponding author upon request.