Extended correlations in the critical superheated solid

TL;DR

This paper demonstrates the existence of long-range energy correlations in a superheated Lennard-Jones crystal, revealing collective atomic behavior that precedes melting in metastable states.

Contribution

It uncovers long-range potential energy correlations in a superheated solid, highlighting the role of vacancy-interstitial pairs in inducing effective long-range forces.

Findings

Long-range energy correlations exist in superheated crystals.

Vacancy-interstitial pairs promote cooperative atomic motion.

Effective long-range forces extend beyond fourth neighbors.

Abstract

Metastable states in first-order phase transitions reveal interesting behavior about a wide range of systems in statistical mechanics, including spin systems, cellular automata and condensed matter systems. These metastable states are often observed in a microcanonical setting, where they manifest long-range correlations due to collective effects. In this work we show the existence of long-range potential energy correlations between atoms in a microcanonical superheated Lennard-Jones crystal prior to homogeneous melting. Our results suggest that the cooperative motion made possible by the presence of vacancy-interstitial pairs above the melting temperature induces effective long-range interatomic forces even beyond the fourth neighboring layer.

Equations

$$E=\frac{3}{2}N{k}_B{T}_0+{\Phi }_0,$$

$$\Phi ({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)$$

$$=\sum _i\left[\frac{1}{2}\sum _{j\ne i}\phi (|{\mathbf{r}}_i-{\mathbf{r}}_j|)\right],$$

$$\varphi(r)=4\epsilon\left[\Big{(}\frac{\sigma}{r}\Big{)}^{12}-\Big{(}\frac{\sigma}{r}\Big{)}^{6}\right].$$

$$P({\varphi }_i|E)=\frac{1}{Z_1(\beta )}\exp (-\beta {\varphi }_i){\Omega }_1({\varphi }_i),$$

$$P(\varphi |E)=P(\varphi |k,\theta )=\frac{\exp (-\varphi /\theta ){\varphi }^{k-1}}{\Gamma (k){\theta }^k},$$

$$P({\varphi }_1|E)=\frac{1}{Z(\beta )}\exp (-\beta {\varphi }_1){\omega }_1({\varphi }_1;E)$$

$$I=\int dxdyP(x,y|K)\ln \left[\frac{P(x,y|K)}{P(x|K)P(y|K)}\right],$$

$$P(\bm{R},\bm{P}|E)=\frac{\delta\Big{(}E-\sum_{i}\frac{p_{i}^{2}}{2m_{i}}-\Phi(\bm{R})\Big{)}}{\Omega(E)},$$

$$P({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N|E)=\frac{1}{\eta (E)}{\left[E-\Phi ({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)\right]}_{+}^{\frac{3N}{2}-1},$$

$$P({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N|E)=\frac{1}{Z(\beta )}\exp (-\beta \Phi ({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)),$$

$$P({\varphi }_1,\dots ,{\varphi }_N|E)=\frac{1}{Z(\beta )}\exp$$

$$\times {\Omega }_N({\varphi }_1,\dots ,{\varphi }_N),$$

$${\Omega }_N({\varphi }_1,\dots ,{\varphi }_N)=\int d{\mathbf{r}}_1\dots d{\mathbf{r}}_N\prod _{i=1}^N\delta ({\Phi }_i-{\varphi }_i)$$

$${\Omega }_N({\varphi }_1,\dots ,{\varphi }_N)$$

$$=\prod _{i=1}^N{\Omega }_1({\varphi }_i),$$

$$P({\varphi }_i|E)=\frac{1}{Z_1(\beta )}\exp (-\beta {\varphi }_i){\Omega }_1({\varphi }_i).$$

$$P({\varphi }_1|E)=\frac{1}{Z(\beta )}\exp (-\beta {\varphi }_1){\omega }_1({\varphi }_1;E)$$

$${\omega }_1({\varphi }_1;E)=\int d{\varphi }_2\dots d{\varphi }_N\exp$$

$$\times {\Omega }_N({\varphi }_1,\dots ,{\varphi }_N).$$

$$P({\varphi }_1,{\varphi }_2|E)=\frac{\exp \left(-\frac{3N}{2E}({\varphi }_1+{\varphi }_2)\right)}{Z(\frac{3N}{2E})}\times {\omega }_2({\varphi }_1,{\varphi }_2;E),$$

$${\omega }_2({\varphi }_1,{\varphi }_2;E)=\int d{\varphi }_3\dots$$

$$\times {\Omega }_N({\varphi }_1,\dots ,{\varphi }_N).$$

Full text

I Introduction

Melting is a common phenomenon in our daily life, and a clear example of a first-order, or discontinuous, phase transition. Although it is understood in thermodynamic (macroscopic) terms, the transition itself has eluded a complete description from the point of view of microscopic dynamics. In recent decades, computer simulation methods such as molecular dynamics (MD) and Monte Carlo methods, together with modern statistical techniques, have become important theoretical tools to address the challenge of understanding the melting mechanism from an atomistic point of view Jin et al. (2001); Luo and Ahrens (2003); Forsblom and Grimvall (2005).

The melting process in most cases turns out to be an heterogeneous process, since it is commonly triggered on the surface of the crystal. Melting under heterogeneous conditions occurs exactly at the melting temperature $T_{m}$ and is well described by the classical nucleation theory Bai and Li (2005); Luo et al. (2007). However, under ideal conditions, i.e. without surfaces or defects, it is possible to carry out this process homogeneously, in such a way that a superheated, metastable solid state can be reached, far above $T_{m}$.

Eventually, as we increase the energy of the crystal, we reach a critical temperature $T_{LS}$, called the limit of superheating, such that the solid structure cannot exist at a higher temperature than this without being spontaneously transformed into liquid. This process of collapse of the crystalline state is believed to originate from the diffusion of mobile vacancy-interstitial pairs that occurs as the solid surpasses $T_{m}$. In the superheated solid state, the atoms in the crystal can temporarily occupy interstitial sites, creating vacancies Belonoshko et al. (2007) and allowing neighboring atoms to diffuse occupying the vacant sites Bai and Li (2008); Davis et al. (2011); Gallington and Bongiorno (2010). Beyond $T_{LS}$, the solid is no longer in the metastable state but reaches an unstable state, that may be referred to as the critical superheated solid.

In stark contrast with the study of metastable states in classical spin Hamiltonians Binder and Müller-Krumbhaar (1974); Tomita and Miyashita (1992); Moreno et al. (2018), cellular automata Wu et al. (2008), glassy systems Monasson (1995) and other models Pluchino et al. (2004), the statistical mechanical description of the microcanonical superheated solid state is clearly lacking. The case of the supercooled liquid, which is a glassy metastable state, is surprisingly similar to the superheated solid state. There, long-range correlations has been observed, both experimentally and via computer simulation Doliwa and Heuer (2000); Bakai (2002); Bakai and Fischer (2004); Weeks et al. (2007); Szamel and Flenner (2011); Flenner and Szamel (2015), associated with cooperative motion of the atoms. Therefore, a natural question is the existence of similar behavior in the superheated solid state.

Traditionally, correlation is measured through the use of the Pearson correlation metric. This, however, is only able to detect certain types of dependencies between variables. On the other hand, the mutual information metric Cover and Thomas (2006) measures the amount of information that one variable contains about to the other, where through its use it is possible to establish the existence of correlation between the two. The mutual information metric has been previously used in studies of phase transitions and metastability in several systems, in particular spin systems and glasses Wicks et al. (2007); Wilms et al. (2011, 2012); Dunleavy et al. (2015); Alba et al. (2016).

A particularly powerful approach to study the superheated solid state in MD simulation is the so called Z-method Belonoshko et al. (2006). Originally proposed as a practical method to compute melting points in atomistic simulation, it has also become a way to explore the behavior of the superheated solid state and explain the physical meaning of the maximum superheating temperature $T_{LS}$. In this method, the perfect lattice configuration is used as the initial state, and the system is simulated entirely in the microcanonical (NVE) ensemble.

In this work, the behavior of a metastable Lennard-Jones crystal was studied from the observation of atoms and their neighbors immediately before the melting process using molecular dynamics simulations. As stipulated in the Z-method, the calculations were carried out under microcanonical conditions, that is, at fixed energy, volume and number of particles. By means of the atomistic study of superheating in the LJ system, we have managed to contribute to the understanding of the metastable solid phase, presenting evidence of correlations at extended distances that we have measured using the mutual information metric.

This paper continues as follows. In Sec. II, the computational procedure is described in detail. Sec. III presents the main simulation and statistical results, while the paper concludes in Sec. IV with some final remarks.

II Computational Procedure

In this work we performed a number of microcanonical (NVE) MD simulations at different initial temperatures $T_{0}$, equivalent to different total energies of the system through the relation

[TABLE]

with $\Phi_{0}$ the potential energy of the ideal crystal. In our case, $\Phi_{0}$= 556.74 eV. We considered as a starting point for these MD simulations a perfect crystal of argon, corresponding to a face-centered-cubic (FCC) structure, and interacting through a standard Lennard-Jones (LJ) potential,

[TABLE]

where $\varphi(r)$ is the pair potential function,

[TABLE]

The crystal is simulated at high density in order to enhance the superheating effect, using a lattice constant $a=4.2$ Å and the usual LJ parameters for argon, namely $\epsilon/k_{B}=119.8$ K and $\sigma=3.41$ Å. All MD simulations were performed using the LPMD package Davis et al. (2010), with a timestep $\Delta t=$1 fs. For every value of initial temperature $T_{0}$ we used a total simulation time of 80 ps, corresponding to 80000 steps.

As shown in Fig.1, through the use of the Z-method we have obtained an isochoric curve with a melting temperature $T_{m}=5231$ K and superheating temperature $T_{LS}=6695$ K. Also a critical energy $E_{LS}=1976$ eV was obtained, such that it is impossible to administer additional energy without the system melting spontaneously.

The regions of thermodynamic stability are given by $E<$ 1600 eV ($T<T_{m}$), where the solid phase is stable, and $E>$ 1976 eV where the liquid phase becomes stable. For energies between 1600 eV and 1850 eV the system is found in a metastable, superheated solid phase, while in the region between $E$=1850 eV and $E$=2050 eV, the system exists in either of two phases, namely superheated solid or liquid.

In the next section we present the results obtained from MD simulation, which we have used to establish the probability distribution describing the potential energy per atom $\phi$ as a function of total energy $E$. Moreover, we present evidence for the existence of increasing correlation between atoms at long-range distances as the temperature increases, as measured through the mutual information metric.

III Results

Fig. 2 shows a typical microcanonical simulation during which melting is triggered spontaneously. The instantaneous temperature of the system, originally prepared at $T_{0}=12600$ K from an ideal crystalline structure, first drops to approximately 6500 K due to the equilibration of potential and kinetic energies. At this point the system remains in a metastable, superheated solid state. After $t_{w}$=9 ps for this particular simulation, the crystal suddenly collapses and the temperature drops to about 5500 K. For several simulations at the same total energy we obtained random melting times ranging from 1 to 80 ps, which is consistent with previous studies of the kinetics of homogeneous melting Alfè et al. (2011); Davis et al. (2019). In terms of the atomistic dynamics, the picture is as follows. As the system begins to evolve from the ideal crystalline structure, the atoms begin to move away from their stable local potential energy minima. For temperatures above $T_{m}$ the energy per atom is such that thermally-activated point defects start to form Balluffi et al. (2005). This increase of the number of defects and their mobility results in a cooperative motion involving several atoms.

In order to reveal the existence of extreme events and interatomic correlations in the superheated phase, we have studied the statistics of the potential energy per atom $\phi$. In the low-temperature, harmonic approximation of a solid, each atom moves independently as an harmonic oscillator anchored to its equilibrium position. It can be shown that for the microcanonical ensemble in the limit $N\rightarrow\infty$ (see the Appendix for more details), the probability distribution function of $\phi$ for non-interacting particles is given by

[TABLE]

with $\beta=3N/2E$. In the superheated solid, however, we have found that $\phi$ follows a gamma distribution,

[TABLE]

with energy-dependent parameters $k=k(E)$ and $\theta=\theta(E)$. Fig. 3 shows an example of this distribution for $T_{0}$=10000 K, having a most probable energy given by the mode of the gamma distribution, $\phi^{*}=(k-1)\theta$, in this case corresponding to $\phi^{*}$=2.217 eV.

Fig. 4 shows the probability density function $P(\phi|E)$ for different initial temperatures ranging from $T_{0}=$2000 K up to $T_{0}=$12500 K. Here we see that the observed $P(\phi|E)$ is incompatible with Eq. 2 for independent oscillators, because $k$ depends on the total energy, and it is in fact described by

[TABLE]

as shown in the Appendix. Therefore, we can conclude at this point that the potential energies of different atoms are correlated up to some distance to be determined.

While at low temperatures the potential energy values are in fact distributed following a Gaussian distribution, clearly indicating independent behavior compatible with Eq. 2, as the temperature of the system increases the atoms in the lattice break their symmetry. In this phase the solid is in superheated phase, producing a correlation of potential energies between nearest neighbors and even larger distances; these results are consistent with Eqs. 4 and 3. Note that in no case we observe a bimodal distribution of potential energies, as might be expected in the case of nucleation of a liquid phase inside a solid matrix.

In Figs. 5, 6 and 7 an energy of superheating between 1640 and 1936 eV can be observed. It is important to notice that the parameters $k$ and $\theta$ present a dependence on the total energy of the system. From a total energy $E=$1947 eV upwards, the system is present in one of two possible phases: critical superheated solid and liquid. Here we can see two branches where the parameters $k$ and $\theta$ are found in the same phase of the superheated solid. The variation of these parameters with the total energy is an important indicator of the existence of correlations between atoms, where possibly every atom depends on the behavior of its nearest neighbors and even beyond third-nearest neighbors.

Fig. 5 shows the behavior of the shape parameter $k$ as a function of the total energy $E$. As the system does not behave like a collection of independent harmonic oscillators, but is correlated, the effective degrees of freedom are now determined by the “density of states” factor in Eq. 3, governed by the exponent $k(E)$. As seen in the figure, for higher energies $k$ decreases, hence a reduced number of configurations are available for each atom at a given potential energy $\phi$ which again is suggestive of a cooperative motion.

On the other hand, Fig. 6 shows that the scale parameter $\theta$ does not present a continuous behavior as a function of energy. The discontinuity is essentially due to the latent heat per atom, which corresponds to the change in potential energy as the system undergoes the phase transition.

Fig. 7 shows the expected value of potential energy per atom as a function of the total energy $E$, that is, \big{<}\phi\big{>}_{E}, which is given by $k\theta$ in terms of the Gamma parameters. We can clearly see that \big{<}\phi\big{>}_{E} increases linearly with $E$ in both the solid and liquid branches (outside the metastable region), as expected in a system with constant specific heat in their stable phases.

From Eq. 14 for $P(\phi_{1},\phi_{2}|E)$ we performed a non-parametric estimation of the joint distribution using the kernel density estimation (KDE) method Bailer-Jones (2017), considering data from the MD simulations where a set of potential energies $\phi_{1},...,\phi_{N}$ are divided into pairs at given distances $r$. This kernel density estimation is shown in Figs. 8 and 9 for a dataset of potential energies at initial temperatures $T_{0}$=2000 K and $T_{0}$=12500 K, respectively at a distance of $r$ = 7.25 Å. For $T_{0}=$2000 K, the potential energy pairs follow a bivariate Gaussian Bailer-Jones (2017) with almost null correlation. However, for $T_{0}=$12500 K a kind of directionality in the shape of the bidimensional distribution can be seen, in a consistent manner with the marginal Gamma distributions, which confirms the presence of correlations between distant pairs of atoms at high temperatures.

As an indicator of correlation between potential energies of distant atoms, we have computed the mutual information Cover and Thomas (2006), which measures statistical dependence between two arbitrary random variables $x$ and $y$, and is defined by

[TABLE]

where $P(x,y|K)$ is the joint probability distribution of $x$ and $y$, and $P(x|K)$, $P(y|K)$ are the marginal probability distributions. The mutual information $I$ is non-negative, and strictly zero only for statistically independent variables, as in this case $P(x,y|K)=P(x|K)P(y|K)$.

Fig. 10 shows the correlation, as measured by the mutual information diagnostic, between potential energies $\phi_{1}$ and $\phi_{2}$ of a pair of atoms as a function of their distance $r$. It can be seen that for low temperatures, around $T_{0}=$2000 K, the system presents only correlations for short distances, where the peaks are simply indicative of the direct interaction of the atoms with their immediate neighbor shells. However, at larger distances this correlation decays quickly, being almost zero beyond $r$ = 7.75 Å. As the temperature is increased, the system gains energy and becomes metastable once it crosses $T_{m}$, with a clear increase in correlation. In Fig. 11 the mutual information $I(r;T_{0})$ is shown as a function of the initial temperature $T_{0}$ at a fixed distance of $r=$7.25 Å. Interestingly, we can see that the mutual information seems to exhibit a marked, almost linear dependence on $T_{0}$ which saturates close to $T_{LS}$.

At this point we have given a characterization of the superheated solid state in terms of some of its statistical properties. The emerging picture is one of increasing correlation between distant atoms and reduced configurational degrees of freedom per atom as the system approaches $T_{LS}$ from below. In addition, note that it is possible to establish similarities with the glassy state (supercooled liquid), a metastable state occurring at temperatures lower than the melting temperature. In both cases the system is located in points outside the equilibrium phase, “locked” in a dynamical state that does not allow it to transition into the “correct” thermodynamical phase, and presenting nontrivial atomic kinetics.

Despite the fact that the case of superheated solid could seem opposite to the glassy state, apparently these transitions occur in a similar way, both involving cooperativity between their atoms, as the movement of one atom depends on the movement of the others. In both cases there is the presence of spatial correlations near the phase transition, in our case near the melting point. This cooperative movement implies a reduction of the effective degrees of freedom, as shown in Figs. 5 and 11. As the energy increases, correlations increase at large distances between the potential energies of the atoms, as shown in Fig. 11, while in turn the number of effective degrees of freedom is reduced, as can be seen in Fig. 5, because these correlations limit the amount of allowed movements of the atoms.

IV Concluding remarks

Using microcanonical molecular dynamics simulations we have characterized the statistics of single-atom potential energies in the metastable solid prior to melting, according to a gamma distribution with energy-dependent parameters $k$ and $\theta$, which is consistent with correlated behavior. Furthermore, we have presented strong evidence confirming the emergence of long-range spatial correlations of potential energy in the critical superheated solid. These correlations may arise due to the cooperative dynamical processes occurring in the crystal after the creation of mobile vacancy-interstitial pairs.

Acknowledgements

This work is supported by FONDECYT grant 1171127 and Anillo ACT-172101. Computations were performed using the supercomputing infrastructure of FENIX, Materials Science Group, UNAB Physics Department (http://www.matbio.cl/fenix).

APPENDIX: Microcanonical distributions of potential energy

In the microcanonical ensemble, the probability density of microstates is given by

[TABLE]

with $\bm{R}=(\bm{r}_{1},\ldots,\bm{r}_{N})$, $\bm{P}=(\bm{p}_{1},\ldots,\bm{p}_{N})$ and where $\Omega(E)$ the density of states. By integration over $\bm{P}$ it can be shown Pearson et al. (1985); Ray (1991); Davis (2011) that

[TABLE]

which in the limit $N\rightarrow\infty$ reduces to a canonical distribution,

[TABLE]

with $\beta=3N/2E$. Therefore, the probability of observing a set of potential energies $\Phi_{1}=\phi_{1}$, $\Phi_{2}=\phi_{2}$, …, $\Phi_{N}=\phi_{N}$ at total energy $E$ is given by

[TABLE]

where

[TABLE]

is the joint configurational density of states. For non-interacting atoms we have that $\Phi_{i}=\Phi_{i}(\bm{r}_{i})$ and therefore

[TABLE]

hence $P(\phi_{1},\ldots,\phi_{N}|E)$ factorizes and $\phi_{i}$ itself follows a canonical distribution

[TABLE]

Here $\Omega_{1}(\phi)=\int d\bm{r}\delta(\phi-\Phi_{1}(\bm{r}))$ plays the role of a “density of states” factor which is independent of the total energy $E$. On the other hand, for interacting atoms the probability density $P(\phi_{1}|E)$ will be given by the marginalization of Eq. 9,

[TABLE]

where

[TABLE]

Now the “density of states” factor $\omega_{1}$ is energy-dependent. Similarly, for the joint distribution $P(\phi_{1},\phi_{2}|E)$ we obtain by marginalization

[TABLE]

where $\omega_{2}(\phi_{1},\phi_{2};E)$ is defined by

[TABLE]