Thermodynamics of freezing and melting

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.

M elting 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 9 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: u LJ (r) ¼ 4e((r/s) À 12 À (r/s) À 6 ). 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 20 . 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-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 3 that the entropy of fusion DS 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 0 ¼ 2.0e/k B , if the structure were invariant.
In order to develop a quantitative theory of freezing and melting, we take as starting point the 'hidden scale invariance' property of systems 38 characterized by a potential-energy function U(R), where R ¼ (r 1 , r 2 ,y,r N ) is the collective coordinate of the system's N particles, which to a good approximation obeys the scaling condition. 39 UðR a ÞoUðR b Þ ) UðlR a ÞoUðlR b Þ: ð1Þ  37 showing results from T ¼ 0.7e/k B , which is close to the triple point, to T ¼ 3.4e/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 0 ¼ 2.0e/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.
Here, l is a scaling factor and it is understood that the sample container undergoes the same scaling as the configuration; thus l41 corresponds to a density decrease and lo1 to a density increase. This form of scale invariance is exact only for systems with Euler-homogeneous interactions (plus a constant) 13 . 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 )3U(lR a ) ¼ U(lR b ), implies invariance of structure and dynamics along the configurational adiabats in the phase diagram 39 . These lines are referred to as isomorphs 42 . It was very recently shown by Maimbourg and Kurchan 43 that in high dimensions all pairpotential 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 glassforming mixtures 47 , nanoflows 48 , molecular models 38 and molecular dynamics (MD) simulations of the dynamics of most metallic elements based on quantum-mechanical, densityfunctional-theory potentials 49 . Isomorphs have also been demonstrated in simulations of out-of-equilibrium situations like zero-temperature shear flows of glasses or nonlinear steadystate 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 41 . 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 38 . 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 16 . The theory presented below makes use of R simple systems' almost onedimensional phase diagrams 38 and gives corrections to the hardsphere picture of melting and freezing calculated by the firstorder Taylor expansions. Figure 2 illustrates the idea.
Along an isomorph the structure is invariant in the reduced-unit system defined 42 by the length unit r À 1/3 (rN/V is the number density and V is the system volume), the energy unit k B T (T is the temperature) and the time unit ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi mr À 2=3 =k B T p (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 32 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 1 0 ,y, R n 0 . Scaling each of these uniformly to density r one obtains configurations representative for the state point with density r and temperature T on the isomorph through the reference state point 39 in which T is identified from the configurational temperature expression 63 k B T ¼ h(rU) 2 i/hr 2 Ui. The average potential energy U and virial W at the state point (r, 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) 39 (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 ¼ h(DU) 2 i/k B T 2 ) and the so-called density-scaling exponent g ð@ ln T=@ ln rÞ S ex also calculated from the fluctuations (g ¼ hDUDWi/h(DU) 2 i), 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 gðrÞ is isomorph invariant ðr ¼ r 1=3 rÞ. Consequently, for pairs of LJ particles at distance r the thermal average hr À n i scales with density as r n/3 along an isomorph. Thus hr À n i p r n/3 with a proportionality constant that only depends on S ex , implying that the average potential energy U is of the form 64 U ¼ A 12 ðS ex Þr 4 þ A 6 ðS ex Þr 2 in whichr is the density relative to the reference state-point density and A 6 (S ex )o0 derives from the attractive term of the LJ potential. 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 A 0 6 S ex ð Þ 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 12 (S ex ) and A 6 (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 B T þ W in which the virial is given 42 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 0 ¼ 2.0e/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 (qG/qp) 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 byr I , which for temperature T is found by solving T ¼ A 0 12 ðS ex Þðr I Þ 4 þ A 0 6 ðS ex Þðr I Þ 2 . 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 1 (T) þ C 2 (T) À C 3 (T) in which C 1 (T) is the difference between U s I (T) À (T/T 0 )U s I (T 0 ) 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 2 (T) is the difference between Nk B T ln ðr I s ðTÞÞ and the analogous liquid term and C 3 (T) is the difference between (T/T 0 )W s I (T 0 ) 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 37 . 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 (r,T) phase diagram based on equation (26) to that obtained by the interfacepinning MD simulations. Finally, Fig. 3d shows the predicted and simulated fusion entropy DS fus and enthalpy DH 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. 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)e/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)e/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, DGD(T m À T)DS fus (DS 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.

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 16,22,[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 67 . 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 68 . 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 49 .

Methods
Computer simulations. We studied a LJ system of N ¼ 5,000 particles with pair interactions truncated and shifted at 6s. Coexistence pressures, p m , are computed with the interface-pinning method 37 in which coexistence state points are determined by computing the thermodynamic driving force on a solid-liquid interface. Table 1 lists the energy U 0 and virial W 0 at the reference temperature T 0 ¼ 2e/k B for both the liquid and crystal states at coexistence. The A 12 and A 6 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, A 0 12 and A 0 6 , are computed from reference state-point data using equation (11) with the g 0 '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, A 00 12 and A 00 6 , are given by equation (15) where the reference state point excess heat capacity C ex V;0 and B 0 @ðT=C ex V Þ=@ ln r À Á Sex 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 DS 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 0 and its temperature derivative along an isochore, (qL/qT) r , 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 Z (Fig. 4c) are found by letting X ¼D ¼ Dr 1=3 ffiffiffiffiffiffiffiffiffiffiffiffiffiffi m=k B T p and X ¼Z ¼ Z=r 2=3 ffiffiffiffiffiffiffiffiffiffiffi ffi mk B T p , respectively. D is determined from the long-time limit of the mean-square displacement; Z 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.
We proceed to describe the theory in detail. The reference state point is selected at coexistence, that is, with known temperature T 0 and pressure p 0 . There are two different reference densities, a solid and a liquid one, below denoted, respectively, by r s,0 and by r 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 r 0 and temperature T 0 (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 1 0 ,R 2 0 ,y,R n 0 are sampled. In order to map out the reference isomorph parametrized by density, one first identifies the temperature T such that (r,T) is on the isomorph through the reference state point (r 0 ,T 0 ). This is done as follows. If the configurations scaled uniformly to density r are denoted by R 1 ,R 2 ,y,R n in which R i ¼ (r 0 /r) 1/3 R i 0 , 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 Tðr I Þ where we define the relative density along the isomorph byr I r I =r 0 with superscript I indicating 'isomorph' (thus T(1) ¼ T 0 ). By averaging the potential energy U(R) and the virial W(R)( À 1/3)R?rU(R) over the scaled configurations one identifies the functions Uðr I Þ and Wðr I Þ. C ex V ðr I Þ is found from the scaled configurations' potential energy via C V ex ¼ h(DU) 2 i/k B T 2 in which T ¼ Tðr I Þ. The density-scaling exponent gðr I Þ ð@lnT=@lnrÞ Sex may be found either via the statistical-mechanical identity 42,69 g ¼ hDUDWi/h(DU) 2 i or simply by taking the derivative of an analytical approximation to the the function TðrÞ.
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 tor I by W 0 ðr I Þ and recalling that W ¼ ð@U=@ lnrÞ Sex and T ¼ ð@U=@S ex Þr (refs 42,69), the three required quantities are given by The entropy of fusion DS 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 À n . The derivation given below applies for any exponents m4n40 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 U ¼ A mr m=3 þ A nr n=3 (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 A 0 m and A 00 m and likewise for A n . The identity for the virial W ¼ ð@U=@ lnrÞ Sex implies At the reference state pointr ¼ 1, so for determining A m and A n from reference state-point data we have the following two equations: This implies From the identity T ¼ ð@U=@S ex Þr and the definition of the density-scaling exponent, g ð@ ln T=@ lnrÞ Sex , we get For determining A 0 m and A 0 n from reference state-point data one has

This implies
In order to arrive at equations for A 00 m and A 00 n , we first note that nr n=3 . If we define a thermodynamic quantity B by one has The two equations for determining A 00 m and A 00 n 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 , A 0 m , A 0 n , A 00 m and A 00 n may be found from reference state-point simulations determining: (1) the potential energy U 0 , (2) the virial W 0 , (3) the temperature T 0 , (4) the excess isochoric specific heat C ex V;0 , (5) the density-scaling exponent g 0 and (6) the derivative of C V ex along the isomorph via the quantity B 0 defined in equation (12). The first three quantities are determined directly. The next two quantities are determined from fluctuations at the reference state point: Finally, the quantity B 0 is most accurately found from simulations along the reference isomorph carried out close to the reference state point, although in principle B 0 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 0 numerically by directly applying equation (12); alternatively, following the methods used in ref. 70 one may rewrite B as B ¼ (gT/C V ex )[1 þ (q ln g/q ln T) r ] and evaluate B 0 from the (rather weak) constant-density temperature variation of g 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 ¼ (qG/qp) 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) If F id is the ideal-gas Helmholtz free energy, the Helmholtz free energy along the liquid isomorph is given by F I l ðTÞ ¼ U I l ðTÞ À TS I ex;l þ F id ðT; r I l ðTÞÞ: ð18Þ An analogous expression applies for the solid isomorph's Helmholtz free energy, F I s T ð Þ, of course. The two constants S I ex;l and S I ex;s are not known, but one needs only their difference. This is determined from the equilibrium condition at the reference state point, G l (T 0 , p 0 ) ¼ G s (T 0 , p 0 ) as expressed in equation (17), leading, since pV ¼ Nk B T þ W and F id (T, r l ) À F id (T, r s ) ¼ Nk B T ln(r l /r s ), to T 0 ðS I ex;l À S I ex;s Þ ¼ ðU l;0 À U s;0 Þ þ Nk B T 0 ln ðr l;0 =r s;0 Þ þ ðW l;0 À W s;0 Þ: ð19Þ The coexistence condition, equation (17), thus becomes (dropping the explicit temperature dependencies) p m ðV I l À V I s Þ ¼ðU I s À U I l Þ À T T 0 ððU s;0 À U l;0 Þ þ Nk B T 0 lnðr s;0 =r l;0 Þ þ ðW s;0 À W l;0 ÞÞ þ Nk B Tlnðr I s =r I l Þ or, in terms of the relative density along the respective isomorphsr I , 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 functionsr I s ðTÞ andr I l ðTÞ. In the case of generalized LJ pair potentials, for a given temperature T the functionsr I l ðTÞ andr I s ðTÞ are found by solving equation (9) (in general numerically, but analytically for the 12-6 LJ system), using the A 0 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 DWW(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 densityr of the liquid/solid at coexistence is calculated from If ð@W=@ lnrÞ I T is known, we can determiner from equation (22). The following standard identity is used @W @lnr 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 @W @ lnr I In both the arbitrary potential case and that of generalized LJ systems, the equation for the density r(T) ¼ N/V(T) is found from equation (22) solved numerically in the form p m ðTÞVðTÞ À Nk B T À W I ðTÞ ¼ @W @ lnr I T lnðr=r I Þ: ð26Þ The entropy of fusion. In this section we calculate the constant-pressure entropy of fusion DS fus . One way to do this is to use the Clausius-Clapeyron equation dp m /dT ¼ DS fus /DV fus in which we now know all quantities except DS fus . An alternative method similar to the above proceeds as follows. Across the melting line one has DG fus ¼ 0, that is, DE fus À TDS fus þ p m (T)DV 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 DE fus ¼ DU fus and thus This equation is used for evaluating DS fus from interface-pinning simulations. It is also used for predicting DS fus (T) by proceeding as follows. We have predictions for p m ¼ p m (T) and for DV fus ¼ V l (T) À V s (T). The missing term is DU fus ¼ DU 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.

ARTICLE
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 @U @lnr I 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 @X @T 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 Dr likewise is the difference between crystal (respectively) liquid density at melting and that of the corresponding reference isomorph, we estimate X via X ffi X 0 þ f 0 ðS ex ÞDS ex ffi X 0 þ f 0 ðS ex Þ @S ex @r T Dr: Equation (4) implies Thus we have This implies X ffi X 0 À @X @T 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.