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.

ARTICLE

Received 25 Apr 2016 | Accepted 27 Jun 2016 | Published 17 Aug 2016

DOI: 10.1038/ncomms12386 OPEN

Thermodynamics of freezing and melting

Ulf R. Pedersen1, Lorenzo Costigliola1, Nicholas P. Bailey1, Thomas B. Schrder1 & Jeppe C. Dyre1

Although the freezing of liquids and melting of crystals are fundamental for many areas of the sciences, even simple properties like the temperaturepressure 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 liquids 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.

1 Glass and Time, IMFUFA, Department of Science and Environment, Roskilde University, P. O. Box 260, Roskilde DK-4000, Denmark. Correspondence and requests for materials should be addressed to U.R.P. (email:mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 7:12386 | DOI: 10.1038/ncomms12386 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12386

Melting is the prototypical rst-order phase transition13. Its qualitative description has been textbook knowledge for a century, but it has proven difcult 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 research48; 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 compress9 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: uLJ(r) 4e((r/s) 12 (r/s) 6). The rst term reects 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 liquids1116. Their main ingredient is the assumption that the structure of a dense, monatomic uid closely resembles that of a collection of hard spheres14,1618. Conrming this, the structure of melts of, for example, metallic elements near freezing is close to that of the hard-sphere system15,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 entropy20. 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 distanceknown as the Lindemann ratiois constant and about0.1 along the melting line; this is the famous Lindemann melting criterion from 1910 (refs 2025). 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 liquids self-diffusion constant and viscosity are invariant along the freezing line26,27, the HansenVerlet rule17,28 that the amplitude of the rst peak of the liquid static structure factor is about 2.85 at freezing, or Richards melting rule3 that the entropy of fusion DSfus is about 1.1kB (which in a more modern and accurate version is the fact that the constant-volume entropy difference across the densitytemperature coexistence region is close to 0.8kB (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 quanties the deviations from the above mentioned hard-sphere predicted melting-line invariants16,22,3032. The theory is validated by computer simulations of the standard 12-6 LJ system.

ResultsGeneral theory. It is well-known that adding a mean-eld attractive term to the hard-sphere model broadens the coexistence region, which on the other hand, narrows if the repulsive part is softened13,16,3336. 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 T0 2.0e/kB, 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 systems38 characterized by a potential-energy function U(R), where R (r1, r2,y,rN) is the collective

coordinate of the systems N particles, which to a good approximation obeys the scaling condition.39

URaoURb ) UlRaoUlRb: 1

a

Freezing line

b

Isomorph

3.2

1.4 2.0 2.4 3.0 3.4

3.2

0.7 1.0 1.4 2.0 2.4 3.0 3.4

4

4

Structure factor

Structure factor

3

3

2

2

1

1

0 0 5 10 15 q (1/[afii9846])

0 0 5 10 15 q (1/[afii9846])

Figure 1 | Structure of the LJ liquid. (a) Liquid structure factor along the freezing line37 showing results from T 0.7e/kB, which is close to the triple

point, to T 3.4e/kB. The hard-sphere model predicts that the height of the rst 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 temperatureT0 2.0e/kB (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.

2 NATURE COMMUNICATIONS | 7:12386 | DOI: 10.1038/ncomms12386 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12386 ARTICLE

Freezing/melting lines Liquid/crystal isomorphs

2

Temperature

p

mr 2=3=kBT

(m is the particle mass). Figure 1b shows the LJ liquids 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 conrms the recent nding of Heyes and Branka32 that the freezing line is not an exact isomorph, although it is close to one.

The melting pressure as a function of temperature, pm(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 Sex

(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 congurations R10,y, Rn0. Scaling each of these uniformly to density r one obtains congurations representative for the state point with density r and temperature T on the isomorph through the reference state point39 in which T is identied from the congurational temperature expression63 kBT h(rU)2i/hr2Ui. The average potential energy U and virial

W at the state point (r, T) are likewise found by averaging over the scaled congurations. The key assumption here is that the canonical probabilities of the scaled congurations are identical to those of the original congurations, 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 specic heat CVex calculated from the potential-energy uctuations of the scaled congurations (CVex h(DU)2i/kBT2)

and the so-called density-scaling exponent g @ ln T=@ ln rSex

also calculated from the uctuations (g hDUDWi/h(DU)2i), 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 r1=3r. Consequently,

for pairs of LJ particles at distance r the thermal average hr ni

scales with density as rn/3 along an isomorph. Thus hr ni p rn/3

with a proportionality constant that only depends on Sex,

implying that the average potential energy U is of the form64 U A12Sex~r4 A6Sex~r2 in which ~r is the density

relative to the reference state-point density and A6(Sex)o0

derives from the attractive term of the LJ potential. Since T (qU/qSex)r, one has T A012Sex~r4 A06Sex~r2. It follows

that if the ve quantities Sex, A12(Sex), A6(Sex), A012 Sex

and

1

0.9 1.0 1.1 Density

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 rst-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 T0 2.0e/kB. Details are given in the Methods section.

Here, l is a scaling factor and it is understood that the sample container undergoes the same scaling as the conguration; 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 invariance3942. Equation (1), which is formally equivalent to the conformal-invariance condition U(Ra) U(Rb)3U(lRa)

U(lRb), implies invariance of structure and dynamics along the congurational adiabats in the phase diagram39. These lines are referred to as isomorphs42. It was very recently shown by Maimbourg and Kurchan43 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 polymers4446. Further evidence for the existence of isomorphs comes from computer simulations of single-component systems40,42 as well as, for example, of glass-forming mixtures47, nanoows48, molecular models38 and molecular dynamics (MD) simulations of the dynamics of most metallic elements based on quantum-mechanical, density-functional-theory potentials49. Isomorphs have also been demonstrated in simulations of out-of-equilibrium situations like zero-temperature shear ows of glasses or nonlinear steady-state liquid ows (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 exception41. 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 not38. 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) simple35,5062 to distinguish them from simple systems traditionally dened as pair-potential systems16. The theory presented below makes use of R simple systems almost one-dimensional phase diagrams38 and gives corrections to the hard-sphere picture of melting and freezing calculated by the rst-order Taylor expansions. Figure 2 illustrates the idea.

Along an isomorph the structure is invariant in the reduced-unit system dened42 by the length unit r 1/3 (r N/V is the number density and V is the system volume),

the energy unit kBT (T is the temperature) and the time unit

NATURE COMMUNICATIONS | 7:12386 | DOI: 10.1038/ncomms12386 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12386

are known, the excess Helmholtz free energy, UTSex, 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-points potential energy and virial give two linear equations for determining A12(Sex) and A6(Sex). 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 NkBT W in which the virial is given42 by

W @U=@ ln ~rSex 4A12Sex~r4 2A6Sex~r2).

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 T0 2.0e/kB, a liquid-phase isomorph and

a crystal-phase isomorph. Gibbs free energy of the liquid phase at coexistence, Gl(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 Gl at coexistence is

GlpmT; T GIlT VIlTpmT pIlT: 2 Here, pm(T) is the coexistence pressure to be determined; GlI(T) is the Gibbs-free energy, VlI(T) the volume and plI(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 ~rI, which for temperature T is found by solving T A012Sex~rI4 A06Sex~rI2.

An analogous expression applies for the crystals 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 pm(T)(VlI(T) VsI(T)) C1(T) C2(T) C3(T) in

which C1(T) is the difference between UsI(T) (T/T0)UsI(T0) and

the analogous term for the liquid reference isomorph (here UsI(T)

is the crystals potential energy along the reference isomorph), C2(T) is the difference between NkBT ln ~rIsT and the

analogous liquid term and C3(T) is the difference between (T/T0)WsI(T0) and the analogous liquid term.

Figure 3a,b compare the theoretically predicted pm(T) to the

coexistence pressure computed numerically by means of the interface-pinning method37. The density of the crystalline and liquid phases may also be computed by means of a rst-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 interface-pinning MD simulations. Finally, Fig. 3d shows the predicted and simulated fusion entropy DSfus and enthalpy DHfus, 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.

A06 Sex

a

Reference point

Theoretical prediction

100

b

50

40

10

p m([afii9830]/[afii9846]3 )

p m([afii9830]/[afii9846]3 )

30

20

1

10

0.1

0 1 2 3 4 T ([afii9830]/kB)

0 0 1 2 3 4 T ([afii9830]/kB)

c d

1.6

3

1.4

T ([afii9830]/k B)

2

S fus(k B)

1.2

1

0.9 1.0 1.1 1.2

[afii9845] (1/[afii9846]3)

1 1 2 3 4T ([afii9830]/kB)

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 gure (T0 2.0e/kB), the MD simulations employed the

interface-pinning method37, see the Methods section. No tting was done in these guresthe only input to the theory is properties of the coexisting liquid and crystal at the reference temperature. (a) Temperaturepressure phase diagram. The green dashed line marks the expectation based on a linear extrapolation of the ClausiusClapeyron relation dpm/dT DSfus/DVfus 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 densitytemperature phase diagram; the coloured area is the coexistence region. (d) Fusion entropy (main panel) and enthalpy (inset).

4 NATURE COMMUNICATIONS | 7:12386 | DOI: 10.1038/ncomms12386 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12386 ARTICLE

a b

0.14

Theoretical prediction Reducedunit invariant

Diffusion constant along

0.05

Lindemann ratio along

metling line

freezing line

0.04

0.13

0.03

0.12

1 2 3 4T ([afii9830]/kB)

1 2 3 4T ([afii9830]/kB)

c

12

Viscosity along freezing line

8

4

1 2 3 4T ([afii9830]/kB)

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.

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 Ttp 0.688(2)e/kB 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 ClausiusClapeyron relation from the reference temperature (the green dashed lines in Fig. 3a,b) predicts a triple point temperature of 0.909(2)e/kB.

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 rst-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 WilsonFrenkel law65,66. This motivated us to use the theory also for calculating the liquids diffusion constant variation along the freezing line (Fig. 4b). Another important component for crystal growth is the thermodynamic driving force on the crystalliquid interface, which is the Gibbs free energy difference between the two phases, DGD(Tm T)DSfus (DSfus 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.

DiscussionThe 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 picture16,22,3032. 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 rm 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 system67. 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 uid of metallic hydrogen. The quantum nature of the proton modies classical melting, for example by increasing the value of the Lindemann ratio68. 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

NATURE COMMUNICATIONS | 7:12386 | DOI: 10.1038/ncomms12386 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12386

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

T0 2.0 e/kB Liquid Crystal

V0/N [s3] 0.9403(2)* 0.8827(2) r0 [s 3] 1.0633(2) 1.1329(2) U0/N [e] 4.7792(2) 5.7774(2)

W0/N [e] 17.5418(7) 16.3628(6) g0 4.9164(8) 4.8704(8) CexV;0/N [kB] 1.323(5) 1.301(7)

B0N [e/kB2] 6.9(5) 7.2(5) L0 0.1398(2) (qL/qT)r [kB/e] 0.041(2)

~D0 0.02921(9)

@~D=@Tr [kB/e] 0.0201(4)

0 5.2487(6)

@=@Tr [kB/e] 2.60(14)

These numbers were used for calculating the theoretical predictions in Figs 3 and 4 (red curves). *Numbers in parenthesis give the estimated statistical uncertainty using a 95% condence interval.

This determines the function T~rI where we dene the relative density

along the isomorph by ~rI rI=r0 with superscript I indicating isomorph

(thus T(1) T0). By averaging the potential energy U(R) and the virial

W(R) ( 1/3)R?rU(R) over the scaled congurations one identies the functions

U~rI and W~rI. CexV~rI is found from the scaled congurations potential

energy via CVex h(DU)2i/kBT2 in which T T~rI. The density-scaling exponent

g~rI @lnT=@lnrS may be found either via the statistical-mechanical

identity42,69 g hDUDWi/h(DU)2i 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, pm(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 ~rI by W0~rI and recalling that W @U=@ ln ~rS and

T @U=@Sex~r (refs 42,69), the three required quantities are given by@W @ ln ~r

IS

~rIW0~rI

@W @Sex

I~r

@2U

@Sex@ ln ~r

@T @ ln ~r

S T

~rIg~rI

4

The entropy of fusion DSfus is calculated by use of equations (2730) below. The three quantities needed here are given by

@U @ ln ~r

IS

W~rI

@U @Sex

I~r

@Sex

@ ln ~r

IT

@T @ ln ~r

S

@T @S

~r

T~rIg~rI T~r

I

=CexV~rI

CexV~rIg~rI:

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

Methods

Computer simulations. We studied a LJ system of N 5,000 particles with pair

interactions truncated and shifted at 6s. Coexistence pressures, pm, are computed with the interface-pinning method37 in which coexistence state points are determined by computing the thermodynamic driving force on a solid-liquid interface. Table 1 lists the energy U0 and virial W0 at the reference temperature T0 2e/kB for both the liquid and

crystal states at coexistence. The A12 and A6 coefcients (for the liquid and the crystal separately) are computed from reference state-point data using equation (8) below. The derivatives of the A coefcients with respect to excess entropy, A012 and A06, are computed from reference state-point data using equation (11) with the g0s 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 coefcients, A0012 and A006, are given by equation (15) where the reference state point excess heat capacity CexV;0 and

B0 @T=CexV=@ ln r

T~rI

5

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 Am~rm=3 An~rn=3 (in this section the superscript I is dropped on the reference

isomorph density) in which the two A coefcients are functions only of the excess entropy Sex. For simplicity of notation we shall not indicate the Sex dependence.

The rst and second order derivatives of Am with respect to Sex are marked by A0m and A00m and likewise for An.

The identity for the virial W @U=@ ln ~rS impliesU Am~rm=3 An~rn=3

W

m

3 Am~rm=3

n

@Sex

@ ln ~r

IT CexV

~rIg~rI:

S 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 DSfus (Fig. 3d) is computed by combining equations (2730).

The value of the Lindemann ratio L of the crystal at the reference temperature, L0 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 Dr1=3

p

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 T0 and pressure p0. There are two different reference densities, a solid and a liquid one, below denoted, respectively, by rs,0 and by rl,0. In the densitytemperature 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 temperaturepressure reference state point denes two reference densitytemperature state points, a liquid and a solid one. Let us focus on one of these with density r0 and temperature T0 (thus dropping in this and the next subsection subscripts s and l). From an NVT MD equilibrium simulation(for example, with a NosHoover thermostat) n congurations R10,R20,y,Rn0 are

sampled. In order to map out the reference isomorph parametrized by density, one rst identies the temperature T such that (r,T) is on the isomorph through the reference state point (r0,T0). This is done as follows. If the congurations scaled uniformly to density r are denoted by R1,R2,y,Rn in which Ri (r0/r)1/3Ri0,

the temperature T is determined from the standard congurational temperature expression (in which the averages are over the n sampled congurations)

kBT hr

m=kBT

and X Z=r2=3

mkBT

3 An~rn=3:

6

At the reference state point ~r 1, so for determining Am and An from reference

state-point data we have the following two equations:

U0 Am An

W0

m

3 Am

n

3 An:

7

This implies

Am

3W0 nU0 m n

An

8

mU0 3W0 m n

:

From the identity T @U=@Sex~r and the denition of the density-scaling

exponent, g @ ln T=@ ln ~rS , we get

T A0m~rm=3 A0n~rn=3

gT

m

3 A0m~rm=3

n

3 A0n~rn=3:

9

URi2ii hr

2URiii

: 3

For determining A0m and A0n from reference state-point data one has

T0 A0m A0n

g0T0

m

3 A0m

n

3 A0n:

10

6 NATURE COMMUNICATIONS | 7:12386 | DOI: 10.1038/ncomms12386 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12386 ARTICLE

This implies

A0m

3g0 n m n

T0

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

pmVIl VIs UIs

TT0 Us;0

UIl

TT0 Ul;0

NkBTln~rIs=~rIl

A0n

m 3g0

m n

T0:

11

TT0 Wl;0 Ws;0:

; 12

one has

TCexV A00m~rm=3 A00n~rn=3

B

m

3 A00m~rm=3

n

3 A00n~rn=3:

21

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 identied from equation (3), resulting by numerical inversion in two functions ~rIsT and ~rIlT. In the case of generalized LJ pair potentials, for a given

temperature T the functions ~rIlT and ~rIsT are found by solving equation (9)

(in general numerically, but analytically for the 12-6 LJ system), using the A0 coefcients 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 dene the isothermal difference DW W(T) WI(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 ~r of the liquid/solid at coexistence is calculated from

DW

@W @ln~r

IT Dln

In order to arrive at equations for A00m and A00n, we rst note that CexV T@Sex=@T~r, that is, @T=@Sex~r T=CexV. This implies that

T=CexV A00m~rm=3 A00n~rn=3. If we dene a thermodynamic quantity B by

B

@T=CexV @ln~r

S

13

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

T0CexV;0 A00m A00n

B0

m

3 A00m

14

n

3 A00n:

~r=~rI: 22

If @W=@ ln ~rIT is known, we can determine ~r from equation (22).

The following standard identity is used

@W @ln~r

T

~r

@W @ln~r

IT ln

This implies

A00m

3B0 nT0=CexV;0

m n A00n

mT0=CexV;0 3B0 m n

:

@W @ln~r

S

@W @Sex

~r

@Sex

@ln~r

T

: 23

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 coefcients by making use of equations (6) and (9), resulting in

@W @ ln ~r

IS

m

3

2

15

In summary, we have shown that for each of the two reference isomorphs the six numbers Am, An, A0m, A0n, A00m and A00n may be found from reference state-point simulations determining: (1) the potential energy U0, (2) the virial W0, (3) the temperature T0, (4) the excess isochoric specic heat CexV;0, (5) the density-scaling exponent g0 and (6) the derivative of CVex along the isomorph via the quantity B0 dened in equation (12). The rst three quantities are determined directly. The next two quantities are determined from uctuations at the reference state point: CexV;0 hDU2i=kBT20 and g0 hDWDUi/h(DU)2i. Finally, the quantity B0 is most

accurately found from simulations along the reference isomorph carried out close to the reference state point, although in principle B0 can be calculated from uctuations at the reference state point (those needed are of third order and consequently of considerable numerical uncertainty). We calculated B0 numerically by directly applying equation (12); alternatively, following the methods used in ref. 70 one may rewrite B as B (gT/CVex)[1 (q ln g/q ln T)r] and evaluate B0 from

the (rather weak) constant-density temperature variation of g at the reference state point.

The melting-line pressure. In the temperaturepressure 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, pm(T), which is determinedby equating the liquid and solid phases Gibbs free energies. Recalling thatV (qG/qp)T we estimate these from the Gibbs free energies along the isomorphs,

GlI(T) and GsI(T), as follows (below FlI(T) is the Helmholtz free energy along the liquid reference isomorph and likewise for the solid)

GlT; pmT GIlT VIlTpmT pIlT FIlT VIlTpmT GsT; pmT GIsT VIsTpmT pIsT FIsT VIsTpmT:

16

Am~rIm=3

n

3

2

An~rIn=3

@W @Sex

I~r

m

3 A0m~rIm=3

n

3 An~rIn=3

24

@Sex

@ ln ~r

IT

@T @ ln ~r

IS @T @S

I~r

m

3 A0m~rIm=3 n3 A0n~rIn=3 A00m~r

I

m=3 A00n~r

I

n=3

:

We thus have in the generalized LJ case

@W @ ln ~r

IT

m

3

2

Am~rIm=3

n

3

2

An~rIn=3

m

3 A0m~rIm=3 n3 A0n~rIn=3

2

A00m~r

I

:

25

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

pmTVT NkBT WIT

@W @ ln ~r

IT ln

m=3 A00n~r

I

n=3

~r=~rI: 26

The entropy of fusion. In this section we calculate the constant-pressure entropy of fusion DSfus. One way to do this is to use the ClausiusClapeyron equation

dpm/dT DSfus/DVfus in which we now know all quantities except DSfus. An

alternative method similar to the above proceeds as follows. Across the melting line one has DGfus 0, that is, DEfus TDSfus pm(T)DVfus 0 (E is the total energy).

Since the kinetic energy is the same for liquid and solid at the given temperature T, one has DEfus DUfus and thus

DSfus

The coexistence condition Gl(T,pm) Gs(T,pm) leads to

pmTVIlT VIsT FIsT FIlT: 17 If Fid is the ideal-gas Helmholtz free energy, the Helmholtz free energy along the liquid isomorph is given by

FIlT UIlT TSIex;l FidT; rIlT: 18 An analogous expression applies for the solid isomorphs Helmholtz free energy, FIs T

, of course. The two constants SIex;l and SIex;s are not known, but one needs

only their difference. This is determined from the equilibrium condition at the reference state point, Gl(T0, p0) Gs(T0, p0) as expressed in equation (17), leading,

since pV NkBT W and Fid(T, rl) Fid(T, rs) NkBT ln(rl/rs), to

T0SIex;l SIex;s Ul;0 Us;0 NkBT0 ln rl;0=rs;0 Wl;0 Ws;0: 19 The coexistence condition, equation (17), thus becomes (dropping the explicit temperature dependencies)

pmVIl VIs UIs UIl

TT0 Us;0 Ul;0 NkBT0lnrs;0=rl;0

Ws;0 Wl;0 NkBTlnrIs=rIl

DUfus pmDVfus

T : 27 This equation is used for evaluating DSfus from interface-pinning simulations. It is also used for predicting DSfus(T) by proceeding as follows. We have predictions for pm pm(T) and for DVfus Vl(T) Vs(T). The missing term is DUfus DUfus(T),

which is estimated via

DUfus UIlT

@U @ln~r

I;l

T

ln~rlT=~rIlT UIsT

I;s

T

@U @ln~r

ln~rsT=~rIsT:

28

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

@U @ln~r

IT

@U @ln~r

IS

20

@U @Sex

I~r

@Sex

@ln~r

IT

: 29

NATURE COMMUNICATIONS | 7:12386 | DOI: 10.1038/ncomms12386 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12386

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 coefcients of the reference isomorph as follows

@U @ln~r

IS

WI

m

3 Am~rIm=3

n

3 An~rIn=3

@U @Sex

I~r

T A0m~rIm=3 A0n~rIn=3

@Sex

@ln~r

IT

30

@T @ln~r

IS @T @S

I~r

m

3 A0m~rIm=3 n3 An ~rIn=3 A00m~r

I

m=3 A00n~r

I

n=3

:

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

Melting-line variation of isomorph invariants. We nally turn to the problem of evaluating how much an isomorph invariant Xin casu the reduced vibrational crystalline mean-square displacement, the reduced liquid-state diffusion constant, and the reduced liquid-state viscosityvaries along the freezing/melting line. The starting point is that

X fSex: 31

On the one hand

@X @T

r f0

Sex

r f0

Sex

CexV

T : 32

@Sex

@T

On the other hand we have the standard uctuation formula

@X @T

r

hDXDUi

kBT2 : 33

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

f0Sex h

DXDUi0

kBT0CexV;0 : 34

Consider next an arbitrary temperature T on the freezing/melting line. If DSex 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 X0 f0SexDSex X0 f0Sex

@Sex

@r

T

Dr: 35

Equation (4) implies

@Sex

@r

IT

gCexV

r : 36

Thus we have

X X0 f0SexgCexV

Drr : 37

This implies

Drr 38

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

X X0 h

X X0

@X @T

rg

T0 CexV

CexV;0

Drr : 39

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.

References

1. Chandler, D. Introduction to Modern Statistical Mechanics (Oxford University Press, 1987).

2. Atkins, P. W. Physical Chemistry 4th edn (Oxford Univ. Press, 1990).3. Glicksman, M. E. Principles of Solidication: An Introduction to Modern Casting and Crystal Growth Concepts (Springer, 2011).

4. Eggert, J. H. et al. Melting temperature of diamond at ultrahigh pressure. Nat. Phys. 6, 4043 (2010).

5. Guillaume, C. L. et al. Cold melting and solid structures of dense lithium. Nat. Phys. 7, 211214 (2011).

6. Peng Tan, N. X. & Xu, L. Visualizing kinetic pathways of homogeneous nucleation in colloidal crystallization. Nat. Phys. 10, 7379 (2014).

7. Deutschlander, S., Puertas, A. M., Maret, G. & Keim, P. Specic heat in two-dimensional melting. Phys. Rev. Lett. 113, 127801 (2014).

8. Statt, A., Virnau, P. & Binder, K. Finite-size effects on liquid-solid phase coexistence and the estimation of crystal nucleation barriers. Phys. Rev. Lett. 114, 026101 (2015).

9. van der Waals, J. D. On the Continuity of the Gaseous and Liquid States. Ph.D. thesis, Universiteit Leiden (1873).

10. Lennard-Jones, J. E. On the determination of molecular elds. I. From the variation of the viscosity of a gas with temperature. Proc. R. Soc. Lond. A 106, 441462 (1924).

11. Barker, J. A. & Henderson, D. What is liquid? Understanding the states of matter. Rev. Mod. Phys. 48, 587671 (1976).

12. Gubbins, K., Smitha, W., Tham, M. & Tiepel, E. Perturbation theory for the radial distribution function. Mol. Phys. 22, 1089 (1971).

13. Hoover, W. G., Gray, S. G. & Johnson, K. W. Thermodynamic propertiesof the uid and solid phases for inverse power potentials. J. Chem. Phys. 55, 11281136 (1971).

14. Weeks, J. D., Chandler, D. & Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys. 54, 52375247 (1971).

15. Hausleitner, C., Kahl, G. & Hafner, J. Liquid structure of transition metals: investigations using molecular dynamics and perturbation- and integral-equation techniques. J. Phys.: Condens. Mat 3, 1589 (1991).

16. Hansen, J.-P. & McDonald, I. R. Theory of Simple Liquids: With Applications to Soft Matter 4th edn (Academic, 2013).

17. Hansen, J.-P. & Verlet, L. Phase transitions of the Lennard-Jones system. Phys. Rev. 184, 151161 (1969).

18. Rosenfeld, Y. Theory of simple classical uids: universality in the short-range structure. Phys. Rev. A 20, 12081235 (1979).

19. Waseda, Y., Yokoyama, K. & Suzuki, K. Structure of molten Mg, Ca, Sr, and Ba by X-ray diffraction. Z. Naturforsch. A 30, 801805 (1975).

20. Ross, M. Generalized Lindemann melting law. Phys. Rev. 184, 233242 (1969).21. Lindemann, F. A.ber die Berechning molekularer Eigenfrequenzen. Phys. Z. 11, 609 (1910).

22. Lawson, A. C. An improved Lindemann melting rule. Philos. Mag. B 81, 255266 (2001).

23. Stishov, S. M. The thermodynamics of melting of simple substances. Sov. Phys. Usp. 17, 625643 (1975).

24. Luo, S.-N., Strachan, A. & Swift, D. C. Vibrational density of states and Lindemann melting law. J. Chem. Phys. 122, 194709 (2005).

25. Chakravarty, C., Debenedetti, P. G. & Stillinger, F. H. Lindemann measures for the solidliquid phase transition. J. Chem. Phys. 126, 204508 (2007).

26. Andrade., E. N. C. A theory of the viscosity of liquidsPart I. Phil. Mag. 17, 497511 (1934).

27. Costigliola, L., Schrder, T. B. & Dyre, J. C. Freezing and melting line invariants of the Lennard-Jones system. Phys. Chem. Chem. Phys. 18, 1467814690 (2016).

28. Hansen, J.-P. Phase transition of the Lennard-Jones system. II. High-temperature limit. Phys. Rev. A 2, 221230 (1970).

29. Tallon, J. L. The entropy change on melting of simple substances. Phys. Lett. A 76, 139142 (1980).

30. Pedersen, U. R., Hummel, F. & Dellago, C. Computing the crystal growth rate by the interface pinning method. J. Chem. Phys. 142, 044104 (2015).

31. Benjamin, R. & Horbach, J. Crystal growth kinetics in Lennard-Jones and WeeksChandlerAndersen systems along the solidliquid coexistence line.J. Chem. Phys. 143, 014702 (2015).32. Heyes, D. M. & Branka, A. C. The Lennard-Jones melting line and isomorphism. J. Chem. Phys. 143, 234504 (2015).

33. Agrawal, R. & Kofke, D. A. Soliduid coexistence for inverse-power potentials. Phys. Rev. Lett. 74, 122 (1995).

34. Ben-Amotz, D. & Stell, G. Hard sphere perturbation theory for uids with soft-repulsive-core potentials. J. Chem. Phys. 120, 4844 (2004).

35. Heyes, D. M., Dini, D. & Branka, A. C. Scaling of Lennard-Jones liquid elastic moduli, viscoelasticity and other properties along uidsolid coexistence. Phys. Status Solidi (b) 252, 15141525 (2015).

36. Ding, Y. & Mittal, J. Equilibrium and nonequilibrium dynamics of soft sphere uids. Soft Matter 11, 52745281 (2015).

37. Pedersen, U. R. Direct calculation of the solidliquid Gibbs free energy difference in a single equilibrium simulation. J. Chem. Phys. 139, 104102 (2013).

38. Dyre, J. C. Hidden scale invariance in condensed matter. J. Phys. Chem. B 118, 1000710024 (2014).

39. Schrder, T. B. & Dyre, J. C. Simplicity of condensed matter at its core: generic denition of a Roskilde-simple system. J. Chem. Phys. 141, 204502 (2014).40. Pedersen, U. R., Bailey, N. P., Schrder, T. B. & Dyre, J. C. Strong pressureenergy correlations in van der Waals liquids. Phys. Rev. Lett. 100, 015701 (2008).

DXDUi0

kBT0 g

CexV

CexV;0

8 NATURE COMMUNICATIONS | 7:12386 | DOI: 10.1038/ncomms12386 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12386 ARTICLE

41. Bailey, N. P., Pedersen, U. R., Gnan, N., Schrder, T. B. & Dyre, J. C. Pressure-energy correlations in liquids. I. Results from computer simulations.J. Chem. Phys. 129, 184507 (2008).42. Gnan, N., Schrder, T. B., Pedersen, U. R., Bailey, N. P. & Dyre, J. C. Pressure-energy correlations in liquids. IV. Isomorphs in liquid phase diagrams. J. Chem. Phys. 131, 234504 (2009).

43. Maimbourg, T. & Kurchan, T. Approximate scale invariance in particle systems: a large-dimensional justication. EPL 114, 60002 (2016).

44. Roland, C. M., Hensel-Bielowka, S., Paluch, M. & Casalini, R. Supercooled dynamics of glass-forming liquids and polymers under hydrostatic pressure. Rep. Prog. Phys. 68, 14051478 (2005).

45. Gundermann, D. et al. Predicting the density-scaling exponent of a glass-forming liquid from prigogine-defay ratio measurements. Nat. Phys. 7, 816821 (2011).

46. Xiao, W., Tofteskov, J., Christensen, T. V., Dyre, J. C. & Niss, K. Isomorph theory prediction for the dielectric loss variation along an isochrone.J. Non-Cryst. Solids 407, 190195 (2015).47. Pedersen, U. R., Schrder, T. B. & Dyre, J. C. Repulsive reference potential reproducing the dynamics of a liquid with attractions. Phys. Rev. Lett. 105, 157801 (2010).

48. Ingebrigtsen, T. S., Errington, J. R., Truskett, T. M. & Dyre, J. C. Predicting how nanoconnement changes the relaxation time of a supercooled liquid. Phys. Rev. Lett. 111, 235901 (2013).

49. Hummel, F., Kresse, G., Dyre, J. C. & Pedersen, U. R. Hidden scale invariance of metals. Phys. Rev. B 92, 174116 (2015).

50. Malins, A., Eggers, J. & Royall, C. P. Investigating isomorphs with the topological cluster classication. J. Chem. Phys 139, 234505 (2013).

51. Abramson, E. H. Viscosity of uid nitrogen to pressures of 10 GPa. J. Phys. Chem. B 118, 1179211796 (2014).

52. Fernandez, J. & Lopez, E. R. Experimental Thermodynamics: Advances in Transport Properties of Fluids ch. 9.3, 307317 (Royal Society of Chemistry, 2014).

53. Flenner, E., Staley, H. & Szamel, G. Universal features of dynamic heterogeneity in supercooled liquids. Phys. Rev. Lett. 112, 097801 (2014).

54. Prasad, S. & Chakravarty, C. Onset of simple liquid behaviour in modied water models. J. Chem. Phys. 140, 164501 (2014).

55. Buchenau, U. Thermodynamics and dynamics of the inherent states at the glass transition. J. Non-Cryst. Solids 407, 179183 (2015).

56. Grzybowska, K., Grzybowski, A., Pawlus, S., Pionteck, J. & Paluch, M.

Role of entropy in the thermodynamic evolution of the time scale of molecular dynamics near the glass transition. Phys. Rev. E 91, 062305 (2015).

57. Harris, K. R. & Kanakubo, M. Self-diffusion, velocity cross-correlation, distinct diffusion and resistance coefcients of the ionic liquid [BMIM][Tf2N] at high pressure. Phys. Chem. Chem. Phys. 17, 2397723993 (2015).

58. Ingebrigtsen, T. S. & Tanaka, H. Effect of size polydispersity on the nature of Lennard-Jones liquids. J. Phys. Chem. B 119, 1105211062 (2015).

59. Kipnusu, W. K. et al. Connement for more space: a larger free volume and enhanced glassy dynamics of 2-ethyl-1-hexanol in nanopores. J. Phys. Chem. Lett. 6, 37083712 (2015).

60. Schmelzer, J. W. P. & Tropin, T. V. Kinetic criteria of glass-formation, pressure dependence of the glass-transition temperature, and the PrigogineDefay ratio.J. Non-Cryst. Solids 407, 170178 (2015).61. Khrapak, S. A., Klumov, B., Couedel, L. & Thomas, H. M. On the long-waves dispersion in Yukawa systems. Phys. Plasmas 23, 023702 (2016).62. Adrjanowicz, K., Paluch, M. & Pionteck, J. Isochronal superpositionand density scaling of the intermolecular dynamics in glass-forming

liquids with varying hydrogen bonding propensity. RSC Adv. 6, 49370 (2016).63. Powles, J. G., Rickayzen, G. & Heyes, D. M. Temperatures: old, new and middle aged. Mol. Phys. 103, 13611373 (2005).

64. Rosenfeld, Y. New method for equation-of-state calculations: linear combinations of basis potentials. Phys. Rev. A 26, 36333645 (1982).

65. Wilson, H. A. On the velocity of solidication and viscosity of supercooled liquids. Philos. Mag. 50, 238250 (1900).

66. Frenkel, J. Kinetic Theory of Liquids (Clarendon, 1946).67. Auer, S. & Frenkel, D. Prediction of absolute crystal-nucleation rate in hard-sphere colloids. Nature 409, 10211023 (2001).

68. Chui, S. T. Solidication instability of quantum uids. Phys. Rev. B 41, 796798 (1990).

69. Bailey, N. P., Pedersen, U. R., Gnan, N., Schrder, T. B. & Dyre, J. C. Pressureenergy correlations in liquids. II. Analysis and consequences. J. Chem. Phys. 129, 184508 (2008).

70. Bailey, N. P., Bhling, L., Veldhorst, A. A., Schrder, T. B. & Dyre, J. C. Statistical mechanics of Roskilde liquids: congurational adiabats, specic heat contours, and density dependence of the scaling exponent. J. Chem. Phys. 139, 184506 (2013).

Acknowledgements

We are indebted to Karolina Adrjanowicz and Kristine Niss for inspiring the present work via their attempts at interpreting crystallization data for van der Waals liquids in terms of the isomorph theory. This work was supported by the Villum Foundation by the YIP Grant VKR-023455 and by the Danish National Research Foundation by Grant DNRF61.

Author contributions

The project was conceived by U.R.P. and N.P.B. Computations were carried out by U.R.P. and L.C. The theory was devised by J.C.D. based on initial works by N.P.B. and T.B.S. The paper was written by U.R.P. and J.C.D.

Additional information

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Pedersen, U.R. et al. Thermodynamics of freezing and melting. Nat. Commun. 7:12386 doi: 10.1038/ncomms12386 (2016).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

r The Author(s) 2016

NATURE COMMUNICATIONS | 7:12386 | DOI: 10.1038/ncomms12386 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9