Document [original]

This version is available at https://doi.org/10.14279/depositonce-10831

right to use is granted. This document is intended solely for

personal, non-commercial use.

This article may be downloaded for personal use only. Any other use requires prior permission of the

author and AIP Publishing. This article appeared in J. Chem. Phys. 153, 104506 (2020) and may be found

at https://doi.org/10.1063/5.0015371.

Mausbach, P., Fingerhut, R., & Vrabec, J. (2020). Structure and dynamics of the Lennard-Jones fcc-solid

focusing on melting precursors. The Journal of Chemical Physics, 153(10), 104506.

https://doi.org/10.1063/5.0015371

Peter Mausbach, Robin Fingerhut, Jadran Vrabec

Structure and dynamics of the Lennard-

Jones fcc-solid focusing on melting

precursors

Accepted manuscript (Postprint)Journal article |

Structure and dynamics of the Lennard-Jones fcc-solid

focusing on melting precursors

Peter Mausbacha, Robin Fingerhutb, and Jadran Vrabecb,*

aPlant and Process Engineering, Technical University of Cologne, 50678 Cologne, Germany

bThermodynamics and Process Engineering, Technical University Berlin, 10587 Berlin,

Germany

*Corresponding author, E-mail address: vrabec@tu-berlin.de

Abstract

The Lennard-Jones potential is taken as a basis to study structure and dynam-

ics of the face centered cubic (fcc) solid along an isochore from low temperatures

up to the solid/fluid transition. The Z method is applied to estimate the melting

point. Molecular dynamics simulations are used to calculate the pair distribution

function, numbers of nearest neighbors and the translational order parameter, an-

alyzing the weakening of the fcc-symmetry due to emerging premelting effects. A

range of dynamic properties, such as mean-squared displacement, non-Gaussian pa-

rameter, Debye-Waller factor and the vibrational density of states, is considered

for the analysis of the solid state. All of these parameters clearly show that bulk

mobility is activated at about 2/3 of the melting temperature, known as Tammann

temperature. This indicates that vibrational motion of atoms is not maintained

exclusively in the entire stable solid state and that collective atomic motion consti-

tutes a precursor of the melting process.

Keywords: Lennard Jones fcc-solid, solid/fluid transition, melting precursors, Z

method, pair distribution function, numbers of nearest neighbors, translational or-

der parameter, mean-squared displacement, non-Gaussian parameter, Debye-Waller

factor, vibrational density of states, Tammann temperature.

1 Introduction

The melting process of crystals is a long-standing, highly debated issue of condensed

matter physics. Its microscopic mechanism is not yet completely understood and a sat-

isfactory theory is still missing. Normal melting is typically initiated at free surfaces or

defects as a heterogeneous process occuring at the melting temperature Tm. A thermody-

namic criterion for determining Tmis the equality of the Gibbs energy of orthobaric solid

and fluid phases. However, under ideal conditions, melting can progress homogeneously

from the interior of a surface-free, perfect crystal. In such a situation, superheating of the

crystal well above its equilibrium melting temperature is possible. Then, the phase tran-

sition to the fluid state starts from a metastable solid state at the limit of superheating

(LS), characterized by the temperature TLS.

Several stability theories have been developed to interpret the melting mechanism.

Lindemann [1] proposed a theory in which the vibration of atoms reaches a critical fraction

of the nearest neighbor distance. Born [2] suggested that melting occurs when one of

the shear moduli of the system vanishes. Jin et al. [3] were able to connect the ideas

of Lindemann and Born by analyzing clusters of defective atoms, so-called Lindemann

particles, which exhibit strong spatial correlations. Furthermore, Fecht and Johnson [4]

proposed an entropy catastrophe theory that is analogous to Kauzmann’s paradox [5]

for the glass transition. Other models deal with dislocation lines composed of defective

atoms [6–8] that are frequently defined as those which do not have 12 nearest neighbors in

a face-centered cubic (fcc) lattice. In several recent studies [9–12] melting was discussed

as a diffusion-mediated process, exhibiting string-like cooperative atomic motion, which

is a form of dynamic heterogeneity.

Many of the above-mentioned studies on the solid-fluid (S/F) transition focus on the

superheated regime because various properties exhibit a specific anomalous behavior as

they approach the LS. However, a study of K¨oster et al. [13] showed that the onset of

accelerated ”excitation” of many thermodynamic properties starts well below the melt-

ing temperature in the stable solid region and that has frequently been suggested as an

indicator for premelting. This issue is briefly reviewed in the supplementary material,

taking the isochoric heat capacity cVas an example (cf. Figs. S1 and S2). The ther-

modynamics of the considered state region was entirely described in Ref. [13], thus, a

complementary study on the structure and dynamics of the Lennard-Jones (LJ) fcc-solid

is of interest and the aim of the present study. It focuses on aspects of diffusion-mediated

models [9–12] with a special consideration of melting precursors, reflecting anomalous

behavior of thermodynamic response functions.

The investigations are supplemented by estimations of the LS, melting point (MP)

and freezing point (FP). The so-called Z method is applied for this purpose, which has

recently become widespread in the literature [10,14–17]. Assessments of finite size effects

are carried out over a large range of atom numbers, which has not yet been done for the

Z method. In addition, the FP is determined with the Hansen-Verlet criterion.

The paper is organized as follows: First, the molecular simulation method for the

analysis and its associated parameter settings is introduced. A discussion of the Z method

estimating the LS and the MP follows. Next, structural changes of the fcc-solid as it

approaches the S/F transition are examined. Subsequently, the dynamic behavior is

elucidated. Conclusions sum up the findings.

2 Molecular simulation method

The present study is limited to atoms interacting via the LJ potential commonly expressed

uLJ = 4 εσ

r12

−σ

r6,(1)

where σand εare size and energy parameters, while ris the distance between two atoms.

The usual conventions were used for the reduced density (ρ∗=ρ σ3), temperature (T∗=

kBT/ε), energy (u∗=u/ε), distance (r∗=r/σ) and time (t∗=tpε/mσ2), with Boltz-

mann’s constant kB. The reduced form of the other discussed properties is: translational

order parameter (τ∗=τ/σ), mean-squared displacement (h∆r∗2(t)i=h∆r2(t)i/σ2),

Debye-Waller factor (hu∗2i=hu2i/σ2), self-diffusion coefficient (D∗=Dpm/ε σ2) and

vibrational density of states ρ∗

V=ρVpε/mσ2). Pair distribution function (g(r)), non-

Gaussian parameter (α2(t)) and normalized velocity autocorrelation function (ψ(t)) are

dimensionless by definition. All quantities quoted in this work are given in terms of these

reduced quantities so that the asterisk superscript will be omitted in the remainder.

Throughout the paper, molecular dynamics (MD) simulations were performed, inte-

grating Newton’s equations of motion with a fifth-order Gear predictor-corrector scheme

by using the molecular simulation tool ms2 [18–20]. The microcanonical (NVE) ensemble

was employed for the Z method, whereas all other simulation data were sampled in the

canonical (NVT) ensemble. Since long-range effects are important in solids, the atomic

interactions were explicitly evaluated up to a large cutoff radius rc= 9.1 in all NVT simu-

lations. Beyond that distance, analytical long-range corrections were applied. To exclude

vacancy effects, all simulations were initiated from a perfect fcc-crystal. Consequently,

with a cubic simulation volume, the atom number must adhere to N= 4 i3, where iis an

integer. Simulations were performed along the isochore ρ= 1.8 for temperatures ranging

from T= 1 to 26. The isochore ρ= 1.8 was chosen in order to cover a large region with

anomalous behavior (cf. Fig. S2 in the supplementary material).

For each state point, the system was sufficiently equilibrated and then sampled dur-

ing the production period. As a result of the finite size investigations in section 3, an

atom number N= 10 976 was used for most simulations. The specific simulation param-

eters (such as the cutoff radius rc, time step ∆tor simulation length) were adapted to

the requirements of the function to be determined and are compiled in Table S1 of the

supplemental material.

3 Z method and finite size effects

Methods for determining the S/F phase transition can roughly be classified in two groups.

The first group comprises phenomenological procedures based on atomistic simulations

revealing a specific behavior of a first-order phase transition. This includes, e.g., the

evaluation of an equation of state (EOS) that is employed to identify discontinuous jumps

or hysteresis loops as a hallmark of phase transition. Other approaches test the phase

equilibrium conditions with direct simulations of an inhomogeneous system containing

the two coexisting phases separated by an interface. In the second group, methods are

applied that explicitly evaluate the Gibbs energy of the system, which allows for a clear

localization of coexistence between the two phases.

In the present study, a phenomenological method was applied which determines the

S/F transition starting from the solid state. The so-called Z method allows for an es-

timation of the LS and the MP. In order to reach the LS, the system has to evolve on

its own without external interference into the dynamics of the melting process as, e.g.,

by constraining the temperature. Therefore, MD simulations in the NVE ensemble have

to be applied [10, 14]. The idea of the Z method is to equilibrate the system to some

temperature at a given total energy. In a certain total energy range, the crystal resides in

a metastable solid state. That range is delimited by uLS, where the crystal spontaneously

melts at TLS. A total energy slightly above uLS leads to a temperature drop to the melting

temperature Tmbecause the kinetic energy has to supply the internal energy of fusion.

At this point, the total energies of the solid at TLS and the fluid at Tmare assumed to be

equal if the volume/density was conserved [10,14]

usolid(ρ, TLS) = ufluid(ρ, Tm).(2)

Fig. 1 depicts time dependence of total energy uand temperature Tat four state

points, calculated with an atom number N= 10 976. These state points are marked

in the (u, T) representation of Fig. 2 by red arrows 1 to 4. Within the scope of usual

simulation accuracies, the total energy remains constant over the entire simulation period

(note the very magnified scale of u). State point 1 starts with u= 84.00025 and resides

in a metastable solid state with constant temperature. State point 2 with u= 87.00025

is just below uLS, being still in the metastable solid. At u=uLS = 87.25025 (state

point 3) the temperature lingers for a certain time at the LS, then drops to the melting

temperature Tmand then remains constant. At state point 4 with a total energy of

u= 89.00025 > uLS, the temperature drops almost immediately towards a fluid state

temperature.

Collecting (u, T) data obtained from a series of NVE simulations results in a Z shaped

curve, where the maximum temperature corresponds to TLS and that at the local minimum

to Tm. Since studies on finite size effects related to the Z method are rare in the literature,

attention is given to this issue. For this purpose, simulations were performed for atom

numbers N= 256,500,1 372,4 000,10 976,32 000 and 108 000, all initially arranged as a

perfect fcc-crystal.

0 10 20 30

83.9995

84.0000

84.0005

86.9995

87.0000

87.0005

87.2495

87.2500

87.2505

88.9995

89.0000

89.0005

state point 1 (u = 84.00025)

state point 2 (u = 87.00025)

state point 3 (u = 87.25025)

state point 4 (u = 89.00025)

TLS

Figure 1: Instantaneous total energy u(black) and temperature T(blue) as a function

of time tobtained from NVE simulations with an atom number N= 10 976. The left

vertical axis shows the scale for uand the right vertical axis shows the scale for T. Four

different phase points are shown, 1 and 2 belong to the metastable solid state, 3 is at the

LS and 4 is in the fluid state.

80 85 90 95

N = 1372

N = 10976

N = 108000

Solid branch

Fluid branch

Figure 2: Temperature Tas a function of total energy uobtained from NVE simulations

with the Z method. Results for different atom numbers N= 1 372,10 976 and 108 000

along the isochore of ρ= 1.8 are shown. The upper curves belong to the solid branch, the

lower curves to the fluid branch. Red arrows 1 to 4 mark the four state points discussed

in Fig. 1. Vertical lines indicate the temperature drop from the LS to the MP according

to Eq. (2).

103104105

Figure 3: Melting temperature Tmas a function of system size Nalong the isochore

ρ= 1.8. Tmwas obtained from the Z method according to Fig. 2 and Fig. S3 in the

supplemental material. The red line is a power-law fit for system sizes between N= 500

and 108 000.

For clarity, Fig. 2 shows typical results for N= 1 372,10 976 and 108 000, other atom

numbers are discussed in the supplemental material (Fig. S3). Eq. (2) is satisfied along

the vertical lines in Fig. 2. Interestingly, the MP and LS show a rather strong dependence

on the atom number, which becomes even clearer when the melting temperature Tmis

shown as a function of N, cf. Fig. 3. For N≥500, the melting temperature decreases

with rising atom number N, roughly following a power law. Obviously, an approach to

an asymptotic limit of the melting temperature cannot be achieved for atom numbers

analyzed in this study. This corresponds to the behavior of the melting temperature

of silicon examined in a study by Nieves and Sinno [21] for a large range of particle

numbers between N= 2 744 and 238 328. Nieves and Sinno concluded that finite size

effects intrinsically extend to rather large system sizes because critical defect clusters are

not well sampled unless the system size is very large. Therefore, additional superheating

must be provided for smaller system sizes in order to initiate crystal melting. Results

that support these findings have been published by Bai and Li [22], who examined the

critical liquid nucleus radius of a superheated LJ crystal. Finite size investigations for

N= 32 000, 108 000 and 256 000 showed a strong size dependence. Aspects like these

have not yet been properly investigated for the Z method.

In another study, Abramo et al. [23] analyzed the system size dependence of the S/F

phase transition for a LJ system between N= 108 and 4 000, using a phenomenological

method comparable to that employed in Ref. [13]. The authors concluded that with small

atom numbers, specifically N= 108, the best agreement with respect to the location

of the freezing line and melting line can be found. However, this statement cannot be

confirmed on the basis of the present results obtained with the Z method.

Although finite size effects are important for the determination of the S/F transition,

their quantitative influence is rather small for the thermodynamic properties when N≥

1 372. This can be recognized, for example, from the almost congruent temperature

curves T=T(u) in Fig. 2 and Fig. S1 in the supplemental material as well as from

other properties [13]. This is an important result since the present study focuses on

melting precursors, appearing well below the MP. The LS serves as an upper bound of

the metastable solid state and reaching it truly is not of central importance here. However,

the atomistic ensemble should be sufficiently large, at least larger than typical structural

inhomogeneities arising when the MP is approached. Therefore, in the remainder of the

present study, N= 10 976 was chosen.

The present estimation of the MP and LS is supplemented by the determination of the

FP at ρ= 1.8. The Hansen-Verlet criterion [24] was applied, where freezing is indicated

if the magnitude of the first peak of the fluid structure factor reaches a value of 2.85. The

criterion is fulfilled for Tf= 24.775, details are discussed in the supplemental material.

In summary, the results of the present study for the temperatures determining the

S/F transition at ρ= 1.8 are Tm≈20.18, TLS ≈24.2 and Tf≈24.775.

4 Structure

The melting process of a crystal can be thought of as an order-disorder transition in which

the crystal loses its long-range translational order. In this context, the number of first

nearest neighbors (1NN) is a basic order parameter. Since for a perfect fcc-crystal 1NN =

12, any deviation from this number indicates structural changes in the crystal symmetry.

For example, a missing neighboring atom leaves free volume in the lattice structure with

important consequences for the melting behavior. Other defects arise when atoms are

located on interstitial sites, pushing the surrounding atoms out of their equilibrium lattice

positions. The key function for examining 1NN, or any other number of neighbors from

which defective behavior may be observed, is the pair distribution function (PDF).

4.1 Pair distribution function

The PDF g(r) was sampled along the isochore ρ= 1.8 in a temperature range from

T= 1 to 26. Fig. 4 depicts a 3D representation of the PDF up to an atom-atom

distance of r≤1.95. The solid phase is characterized by sharp peaks around the lattice

sites, which are increasingly smeared out as the temperature rises. The magnitude of the

peaks decreases drastically with raising temperature. Moreover, the second and fourth

peak almost disappear within the stable solid phase, weakening the fcc geometry. In a

relatively narrow temperature range from T= 22 to 23, the PDF rapidly smoothens,

indicating the transition from the crystal to the fluid state. Within the fluid state, the

second and fourth peaks of the solid phase PDF disappear completely .

In order to localize the transition temperature from the solid to the fluid phase more

precisely, PDF are shown in Fig. 5 ranging from T= 22 to 22.6. The PDF at T= 22

and 22.2 as well as at T= 22.4 and 22.6 are almost congruent in Fig. 5 so that the S/F

transition was identified to be between T= 22.2 and 22.4.

4.2 Running coordination number and nearest neighbors

The running coordination number

n(r)=4πρ Zr

r02g(r0)dr0,(3)

allows for the determination of the number of neighbors in spherical shells around a ref-

erence atom. For the considered fcc-solid, n(r) is depicted in Fig. S7 in the supplemental

material. Defining the jth spherical shell around the jth peak between the (j−1)-th and

jth minima (rj−1, min, rj, min) of the PDF, the number of next neighbors in this shell is

given by jNN = n(rj, min)−n(rj−1, min), with minimum positions rj, min compiled in Table

1 and visualized in Fig. S8 in the supplemental material. Relying on this definition, the

number of atoms in each shell may change with temperature and density, which allows

Figure 4: 3D representation of the of the first four PDF peaks along the isochore ρ= 1.8

for T= 1 to 26. Upon the approach of the MP, the second and fourth peaks disappear.

The S/F transition region is yellow, the fluid phase is red.

0 2 4 6 8

g(r)

3T =

22.2

22.4

22.6

Figure 5: Pair distribution function g(r) along the isochore ρ= 1.8 for T= 22 to 22.6.

Due to the very small numerical differences between the data for T= 22 and 22.2 as well

as for T= 22.4 and 22.6, the black line is behind the red line and the green line behind

the blue line.

Table 1: Minima separating the peaks of the PDF at T= 1 and ρ= 1.8.

j rj, min

0 0

1 1.1108

2 1.4524

3 1.7283

4 1.9512

for statements about the lattice order. Next neighbors 1NN and 2NN are shown in Fig.

6, next neighbors 3NN and 4NN are shown in Fig. S9 in the supplemental material.

Studying the microscopic mechanism of LJ fcc-crystal melting, G´omez et al. [7] in-

vestigated dislocation lines, suggesting that the appearance of these lines constitutes a

precursor of the melting process. At temperatures above ∼0.8Tm, thermal fluctuations

and distortions should create defective atoms, clustering into dislocation loops when these

atoms become neighbors of each other. Knowing the individual coordination numbers,

it is possible to identify under-coordinated and over-coordinated atoms as point defects.

G´omez et al. defined a defective atom when its number of first nearest neighbors 1NN is

not 12.

The course of 1NN in Fig. 6 evidences a systematic and gradual reduction of the

number of first nearest neighbors with increasing temperature, dropping from 11.7 to

11.2 at the S/F transition. Therefore, a large number of atoms must be individually

under-coordinated, which confirms the findings of G´omez et al. [7]. Furthermore, 2NN

increases from 6 to 7.4 near the S/F transition and jumps to 8.9 in the fluid phase.

Thus, in the solid phase, on average more than one atom migrates into the 2NN shell

due to isochoric heating, while its peak structure is lost. Fig. 6 shows the present data

for MP, LS and FP in comparison to MP/FP coordinates from the literature. As a

general trend, the FP agrees well for a number of studies. The FP coordinates of K¨oster

0 5 10 15 20 25

1NN, 2NN

1NN

2NN

LS (Z method)

MP (Z method)

FP (Köster et al.)

MP (Köster et al.)

FP (Schultz & Kofke)

MP (Schultz & Kofke)

FP (Gottschalk)

FP (Ahmed & Sadus)

FP (Hansen-Verlet criterion)

Figure 6: Temperature dependence of next nearest neighbors 1NN and 2NN along the

isochore ρ= 1.8 for T= 1 to 26. The vertical blue dashed line marks the MP, the blue

solid line the LS calculated with the Z method and the red vertical dashed dotted line

marks the FP according to the Hansen-Verlet criterion (cf. supplemental material). All

data were generated with N= 10 976. Additionally shown is the FP from Gottschalk [25]

and Ahmed and Sadus [26] as well as the FP/MP from K¨oster et al. [13] and Schultz and

Kofke [27].

et al. [13], Gottschalk [25] and the results according to the Hansen-Verlet criterion are

almost identical at T≈25. The FP of Ahmed and Sadus [26] is also close to this value.

The situation is different for the MP, where the associated values vary over a larger

temperature range. Note that the S/F transition drop of next nearest neighbors 1NN

and 2NN occurs before the LS (blue vertical solid line in Fig. 6 with TLS ≈24.2 for

N= 10 976). This indicates that the LS can not be reached with the NVT ensemble and

confirms that the NVE ensemble needs to be applied to estimate the LS [10]. For clarity,

only MP, LS and FP data obtained in this study are shown in the following figures.

Another interesting aspect can be seen in Fig. 7 that depicts the temperature depen-

dence of g(r) for the first four peaks. The arrow at the top left marks the S/F transition

temperature. The rapid disappearance of the second and fourth peaks can clearly be seen.

More interesting is the fact that the 1NN shell is compressed during the isochoric tem-

perature rise (the first peak moves to smaller interatomic distances), while the position of

the other three peaks is fairly independent on temperature. Considering the interstitial,

atom-free interval 1.034 < r < 1.189 between the first and second peak at T= 1, it was

observed that this shell is filled up with atoms when the temperature rises. Fig. 8 shows

the temperature dependence for nint(1p, 2p) = n(1.189) −n(1.034). Near the MP, on

average about two atoms migrate into this interstitial shell, compressing the 1NN shell.

4.3 Translational order parameter

The strong increase of 2NN with a simultaneous loss of the peak structure shows that the

number of nearest neighbors jNN is not a good quantity to assess the translational order

of a system. Errington et al. [28] introduced a metric that quantifies the deviation of an

actual structure from a reference arrangement by

Figure 7: Contours of the temperature dependence of g(r) along the isochore ρ= 1.8 for

T= 1 to 26. The arrow marks the S/F transition.

0 5 10 15 20 25

nint (1p, 2p)

LS (Z method)

MP (Z method)

FP (Hansen-Verlet criterion)

Figure 8: Temperature dependence of the atom number nint within the interstitial spher-

ical shell 1.034 < r < 1.189 between the first and second peak of the PDF. Symbols and

conditions are the same as in Fig. 6.

τ=Zrc

|g(r)−1|dr , (4)

where rccorresponds to a large cutoff radius. A completely uncorrelated system (g(r) = 1)

gives τ= 0, while a crystal with long-range order exhibits large τvalues. The translational

order parameter τcan be decomposed into τ1, τ2, τ3, ... by integrating |g(r)−1|over the

sections 0 < r < r1, min, r1, min < r < r2, min, r2, min < r < r3, min, ... so that τ=τ1+τ2+

τ3+... . Fig. 9 shows the temperature dependence of τ1, ..., τ4and τ=τ1+τ2+τ3+τ4.

As expected, all parameters decrease with raising temperature, indicating a gradual loss

of translational order. Similar to other properties, a strong drop can be observed at the

transition to the fluid state. The order parameters of the second and fourth shell, τ2

and τ4, decrease rapidly even at relatively low temperatures so that they are only weakly

dependent on temperature for T > 10. In contrast, τ1, τ3and τdecrease without this

characteristic.

In summary, restructuring within various shells, starting in the stable solid region far

away from the melting point, was observed. It is caused by activated diffusion as discussed

in the following section.

5 Dynamic properties

Several recent studies examined correlated diffusion processes near the S/F transition, fo-

cusing on the region close to the LS. Delogu [9] showed that melting of a superheated LJ

fcc-crystal originates from defectively coordinated atoms arranged in pseudolinear clus-

ters. The superheated crystal is, in this picture, dynamically heterogeneous, whereby

defective atoms have a greater mobility than normally (12-fold) coordinated atoms. Be-

lonoshko et al. [10] also studied the occurrence of diffusion and defects in a fcc-solid,

concluding that melting is likely to be associated with the formation of elongated struc-

Figure 9: Temperature dependence of translational order parameters τand τ1(a)) as well

as τ2,τ3and τ4(b)). Symbols and conditions are the same as in Fig. 6.

tures of dislocations. Common to these dynamical approaches is that string-like collective

atomic motion occurs where groups of atoms follow each other along one-dimensional

paths. Bai and Li [11] observed rapid self-diffusive motion as a precursor of melting when

the system is kept at a constant temperature below the LS. They suggest that two types of

string-like collective motion drive the homogeneous melting process, i.e. open and closed

strings forming linear or ring (polymeric) structures. Qualitative evidence was provided

that open diffusion loops, rather than closed ones, are related to the homogeneous nu-

cleation of the fluid phase. In another study, Zhang et al. [12] examined similarities

of highly correlated motion in a superheated Ni crystal to those found in glass-forming

liquids. They identified a topological transition at higher temperatures when ring-like

atomic exchanges ”open” to form linear chains of atomic motions, similar to strings in

glass-forming liquids. Ring-like diffusion loops rather preserve the crystalline structure

because of their permutational character. However, the open strings result in a local

symmetry breaking effect, triggering the crystal to melt.

Usually, atoms hardly diffuse in a crystal without vacancies. The main mechanism

that drives diffusion in a perfect crystal is the occurrence of thermally activated jumps

between short-time vacant equilibrium lattice sites [29] arising, e.g., when a neighboring

atom occupies an interstitial site [10]. The behavior of ∆g(r, T) = g(r, T)−g(r, T = 1)

shows how these movements become possible, as discussed in Fig. S10 in the supplemental

material.

Literature studies on diffusion mediated melting focused on regions close to the LS.

Nevertheless, it is also of interest to investigate whether such effects are already initiated

in the stable solid state, as it is the case with premelting effects in terms of structure and

thermodynamic properties [13].

5.1 Mean-squared displacement

A common quantity describing the dynamics of a system is the single atom mean-squared

displacement (MSD) h∆r2(t)i, where ∆r2(t) is the squared distance propagated by a given

atom at time tand the angular brackets denote the ensemble average. Fig. 10.a shows

typical MSD for the solid state. The well-known ballistic regime with h∆r2(t)i ∼ t2can be

recognized at very short times. A crossover from ballistic to vibrational motion follows,

where the MSD exhibits a well-defined plateau, reflecting oscillations of atoms around

their equilibrium lattice positions.

For higher temperatures, an additional regime appears at longer times, where atoms

escape from their equilibrium lattice sites due to diffusion, cf. Fig. 10.b. By analyzing

distributions of the atom displacements, Bai et al. [11] and Zhang et al. [12] showed that

atoms migrate over a distance corresponding to their first neighboring lattice sites. At

longer times, the MSD approaches the Einstein relation of h∆r2(t)i= 6 D t, where Dis

the self-diffusion coefficient.

It was frequently observed that many crystals acquire sufficient energy for their bulk

mobility above an onset temperature TT≈2Tm/3, known as Tammann temperature

[30, 31]. For example, Zhang et al. [32] consider TTas a premelting onset temperature

at which accelerated dynamics can be observed. On the time scale shown in Fig. 10.b,

diffusive motion is initiated at approximately T≈17 ≈0.75 Tm, indicating that pure

vibrational motion is not maintained in entire stable solid state, as discussed in the more

recent literature.

The intermediate vibrational regime disappears completely at a temperature of T=

22.4 and the MSD exhibits fluid behavior for all temperatures above. Note that this is

exactly the temperature where S/F phase transition was observed in the PDF, cf. Fig. 5.

Figure 10: Mean-squared displacement h∆r2(t)ias a function of time talong the isochore

ρ= 1.8 for T= 1 to 26. Typical solid behavior can be observed up to a temperature of

T≈16, consisting of a ballistic regime at very short times and a vibrational regime with

a well-defined plateau for longer times (a)). An additional diffusive regime appears for

T≥17 at longer times (b)).

5.2 Non-Gaussian parameter and Debye-Waller factor

Dynamic heterogeneity caused by mobile and immobile atoms is commonly quantified by

the non-Gaussian parameter (NGP)

α2(t) = 3h∆r4(t)i

5h∆r2(t)i2−1.(5)

It is known that the NGP mainly receives contributions from those atoms that move

further than the Gaussian distribution of particle displacements. Consequently, α2(t) has

been suggested as an indicator for dynamical heterogeneity. Investigations based on the

NGP should therefore allow for conclusions about mobile atoms, diffusing through a crys-

tal consisting of immobile atoms. In case of ballistic and diffusive motion, displacements

are known to follow a Gaussian distribution and consequently the NGP is zero in these

regimes. This behavior can be recognized in Fig. 11.a for the fluid phase in the temper-

ature range from T= 23 to 26. At the crossover between the two regimes, α2develops

a pronounced maximum. For the stable solid phase in a temperature range from T= 1

to 16, the behavior is different. The NGP starts at zero in the ballistic regime and then

peaks at a time t= 0.1, corresponding to the appearance of vibrational motion, followed

by an approach to an almost constant plateau with a relatively small magnitude. At a

temperature of T= 17, the NGP rises again at longer times, comparable to the behavior

of the MSD at this temperature.

Of greater interest is the behavior of α2in the vicinity of the S/F phase transition,

cf. Fig. 11.b. The NGP is depicted there for temperatures from T= 17 to 26. A

dramatic increase of α2can be observed when the system is heated from T= 17 to 19,

which signals the appearance of strongly heterogeneous dynamics. The maximum of α2

at T= 19 corresponds to a characteristic time tmax when displacements are mostly non-

Figure 11: Non-Gaussian parameter α2(t) as a function of time talong the isochore

ρ= 1.8 for T= 1 to 26 (a)). In b), the occurrence of the diffusive regime is characterized

by pronounced peaks of α2(t) between T= 19 and 22.2 before the α2(t) peak suddenly

recedes at T= 23, indicating the S/F transition marked by the arrow.

Figure 12: Mean-squared displacement h∆r2(t)i(black) and non-Gaussian parameter

α2(t) (red) at T= 18.6 and ρ= 1.768. The scale for h∆r2(t)iis given on the left axis,

the scale for α2(t) on the right axis.

Gaussian. Furthermore, the magnitude of α2(tmax) and tmax decreases from T= 19 to

22.2. This is an indication that the transient vibrational regime becomes less pronounced

as the atoms need less and less time to enter the diffusive regime.

The interplay between the MSD and α2(t) is compared in Fig. 12 for the state point

T= 18.6 and ρ= 1.768, which is close to the MP. The NGP peaks strongly at the

transition from the vibrational to the diffusive regime, appearing roughly at tmax in the

MSD. Consequently, the position of the peak tmax can be considered as an activation time

of diffusion, where a single atom escapes from its equilibrium lattice position. Zhang et

al. [12] argued that tmax has the significance of a characteristic diffusive relaxation time.

From the behavior of h∆r2(t)i= 6 D t for long times, the self-diffusion coefficient D

1 / tmax

0.005 0.010 0.015 0.020 0.025

D / 106T

Figure 13: Self-diffusion coefficient Ddivided by temperature Tas a function of inverse

activation time of diffusion tmax.

1 / T

10-1 100

<u2>

10-3

10-2

10-1

TTTcv,min

Figure 14: Debye-Waller factor hu2ias a function of inverse temperature 1/T. The red

line represents a power-law fit for T < TT. The Tammann temperature TT≈13.5 is

marked by an arrow. TcV,min ≈12 marks the crossover temperature of the cV,min line at

ρ= 1.8.

can be calculated at least for high temperatures close to the LS. Fig. 13 shows that D/T

scales linearly with inverse tmax.

Another interesting quantity that can be estimated directly from the MSD is the

Debye-Waller factor (DWF) hu2i, measuring oscillations of atoms around their equilibrium

lattice positions. Starr et al. [33] and Larini et al. [34] defined the DWF as the MSD after

the crossover from the ballistic to the vibrational regime. Following this definition, the

DWF was calculated and is shown in Fig. 14 as a function of inverse temperature in

a double logarithmic plot. Raising the temperature causes an increase of the DWF,

which can be described by a power-law fit. For temperatures T > TT, an accelerated

increase of the DWF can be observed, deviating from the power-law. This is an interesting

result, since the behavior of the DWF reflects different dynamic regimes. When the

Tammann temperature TTis exceeded, increasing oscillation amplitudes grow faster than

the power-law, allowing the atoms to diffuse. Obversely, the power-law behavior of the

DWF indicates purely vibrational motion. In addition, a comparison of the temperature

dependence of the DWF and that of the isochoric heat capacity cVis of interest. As

discussed in the supplemental material (cf. Fig. S2), the line of the minimum cV,min bounds

a phase region with anomalous thermodynamic behavior, characterizing premelting. The

Tammann temperature of TT≈13.5 roughly corresponds to the intersection of the cV,min

curve with the isochore ρ= 1.8 at T≈12 (cf. Fig. 14). This indicates a close connection

between the occurrence of anomalous thermodynamic behavior, reported in Ref. [13], and

the occurrence of bulk mobility in the stable fcc-solid.

5.3 Vibrational density of states

A quantity describing the distribution of vibrational frequencies of an atomic system is

given by the vibrational density of states (VDOS). The VDOS ρV(ω) can be calculated

from the normalized velocity autocorrelation function (VACF) ψ(t) = hv(t)v(0)i/hv(0)2i

using the Fourier transformation

ρV(ω) = 2 Z∞

ψ(t) cos(ω t)dt , (6)

with the frequency ω. Near the S/F transition, a typical behavior of both phases should

clearly emerge in ψ(t) and ρV(ω). VACF and VDOS were calculated for temperatures in

a range from T= 19 to 26. The VACF was sampled over a span of 15 000 time steps,

where its decay is sufficient to execute the integral in Eq. (6).

Fig. 15.a depicts the normalized VACF. As expected, backscattering of atoms is much

Figure 15: a) Normalized velocity autocorrelation function and b) vibrational density of

states along the isochore ρ= 1.8 for T= 19 to 26.

more pronounced in the solid state and the transition to the fluid state can be clearly seen

at a temperature of T= 24. A corresponding behavior was observed for the VDOS as

shown in Fig. 15.b. The position of the main solid state peak shifts slightly towards lower

frequencies upon heating, an effect known as softening a solid. The magnitude of this

peak decreases with rising temperature. A remarkable jump can be observed at T= 23,

widening the vibrational peak of the fluid state. The difference between the solid and

fluid phase is most evident at ω= 0, where ρV(0) is finite for a fluid, while ρV(0) = 0 for

a solid. Since ρV(0) is related to diffusivity, a finite value indicates incipient diffusion.

Zhang et al. [12] examined the reduced VDOS ρV(ω)/ω2[35] of a superheated Ni

crystal and found a boson peak at low ω, whose occurrence was interpreted in terms of

collective atomic rearrangement motions. Regardless whether the low-frequency behavior

appears in the sense of a boson-like peak, an investigation of the reduced VDOS of the

superheated LJ fcc-crystal is of general interest, as discussed in the supplemental material.

6 Conclusions

Homogeneous melting is often accompanied by an anomalous behavior of thermodynamic

properties [13], frequently suggested as an indicator for premelting. A question arising in

this context is: To what extent can corresponding effects be found in the structure and

dynamics of a crystal? In the present study, this issue was analyzed for the LJ fcc-solid

with molecular dynamics simulations. In order to specify the region of interest, the limit

of superheating as an upper bound of the metastable solid state and the melting point

were estimated with the Z method. If melting precursors become apparent in a solid,

a weakening of the crystalline order has to be expected. Therefore, translational order

was investigated by means of the number of nearest neighbors and a translational order

metric based on the PDF. Strong restructuring effects were found especially in the first

and second shell of the PDF, when approaching the S/F transition. The first peak of the

PDF moves to smaller interatomic distances and the number of first nearest neighbors

is reduced. The second shell of the PDF already loses its translation order at about

T≈0.5Tm, resulting in a broad distribution of atoms in this shell. Both effects favor the

formation of thermal vacancies [29], which in turn allow for an activation of self-diffusion

in the crystal. Investigations of the mean-squared displacement showed that diffusion

processes are activated at the Tammann temperature TT≈2Tm/3. An analysis of the

non-Gaussian parameter revealed a remarkable peak when approaching Tm. This indicates

strong heterogeneous dynamics, which is generally attributed to the appearance of string-

like collective atomic motion within the crystal. In addition, a changing temperature

dependence of the Debye-Waller factor at TTsupports the findings of increasing mobility

in the fcc-solid. Furthermore, a pronounced peak was detected in the vibrational density of

states at low frequencies. In summary, melting precursors are present in the structure and

dynamics of the LJ fcc-solid and thus reflect the anomalous behavior of thermodynamic

properties in a large portion of the approach to the melting point.

Acknowledgments

Computational support was given by the High Performance Computing Center Stuttgart

(HLRS) under the grant MMHBF2. Furthermore, we gratefully acknowledge the Pader-

born Center for Parallel Computing (PC2) for the generous allocation of computer time

on the OCuLUS and Noctua cluster. We thank Gabriela Guevara-Carrion for her valuable

support.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding

author upon reasonable request.

Supplemental material

The supplemental material contains a brief review of the thermodynamic results obtained

by K¨oster et al. [13] as well as additional information and tests to complement the main

body of this work. The FP was determined on the basis of the Hansen-Verlet criterion

and was additionally checked by the Ravech´e-Mountain-Streett criterion. Furthermore,

the VDOS at low frequencies is discussed.

References

[1] F. A. Lindemann, ”¨

Uber die Berechnung molekularer Eigenfrequenzen,” Phys. Z. 11,

609 (1919).

[2] M. Born, ”Thermodynamics of Crystals and Melting,” J. Chem. Phys. 7, 591 (1939).

[3] Z. H. Jin, P. Gumbsch, K. Lu, and E. Ma, ”Melting mechanisms at the limit of

superheating,” Phys. Rev. Lett. 87, 055703 (2001).

[4] H. J. Fecht and W. L. Johnson, ”Entropy and enthalpy catastrophe as a stability limit

for crystalline material,” Nature 334, 50 (1988).

[5] W. Kauzmann, ”The nature of the glassy state and the behavior of liquids at low

temperatures,” Chem. Rev. 43, 219 (1948).

[6] L. Burakovsky, D. L. Preston, and R. Silbar, ”Melting as a dislocation-mediated phase

transition,” Phys. Rev. B 61, 15011 (2000).

[7] L. Gomez, A. Dobry, Ch. Geuting, H. T. Diep, and L. Burakovsky, ”Dislocation lines as

the precursor of the melting of crystalline solids observed in Monte Carlo simulations,”

Phys. Rev. Lett. 90, 095701 (2003).

[8] F. Delogu, ”Cooperative dynamics and self-diffusion in superheated crystals,” J. Phys.

Chem. B 109, 15291 (2005).

[9] F. Delogu, ”Cooperative atomic displacements and melting at the limit of superheat-

ing,” J. Phys. Chem. B 110, 3281 (2006).

[10] A. B. Belonoshko, S. Davis, N. V. Skorodumova, P. H. Lundow, A. Rosengren, and

B. Johansson, ”Properties of the fcc Lennard-Jones crystal model at the limit of

superheating,” Phys. Rev. B 76, 064121 (2007).

[11] X.-M. Bai and M. Li , ”Ring-diffusion mediated homogeneous melting in the super-

heating regime,” Phys. Rev. B 77, 134109 (2008).

[12] H. Zhang, M. Khalkhali, Q. Liu, and J. F. Douglas, ”String-like cooperative motion

in homogeneous melting,” J. Chem. Phys. 138, 12A538 (2013).

[13] A. K¨oster, P. Mausbach, and J. Vrabec, ”Premelting, solid-fluid equilibria, and ther-

modynamic properties in the high density region based on the Lennard-Jones poten-

tial,” J. Chem. Phys. 147, 144502 (2017).

[14] A. B. Belonoshko, N. V. Skorodumova, A. Rosengren, and B. Johansson, ”Melting

and critical superheating,” Phys. Rev. B 73, 012201 (2006).

[15] F. Gonz´alez-Cataldo, S. Davis, and G. Guti´errez, ”Z method calculations to deter-

mine the melting curve of silica at high pressures,” J. Phys.: Conf. Ser. 720 012032

(2016).

[16] K. Yin, X. Lu, H. Zhou, and Y. Sun, ”Thermodynamic stability limit of the crystalline

state from the Gibbs perspective,” Phys. Rev. B 98, 144113 (2018).

[17] V. Olgu´ın-Arias, S. Davis, and G. Guti´errez, ”Extended correlations in the critical

superheated solid,” J. Chem. Phys. 151, 064507 (2019).

[18] S. Deublein, B. Eckl, J. Stoll, S. V. Lishchuk, G. Guevara-Carrion, C. W. Glass,

T. Merker, M. Bernreuther, H. Hasse, and J. Vrabec, ”ms2: A molecular simulation

tool for thermodynamic properties,” Comp. Phys. Commun. 182, 2350-2367, (2011).

[19] C. W. Glass, S. Reiser, G. Rutkai, S. Deublein, A. K¨oster, G. Guevara-Carrion,

A. Wafai, M. Horsch, M. Bernreuther, T. Windmann, H. Hasse, and J. Vrabec,

”ms2: A molecular simulation tool for thermodynamic properties, new version re-

lease,” Comp. Phys. Commun. 185, 3302-3306, (2014).

[20] G. Rutkai, A. K¨oster, G. Guevara-Carriona, T. Janzen, M. Schappals, C.W. Glass,

M. Bernreuther, A. Wafai, S. Stephan, M. Kohns, S. Reiser, S. Deublein, M. Horsch,

H. Hasse, and J. Vrabec, ”ms2: A molecular simulation tool for thermodynamic prop-

erties, release 3.0,” Comp. Phys. Commun. 221 , 343-351 (2017).

[21] A. M. Nieves and T. Sinnoa, ”An enthalpy landscape view of homogeneous melting

in crystals,” J. Chem. Phys. 135, 074504 (2011).

[22] X.-M. Bai and M. Li, ”Differences between solid superheating and liquid supercool-

ing,” J. Chem. Phys. 123, 151102 (2005).

[23] M. C. Abramo, C. Caccamo, D. Costa, P. V. Giaquinta, G. Malescio, G. Muna`o,

and S. Prestipino, ”On the determination of phase boundaries via thermodynamic

integration across coexistence regions,” J. Chem. Phys. 142, 214502 (2015).

[24] J.-P. Hansen and L. Verlet, ”Phase Transitions of the Lennard-Jones System,” Phys.

Rev. 184, 151 (1969).

[25] M. Gottschalk, ”An EOS for the Lennard-Jones fluid: A virial expansion approach,”

AIP Adv. 9, 125206 (2019).

[26] A. Ahmed and R. J. Sadus, ”Solid-liquid equilibria and triple points of n-6 Lennard-

Jones fluids,” J. Chem. Phys. 131, 174504 (2009); Erratum 133, 229902 (2010).

[27] A. J. Schultz and D. A. Kofke, ”Comprehensive high-precision high-accuracy equa-

tion of state and coexistence properties for classical Lennard-Jones crystals and low-

temperature fluid phases,” J. Chem. Phys. 149, 204508 (2018).

[28] J. R. Errington, P. G. Debenedetti, and S. Torquato, ”Quantification of order in the

Lennard-Jones system,” J. Chem. Phys. 118, 2256 (2003).

[29] M. J. Pozo, S. Davis, and J. Peralta, ”Statistical distribution of thermal vacancies

close to the melting point,” Physica B 457, 310 (2015).

[30] G. Tammann, ”Die Temperatur des Beginns innerer Diffusion in Kristallen,” Z.

Anorg. Allg. Chem. 157, 321 (1926).

[31] G. Tammann ,”Lehrbuch der Metallkunde. Chemie und Physik der Metalle und ihrer

Legierungen,” 4. erweiterte Auflage, Leopold Voß, Leipzig (1932).

[32] H. Zhang, X. Wang, A. Chremos, and J. F. Douglas ,”Superionic UO2: A model

anharmonic crystalline material,” J. Chem. Phys. 150, 174506 (2019).

[33] F. W. Starr, S. Sastry, J. F. Douglas, and S. C. Glotzer, ”What do we learn from

the local geometry of glass-forming liquids?,” Phys. Rev. Lett. 89, 125501 (2002).

[34] L. Larini, A. Ottochian, C. De Michele, and D. Leporini, ”Universal scaling between

structural relaxation and vibrational dynamics in glass-forming liquids and polymers,”

Nat. Phys. 4, 42 (2008).

[35] T. S. Grigera, V. Martin-Mayor, G. Parisi, and P. Verrocchio, ”Phonon interpretation

of the ’boson peak’ in supercooled liquids,” Nature 422, 289 (2003).

Supplemental Material to:

Structure and dynamics of the Lennard-Jones fcc-solid

focusing on melting precursors

Peter Mausbacha, Robin Fingerhutb, and Jadran Vrabecb,*

aPlant and Process Engineering, Technical University of Cologne, 50678 Cologne, Germany

bThermodynamics and Process Engineering, Technical University Berlin, 10587 Berlin,

Germany

*Corresponding author, E-mail address: vrabec@tu-berlin.de

This supplementary material contains a brief review of the current status of thermo-

dynamic results obtained from K¨oster et al. [1] as well as additional information and tests

to complement the main text.

1 Review on thermodynamic anomalies near melting

K¨oster et al. [1] determined the entire set of thermodynamic properties of the Lennard-

Jones fcc-solid over a large phase region. Close to the solid-fluid (S/F) phase transition,

several properties showed an accelerated anomalous behavior. In this section, this issue

is reviewed exemplarily for the isochoric heat capacity cV.

Fig. 1.a depicts density dependent cVdata from simulations of K¨oster et al. [1] along

three isotherms T= 1.3, 6 and 22 in comparison with results from an EOS developed by

Schultz and Kofke [2], supplemented by an EOS taking vacancy effects into account. In

general, cVvalues of both studies [1,2] show a good agreement in the solid region. This

is a challenging comparison since the isochoric heat capacity is a second order derivative

of the Helmholtz energy so that the consistency of the cVdata of K¨oster et al. [1] can be

deduced.

For the temperature T= 1.3, the freezing point (FP) and melting point (MP) coor-

dinates of Refs. [1,2] are almost identical. Here, the isochoric heat capacity cVexhibits

a pronounced minimum when approaching the MP, with a position roughly at the MP.

Entering the metastable two-phase coexistence region by reducing the density causes an

anomalous increase of cV, followed by a sudden drop close to the FP which signals that

the limit of superheating (LS) is reached. The simulation data of the fluid phase together

with the metastable data exhibit a λ-like shape at the S/F phase transition as discussed

in Ref. [1]. Increasing the temperature to T= 6 and 22 leads to quantitatively somewhat

diverging coexistence lines of Refs. [1] and [2]. However, the minimum of cVshifts to the

stable solid region, regardless of differing MP locations in Refs. [1] and [2].

A similar behavior can be observed for the temperature dependence of the isochoric

heat capacity along the isochore ρ= 1.3 as shown in Fig. 1.b. Therein, the temperature

was scaled with the melting temperature Tmof Schultz and Kofke at ρ= 1.3 to allow for

a direct comparison. Note that the formula for the ML in Ref. [2] contains a misprint.

Figure S1: Isochoric heat capacity cV, a) of the fluid and fcc solid phase along three

isotherms T= 1.3, 6 and 22 and b) in the solid phase along the isochore ρ= 1.3. Squares

present simulation data of K¨oster et al. [1] in the solid, circles in the fluid phase. cV

values obtained from the EOS of Schultz and Kofke [2] are shown as black solid lines,

values from an EOS with additional consideration of vacancies are shown as black dashed

lines. Vertical lines indicate the FP and the MP: blue - K¨oster et al. [1], red - Schultz

and Kofke [2]. Vertical dashed lines indicate the MP, vertical solid lines indicate the FP.

The dotted line in a) indicates the Dulong-Petit value cV= 3. Additionally shown in b)

is cVfor small Taccording to the LJD theory of Lustig [3]. All temperatures in b) are

reduced by the melting temperature Tmof Schultz and Kofke at ρ= 1.3.

The pre-factor for ρfcc

melt, Eq. (2), Table I, should be β−1/4instead of β−1/2[4]. Again,

both studies show a good agreement. An accelerated increase of cVcan be observed

within the two-phase coexistence region, however, the onset of the anomalous increase of

cVstarts below the melting temperature at T≈0.8Tm. Such a behavior has frequently

been suggested as an indicator for premelting.

The minimum cV,min opens a possibility to narrow down the phase region with premelt-

ing effects. In Fig. S2, the course of the minimum cV,min at constant density estimated

from the simulation data of K¨oster et al. [1] is shown. A remarkable premelting zone is

visible within the stable solid region. Fig. S2 additionally shows the courses of the FL of

Gottschalk [5] as well as Ahmed and Sadus [6] which are almost identical with the FL of

K¨oster et al. [1]. The MP and LS (both with the Z method, cf. section 3 main-text) as

well as the FP (with the Hansen-Verlet criterion, cf. section 4) are depicted by symbols

in Fig. S2.

Figure S2: Density ρas a function of temperature Talong the freezing (FL) and melting

lines (ML) of K¨oster et al. [1] and Schultz and Kofke [2]. The FL by Gottschalk [5]

and Ahmed and Sadus [6] is shown as well. Furthermore, the dash-dotted line indicates

the minimum cV,min at constant density estimated from the simulation data of K¨oster et

al. [1]. The MP and LS, calculated with the Z method, and the FP, calculated with the

Hansen-Verlet criterion for N= 10 976, are shown by symbols.

2 Molecular simulation parameters

Table S1 compiles the simulation parameters used in this study.

Table S1: Simulation parameters for the Z method and functions sampled in this study.

Function length, equilibration and production period are given in time steps ∆t.

Method/ Ensemble N rc∆tFunction Equilibration Production

Function length steps steps

Z NVE 256 2.6 0.000225 106

500 3.2 106

1 372 4.5 3.2×105

4 000 6.0 5 ×105

10 976 6.5 3.2×105

32 000 6.5 5 ×105

108 000 6.5 106

PDF NVT 10 976 9.1 0.001 2 ×105107

MSD NVT 10 976 9.1 0.001 106- 1072×1051 - 2 ×107

NGP

VACF NVT 10 976 9.1 0.0005 15 000 5 ×1058.5×105

3 Finite size effects of the Z method

Fig. S3 depicts (u, T) curves obtained from the Z method for different atom numbers N.

86 88 90 92

N = 256

N = 500

N = 4000

N = 32000

Solid branch

Fluid branch

Figure S3: Temperature Tas a function of total energy uobtained from microcanonical

simulations with the Z method. Results for different atom numbers N= 256,500,4 000

and 32 000 are shown along the isochore ρ= 1.8. Vertical lines indicate the temperature

drop from the LS to the MP.

4 Phenomenological freezing criteria

Several approximate approaches based on structural information have been proposed to

locate the freezing point (FP) coordinates. Perhaps the most successful phenomenological

criterion for determining the freezing transition was introduced by Hansen and Verlet [7]

on the basis of the static structure factor S(q) of a uniform system. It is defined by the

Fourier transformation

S(q) = 1 + 4 πρ Z∞

(g(r)−1) sin(qr)

qr r2dr , (1)

where qis the modulus of the wave vector q. Hansen and Verlet noticed that the mag-

nitude of the first peak of the fluid structure factor S(q) is nearly 2.85 at the FP, which

seems to be a universal feature that has been verified for many systems. Fig. S4 shows

the amplitude of the first peak of S(q) in the temperature range from T= 24 to 26. The

maximum of the first peak exceeds the value 2.85 roughly at T= 24.75. Fig. S5 depicts

the temperature dependence of the intersection with 2.85 at T= 24.775, which is close

to the freezing temperature Tf= 24.972 according to K¨oster et al. [1].

Another approximate rule for locating freezing is the Ravech´e-Mountain-Streett crite-

rion [8], where the Ravech´e parameter of R=gmin/gmax = 0.2±0.02 indicates freezing.

Ris the ratio of the first non-zero minimum magnitude of g(r) to the the first maximum

magnitude of g(r). The temperature dependence of the inverted Ravech´e parameter R−1

is shown in Fig. S6, exhibiting the value R−1= 5 almost exactly on the FP according to

the Hansen-Verlet criterion.

8.45 8.50 8.55 8.60 8.65 8.70 8.75

S(q)1p

2.65

2.70

2.75

2.80

2.85

2.90

T =

24.0

24.25

24.5

24.75

25.0

25.25

25.5

25.75

26.0

Hansen-Verlet

criterion

Figure S4: Temperature dependence of the magnitude of the first peak of S(q) along the

isochore ρ= 1.8 for T= 24 to 26. Green curves represent the fluid region, black curves

the two-phase coexistence region. The red, dash-dotted line marks the Hansen-Verlet

criterion.

24.0 24.5 25.0 25.5 26.0

max(S(q)1p)

2.80

2.82

2.84

2.86

2.88

Hansen-Verlet criterion

FP (Köster et al.)

FP (Gottschalk)

FP (Ahmed & Sadus)

Figure S5: Temperature dependence of the magnitude of the first peak of S(q) depicted

by squares. The vertical lines represent the FP according to K¨oster et al. [1] (pink),

Gottschalk [5] (black) as well as Ahmed and Sadus [6] (green). The red, dash-dotted line

marks the Hansen-Verlet criterion.

14 16 18 20 22 24 26

R -1 = gmax / gmin

LS (Z method)

MP (Z method)

FP (Hansen-Verlet criterion)

Raveché criterion:

4.55 < R-1 < 5.56

Figure S6: Temperature dependence of the inverse Ravech´e parameter R−1. The horizon-

tal dash-dotted line corresponds to R= 0.2, dashed lines mark upper and lower bounds

(all in green). The vertical dashed line marks the MP, the solid line the LS calculated

with the Z method (both in blue). The vertical dashed dotted line (in red) marks the FP

obtained from the Hansen-Verlet criterion. All data were generated with N= 10 976.

5 Running coordination number and nearest neigh-

bors

The running coordination number n(r) is shown in Fig. S7 for selected temperatures be-

tween T= 1 and 24. The pronounced plateaus correspond to the minima (r1, min , r2, min , ...)

separating the peaks of the PDF shown in Fig. S8. The numbers on the right in Fig. S7

indicate the cumulated number of next nearest neighbors for a perfect fcc-crystal, i.e.

1NN = 12, 2NN = 6, 3NN = 24, ... .

Fig. S9 depicts the temperature dependence of 3NN and 4NN, where the numbers 3NN

= 24 and 4NN = 12 of a perfect fcc-crystal change quickly with increasing temperature.

Jumps in the course of 3NN and 4NN indicate the solid/fluid transition.

0.5 1.0 1.5 2.0 2.5

n(r)

100

150 T =

134

140

Figure S7: Running coordination number n(r) along the isochore ρ= 1.8 for T= 1 to 24.

0 1 2 3

g(r)

T =

Figure S8: Pair distribution function g(r) along the isochore ρ= 1.8 for selected tempera-

tures between T= 1 and 26, depicted up to an atom-atom distance of r≤3.45. Bold lines

on the lower horizontal axis mark the minima (r1, min , r2, min , ...) separating the peaks of

the PDF.

0 5 10 15 20 25

3NN, 4NN

3NN

4NN

LS (Z method)

MP (Z method)

FP (Hansen-Verlet criterion)

Figure S9: Temperature dependence of next nearest neighbors 3NN and 4NN along the

isochore ρ= 1.8 for T= 1 to 26. Symbols and conditions are the same as in Fig. S6.

6 Diffusion mechanism in a perfect crystal

As described in the accompanied paper, the main mechanism for the occurrence of dif-

fusion in a perfect crystal are thermally activated jumps among short-time vacant equi-

librium lattice sites [9]. A close inspection of the temperature dependence of ∆g(r, T) =

g(r, T)−g(r, T0) may explain this. ∆g(r, T) with T0= 1 depicted in Fig. S10 shows that

the interstitial space is occupied upon heating, although the long-range order is well visible

as long as the sample remains in the solid state. This creates vacant lattice sites, at least

in the short term [9]. The arrows mark rising temperature, with ∆g(r, T) increasingly

resembling the curves of the fluid phase.

0.5 1.0 1.5 2.0 2.5

∆g(r)

-8

-4

2p 4p

Figure S10: Temperature dependence of ∆g(r, T) = g(r, T)−g(r, T = 1). 1p, 2p, ...

mark the positions of the first, second, ... peaks of the PDF, arrows mark increasing

temperature. Black curves represent the solid, red curves the fluid phase.

7 Low frequency behavior of the vibrational density

of states

Zhang et al. [10] discussed the potential occurrence of a boson peak for a superheated

Ni crystal at low frequencies. In their analysis, they normalized the VDOS with the

Debye vibrational density of states (∼ω2), which is usually used for investigations of

glass-forming liquids [11]. Regardless whether such a peak can really be interpreted

as a boson peak, the low frequency behavior of the VDOS is also of interest for the

superheated LJ crystal. Therefore, a procedure similar to that of Ref. [10] was applied

in the present study. The reduced VDOS ρV(ω)/ω2is depicted in Fig. S11 that allows a

clear identification of a pronounced peak appearing between ω= 4.25 and 5. This peak is

accompanied by smaller oscillations at lower frequencies and larger oscillations at higher

frequencies, which, however, are increasingly smeared out as the temperature rises. The

peak position shifts to smaller frequencies upon heating. The present results for the LJ

fcc-solid correspond to the findings of Zhang et al. [10]. Nevertheless, for an unequivocal

postulation of the existence of a boson peak, further examinations are definitely necessary,

since a rigorous physical theory describing such a peak is missing.

Figure S11: Reduced vibrational density of states ρV(ω)/ω2along the isochore ρ= 1.8

for T= 19 to 23.

References

[1] A. K¨oster, P. Mausbach, and J. Vrabec, ”Premelting, solid-fluid equilibria, and ther-

modynamic properties in the high density region based on the Lennard-Jones poten-

tial,” J. Chem. Phys. 147, 144502 (2017).

[2] A. J. Schultz and D. A. Kofke, ”Comprehensive high-precision high-accuracy equa-

tion of state and coexistence properties for classical Lennard-Jones crystals and low-

temperature fluid phases,” J. Chem. Phys. 149, 204508 (2018).

[3] R. Lustig, ”On the Lennard-Jones and Devonshire theory for solid state thermody-

namics,” Mol. Phys. 115, 1362 (2017).

[4] A. J. Schultz and D. A. Kofke, Erratum to ”Comprehensive high-precision high-

accuracy equation of state and coexistence properties for classical Lennard-Jones

crystals and low-temperature fluid phases,” J. Chem. Phys. 149, 204508 (2018),

submitted (2020).

[5] M. Gottschalk, ”An EOS for the Lennard-Jones fluid: A virial expansion approach,”

AIP Adv. 9, 125206 (2019).

[6] A. Ahmed and R. J. Sadus, ”Solid-liquid equilibria and triple points of n-6 Lennard-

Jones fluids,” J. Chem. Phys. 131, 174504 (2009); Erratum 133, 229902 (2010).

[7] J.-P. Hansen and L. Verlet, ”Phase Transitions of the Lennard-Jones System,” Phys.

Rev. 184, 151 (1969).

[8] H. J. Ravech´e, R. D. Mountain, and W. B. Streett, ”Freezing and melting properties

of the Lennard-Jones system,” J. Chem. Phys. 61, 1970 (1974).

[9] M. J. Pozo, S. Davis, and J. Peralta, ”Statistical distribution of thermal vacancies

close to the melting point,” Physica B 457, 310 (2015).

[10] H. Zhang, M. Khalkhali, Q. Liu, and J. F. Douglas, ”String-like cooperative motion

in homogeneous melting,” J. Chem. Phys. 138, 12A538 (2013).

[11] T. S. Grigera, V. Martin-Mayor, G. Parisi, and P. Verrocchio, ”Phonon interpretation

of the ’boson peak’ in supercooled liquids,” Nature 422, 289 (2003).