Molecular dynamics simulation of homogeneous nucleation in a superheated Lennard-Jones crystal

Molecular dynamics simulation of homogeneous nucleation in a superheated Lennard-Jones crystal

Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx Contents lists available at ScienceDirect Journal of Non-Crystalline Solids journal homepage: w...

920KB Sizes 0 Downloads 13 Views

Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx

Contents lists available at ScienceDirect

Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/locate/jnoncrysol

Molecular dynamics simulation of homogeneous nucleation in a superheated Lennard-Jones crystal V.G. Baidakov , A.O. Tipeev ⁎

Institute of Thermal Physics, Ural Branch of the Russian Academy of Sciences, Ekaterinburg 620016, Russia

ARTICLE INFO

ABSTRACT

Keywords: Crystal superheating Homogeneous nucleation Classical theory of nucleation Molecular dynamics simulation Lennard-Jones crystal

We present the results of a molecular dynamics study of homogeneous nucleation of liquid drops in a superheated Lennard-Jones (LJ) crystal at positive and negative pressures. The mean lifetime method has been used to determine the nucleation rate, calculate the Zeldovich nonequilibrium factor and the diffusion coefficient of nuclei, and determine the contribution of elastic stresses to the work of formation of a critical nucleus. The data obtained are compared with classical nucleation theory (CNT). A considerable discrepancy between data of MD simulation and CNT (up to 70 orders in the nucleation rate) has been established. With a satisfactory agreement between theory and simulation in the rate of the nucleus transition through the critical size, an essential discrepancy is observed in the work of formation of a critical nucleus. The neglect of elastic stresses in the nucleation work improves the agreement between the results of theory and simulation. The effect of the size dependence of the surface free energy of a critical nucleus on the work of its formation is estimated.

1. Introduction A first-order phase transition presupposes the existence of metastable states on both sides of the line of phase equilibrium. The capacity of liquids for high superheatings is well known [1,2]. It has proved to be much more difficult to register a crystal superheating. The asymmetry of a crystal–liquid transition is explained by the surface effect. Since the surface energy of a crystal at the boundary with vapor (or vacuum) always exceeds that of a liquid at the boundary with vapor, in melting the formation of a liquid layer at the crystal surface proves to be more advantageous thermodynamically. Besides, in melting processes the superheating of a crystal phase is retarded by elastic stresses, nonequilibrium structure defects and inner boundaries. As distinct from melting, spontaneous crystallization proceeds at considerable supercoolings of a liquid phase. The superheating of a crystal phase was registered by Haykin and Bene [3] under heating of a lead monocrystal by an electric current and simultaneous cooling from the surface. Another way of obtaining a metastable crystal is shock-wave experiments [4]. In this case supersaturation is created by stretching a crystal phase. The small volumes of the samples under investigation and the characteristic time intervals of action make it possible to hope that spontaneous nucleation can be achieved in such experiments. In shock-wave experiments the loss of stability of a crystal phase is



possible with respect to the formation of both a liquid and a vapor phase. The melting line has a metastable extension beyond the triple point into the region of negative pressures [5]. Here, in equilibrium are a crystal and a liquid phase, each being metastable against the disturbance of their continuity. The metastable equilibrium of a liquid and a crystal phase ends at the point of meeting of the melting line and the spinodal of a stretched liquid (at the endpoint of the melting line). This means that at a negative pressure there is a certain boundary temperature, below which the appearance in a crystal of a supercritical liquid-phase nucleus becomes impossible [6]. Methods of MD simulation open up new possibilities in studying spontaneous nucleation in a superheated (stretched) crystal. The periodic boundary conditions of MD models make it possible to get rid of the surface melting of a crystal, and the small sizes of models ensure the absence of considerable defects and the possibility of high superheatings. Limiting crystal superheatings at pressures close to atmospheric were investigated in Refs. [7–9]. By the data of Lu and Li [7], the mass formation of liquid-phase nuclei in an LJ crystal is observed at a temperature Tn ≈ 1.2Tm, where Tm is the melting temperature. For pure metals and simple compounds Tn varies within the limits (1.08–1.43) Tm and has a weak dependence on the heating rate [8]. In a wide range of pressures, including negative ones, information on limiting superheatings of a crystal phase has been obtained in Refs. [10–13].

Corresponding author at: Institute of Thermal Physics, Ural Branch of the Russian Academy of Sciences, Amundsen st, 107a, Ekaterinburg 620016, Russia. E-mail address: [email protected] (V.G. Baidakov).

https://doi.org/10.1016/j.jnoncrysol.2018.10.007 Received 15 August 2018; Received in revised form 28 September 2018; Accepted 7 October 2018 0022-3093/ © 2018 Elsevier B.V. All rights reserved.

Please cite this article as: Baidakov, V.G., Journal of Non-Crystalline Solids, https://doi.org/10.1016/j.jnoncrysol.2018.10.007

Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx

V.G. Baidakov, A.O. Tipeev

In the present paper MD simulation is used to study the kinetics of nucleation of a liquid phase in a superheated LJ crystal at positive and negative pressures. The data obtained are analyzed in the framework of CNT. The effect of elastic stresses, anisotropy, etc. on the process of nucleation is considered. The paper is constructed as follows. Section 2 describes the fundamentals of CNT with regard to a crystal–liquid phase transition. The simulation details are described in section 3. In section 4 the results of simulation are compared with CNT. Section 5 is a conclusion.

W = A + [(p

(p

2Kl + 4µs

0 s0

3 (p

pm )(

m / sm ) Vl

(5)

3

pm )

m sm

+e

2

(6)

According to CNT, the kinetic factor in Eq. (1) is

B=Z D

dn dR

(7)

where Z∗ is the Zeldovich nonequilibrium factor, D∗ is the nucleus diffusion coefficient at the top of the activation barrier, and n is the number of particles in the nucleus. The Zeldovich nonequilibrium factor is determined by the curvature of the potential barrier and in a first approximation may be taken equal to its value in an isotropic medium, i.e.

d 2W dR2

Z =

1 2 kB T

1/2

1/2

=2

kB T

(8)

To find D∗, it is necessary to adopt a model of the determining process close to the surface of a growing nucleus. Let us believe that the factor which retards the growth of a liquid drop in an elastic medium is the kinetics of the exchange of particles between the solid and the liquid phase. Assuming the melting mechanism (the growth of a liquid drop) to be normal, by analogy with crystallization theory [1], for the coefficient of diffusion of nuclei along the axis of its size R we have

D =

3 kB T 2 R2

(9)

where κ is the kinetic coefficient that characterizes the rate of particle exchange between the solid and the liquid phase. After substitution of (8) and (9) into (7) for the kinetic factor we obtain

( kB T )1/2

B = 12

(10)

The dependence of κ on the thermodynamic state parameters will be found with the use of the theory of absolute reaction rates [1].

=

2

0

kB T

4/3

exp( U / kB T )

(11)

Here v0 is the particle vibration frequency at the nucleus surface and U is the activation energy, which is taken equal to the melting heat calculated per particle. The vibration frequency v0 will be evaluated from the relation hv0 ≈ kBT, where h is the Planck constant. An alternative to CNT is the theory of transient state [20], which assumes the inertial motion of a nucleus over the activation barrier. This theory gives an estimate from above for the value of the kinetic factor in Expression (1) and according to [21].

(2)

where μs is the shear modulus of a solid, Kl and Ks are the bulk moduli, and Δρ0 = ρs0 − ρl0 is the difference of the densities of a solid and a liquid phase. All the quantities in Formula (2) are taken at p = 0 and temperature T. Taking into account the small compressibility of the liquid and solid phases, for the sum of the volume and the molecular work we have

WV + Wµ = (p

+e

sm

16

W =

where ρ is the number of particles in a crystal volume unit, B is the kinetic factor, W∗ is the work of formation of a critical nucleus, kB is the Boltzmann constant, and T is the temperature. The factor ρ exp (−W∗/kBT) in the right-hand side of (1) is the average number of critical nuclei in a unit of volume, and the factor B corresponds to the average rate of transition of an incipient phase through the critical size. The work of formation of a liquid-phase nucleus W in an isotropic elastic medium is determined by the sum of the surface Wγ, the volume WV, the molecular Wμ, and the deformation We work. At given T, p and nucleus volume Vl the elastic energy, as well as the surface one, depends on the shape of a liquid drop and in the general case must be determined from the condition of minimum W. As distinct from a fluid phase, the spinodal shape of a nucleus in an elastic medium is not always most advantageous in terms of energy [15–17]. It is optimum at high superheatings (stretchings) [17], but at a small degree of metastability leads to an unphysical solution, namely, elastic thermodynamic hysteresis [18]. The present paper considers superheatings (stretchings) close to the boundary of essential instability of a crystal phase, that is the spinodal. Here, as is shown in what follows, the results of molecular dynamics simulation point to the fact [19] that the shape of nuclei is close to spherical. The appearance of a liquid-phase nucleus is accompanied by deformations of a solid. The calculation of elastic energy is connected with the determination of the fields of shears and stresses in such a twophase system. In elasticity theory the initial state is usually an undeformed solid (p = 0), which at a certain temperature may prove to be thermodynamically unstable. According to [17,18], for specific elastic energy of internal stresses arising at a contact of two phases which are separated by a spherical boundary and differ in intrinsic deformation we have

e=

m

pm )

Substitution of (5) into (4) gives the work of formation of a critical nucleus

(1)

K Kl + s p Ks Kl

(4)

+ e] Vl

2

R =

CNT considers the fluctuation growth of the fragments of an incipient phase as the process of their diffusion along the size axis R in the field of thermodynamic forces [14]. The average number of viable nuclei that form in a unit of time in a liquid volume unit, i.e. the nucleation rate, may be presented as follows:

2Kl µs

m / sm )

A critical nucleus is in unstable equilibrium with the environment. From the condition of such equilibrium and (4), for the radius of a critical nucleus we have

2. Classical nucleation theory

J = B exp( W / kB T )

pm )(

(3)

B=

Here the values of pm, ρsm and Δρm refer to the melting line and temperature T. With allowance for (2) and (3), the minimum work of formation of a liquid-phase nucleus in an elastic medium is written as follows:

2 R kB T m

1/2

(12)

where m is the particle mass. In the absence of friction, the inertial motion of a nucleus leads to the independence of a flow over the barrier from the shape of the barrier top. 2

Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx

V.G. Baidakov, A.O. Tipeev

3. Molecular dynamics simulation To study the kinetics of homogeneous nucleation in a superheated LJ crystal, use was made of the mean lifetime method [22]. This method was used previously in experimental investigations [22,23] and computer simulation [24] of bubble nucleation. The systems under investigation contained from 2 048 to 500 000 interacting particles located at the sites of a face-centered cubic (FCC) lattice. The particles interacted with the help of the LJ potential cut off at a distance of 6.58 σ, where σ is a potential parameter. Periodic boundary conditions were imposed on the boundaries of a system which had the shape of a cube. MD calculations were carried out in an NVT ensemble, and the phase decay proceeded in NVE conditions. The temperature in the NVT ensemble was controlled with the help of a Nose-Hoover thermostat [25]. In integrating the particle motion equations use was made of the accelerated Verlet algorithm [26]. An integration time step Δ t was 0.002 τLJ, where τLJ = (mσ2/ε)1/2, and ε is a parameter of the LJ potential. Further, unless otherwise specified, all the numerical values of quantities are given in dimensionless units. The reduction parameters are the parameters of the LJ potential σ and ε, the particle mass m, and the Boltzmann constant kB. The unit of reduced temperature is ε/kB, pressure ε/σ3, density σ−3, and nucleation rate σ−4(m/ε)−1/2. In MD simulation the entry into a metastable region was realized by an isothermal stretching of a crystal. A transition to the state under investigation (p, T = const), at which the mean expectation time for the phase decay was determined, was carried out in two steps. Initially, the system was transferred to a certain intermediate, slightly metastable state (p', T = const), in which the crystal lifetime exceeded considerably the characteristic times of MD simulation (104–106 MD steps). Then an ensemble of initial particle configurations was built. For this purpose, the coordinates and the velocities of all particles in the cell were recorded on the phase trajectory of the system every 2∙104 MD steps. The configurations selected were initial microscopic states, from which a transition to the state under investigation (p, T) was realized by means of linear scaling of the coordinates of all particles in the system (by decreasing the density of the system). To reduce the relaxation time and retain the thermodynamic parameters, a temperature correction was carried out within 5∙103 MD steps after the transition. The spontaneous formation of a new-phase nucleus in a homogeneous metastable system is a random event. To describe a flow of such events, it is necessary to determine the distribution function and its first instants. The time before the beginning of a phase decay τR was registered in MD experiments. The instant of the beginning of the destruction of a metastable phase was registered by changes in the temperature and pressure caused by nucleation. A deviation of 2 % of p and T from the registered average values was taken as the boundary value. The time τR is the sum of the times of expectation for the first viable nucleus τ and its subsequent growth τp to the instant of registration of a phase transition. Simulation has shown that a phase decay in the LJ systems under investigations is always initiated by one nucleus and τp < < τ. Therefore we assumed τ ≈ τR. The kinetics of nucleation of a liquid phase in an LJ crystal has been investigated at reduced temperatures T = 0.55, 0.7, 0.85 and 1.0. A temperature T = 0.7 is close to that of the triple point (Tt= 0.692) of an LJ system, and T = 0.55 to the temperature of the melting line endpoint (TK= 0.529). At these two temperatures spontaneous nucleation in a crystal took place in the region of negative pressures. From N = 50 to 200 independent events of spontaneous crystal melting were registered in every given thermodynamic state (p, T = const). All in all, 20 thermodynamic states have been investigated. The maximum expectation time for the appearance of the first supercritical nucleus reached 107 MD steps, which was equal to 0.05 μs. Fig. 1 shows the distribution of crystal melting events N in time τ. To exclude the effect of a transient process caused by a pressure decrease, a certain value τ′, which wittingly exceeded the duration of the

Fig. 1. Histogram of distribution of melting events for an LJ crystal: T = 0.7, ρ= 0.870, N = 204, = 37.4, and N = 2048. A solid line shows the Poisson distribution.

transient process, was subtracted from all the values of the expectation time τ obtained. The results of simulation with τ < τ′ were discarded. As is evident from Fig. 1, the distribution of expectation times τ satisfies the Poisson distribution (the smooth curve). This justifies the exclusion from consideration time intervals smaller than τ′ because a distribution of the Poisson type has the property of independence of the probability of the onset of an individual event from the beginning of the time reading. In every thermodynamic state, data on τ are used to determine the for spontaneous melting and calculate the mean expectation time nucleation rate J = ( V ) 1, where V is the crystal volume (Fig. 2). In reduced units the minimum and the maximum value of J were 10−9 and 10−5, in dimensional units, respectively, 1031 and 1035 s−1 m−3. An advance into the region of higher nucleation rates (J > 10−5) was limited by the conditions of the process stationarity and the minimum possible number of particles in the model (N = 2048). The latter limitation was connected with the selected cutoff radius of the LJ potential equal to rl = 6.58. Calculations of the nucleation rate at J < 10−9 are hampered by a considerable increase in the computer time.

Fig. 2. Logarithm of nucleation rate as a function of density: 1 – T = 0.55, 2–0.7, 3–0.85, and 4–1.0. 3

Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx

V.G. Baidakov, A.O. Tipeev

Fig. 4. Expectation time for the appearance of a liquid-phase nucleus as a function of the number of particles contained in it: T = 0.7, ρ = 0.885, =472, n∗= 295, and N = 32 000. A solid line shows calculation by Eq. (14).

Fig. 3. Root-mean-square variation of the number of liquid-phase particles in a nucleus as a function of time.

MD simulation makes it possible to determine the kinetic factor B (Eq. (7)). The nucleus diffusion coefficient close to the critical size was calculated as follows [27].

D =

1 2

n2 ( )

(13)

where Δ n(τ) = [n(τ) − n(0)] is the variation of the number of liquidphase particles in a nucleus during the timeτ. In every thermodynamic state 5 random configurations with n close to n∗, n = n∗ ± 20, were chosen and acted as the starting ones for new MD configurations with changed values of the initial particle velocities. The beginning of the time reading was redetermined whenever n was equal to n∗. The calculations ceased at n < 100 or n > 2n∗. Fig. 3 illustrates the dependence 〈Δn2(τ)〉 close to the critical size. The uncertainty of determination of the nucleus diffusion coefficient by the data on 〈Δn2(τ)〉 from Eq. (13) essentially depends on the temperature. For T = 1.0 it reaches 40%, which is connected with considerable fluctuations of the number of particles in a nucleus. The number of particles in a critical nucleus n∗ and the Zeldovich nonequilibrium factor Z∗ were determined by the method of the mean time of the first passage across a given boundary [28]. For this purpose, the coordinates of all the particles in the cell were retained in every act of fluctuation melting with an interval of 103 MD steps. A search for crystal-like particles (particles that may be referred to a crystal phase) was realized in the configurations retained with the use of the algorithm q6 [29]. The parameters of the algorithm q6 were chosen the same as in searching for crystal structures in a supercooled liquids [30]. In every configuration the number of particles belonging to the liquid phase n was determined as the difference of the total number of particles in the system N and the number of particles belonging to the crystal phase Ns. In all thermodynamic states expectation times for the appearance of the first liquid-phase nucleus n were averaged over 30 melting events. The average value of the nucleus expectation time 〈τ〉 with n particles was approximated by the function

(n ) =

2

{1

erf[

1/2Z

(n

n )]}

Fig. 5. Isotherms of an LJ crystal in pressure–density coordinates: 1 – T = 0.55, 2–0.7, 3–0.85, and 4–1.0. Dashed lines show branches of the melting line and dash-dot lines the spinodal of a superheated (stretched) liquid (ab) and a superheated (stretched) crystal (cd). Tt stands for the triple point and Km for the endpoint of the melting line.

4. Results and their comparison with CNT Fig. 5 shows in pressure-density coordinates the isotherms on which the liquid-phase nucleation rate, as a function of the pressure, is shown in Fig. 3. The uncertainty of determination of J is comparable with the size of the dots in the figure. In comparing the results of MD simulation and CNT, it is necessary to possess a complex of thermodynamic and kinetic properties of the system under investigation. In our calculations the surface free energy at a crystal–liquid drop interface was taken equal to its value at a plane crystal–liquid interface [31], and the parameters pm, ρsm, ρlm and ρs0, ρl0 are determined in Refs. [5, 32]. For a defect-free LJ FCC crystal under a hydrostatic pressure the shear μs and bulk Ks moduli are calculated in Ref. [13]. The bulk modulus of the liquid phase Kl = ρ(∂p/∂ρ)T has been found from the equation of state of an LJ fluid [33]. The elastic energy e has been determined from Eq. (2) by data of MD simulation. With an isothermal increase in the crystal volume of 1 or 2% the value of e decreases by 6 to 19%. The contribution of the elastic energy to the motive force of the phase decay process (μl − μ) decreases with decreasing temperature. This contribution is 30 to 36% at T = 1.0 and 4–5% at T = 0.55. According to Eq. (5), in the range of state parameters under investigation the radii of critical liquid-phase nuclei are R∗ = 4.5 − 9.2. A smaller value of R∗ corresponds to a lower temperature. On the

(14)

where erf(x) is a function of errors. For one state the dependence of 〈τ〉 on the number of particles in a nucleus is presented in Fig. 4. An increase in the temperature leads to an increase in the number of particles in a critical nucleus. If at T = 0.7 a nucleus contains n∗ = 295 particles, then at T = 1.0 it already contains 510 particles (J = 10−8).

4

Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx

V.G. Baidakov, A.O. Tipeev

isotherm the value of R∗ decreases with an increase in the stretching. If it is assumed that the density of the liquid phase in a critical nucleus is equal to that on the melting line, then for the number of particles in a liquid drop we have n∗ = 500 − 2000. These values are 3 or 4 times larger than the values obtained in MD simulation by the method of first passage beyond a given boundary. The neglect of the elastic energy in Eq. (5) reduces this discrepancy 1.3 to 1.7 times. The values of the Zeldovich nonequilibrium factor Z∗ and the nucleus diffusion coefficient D∗ differ from those calculated from CNT formulas by (5–7) % and (20−30) %, respectively. According to Eqs. (10) and (12), in the range of states under investigation the values of the kinetic factor are 1.7 to 3.9 and 2.6 to 9.5, respectively. The neglect of the elastic energy in the equations mentioned above leads to an increase in B of (5–10) %. The results of a direct MD calculation of B differ from the data of Eqs. (10) and (12) by approximately an order of magnitude. Since the pressure (density) dependence of the nucleation rate is very strong, the differences indicated have only a slight effect on the value of the crystal superheating (stretching). The most highly elastic stresses tell on the work of formation of a critical nucleus. Taking into account elastic stresses increases the stability of metastable crystals against phase decay. The reduced value of the activation barrier W∗/kBT (the Gibbs number) calculated from Eq. (6) is 52–57 at T = 0.55 and increases to 147–204 at T = 1.0. The neglect of the elastic energy in Eq. (6) reduces the value of W∗/kBT to 48–83. In the systems investigated (N = 32 000–500 000) so high activation barriers cannot be overcome in characteristic times of MD simulation (107 MD steps). In MD simulation of nucleation in superheated and supercooled liquids, and in supersaturated vapor, the value of W∗/kBT does not usually exceed 20–25 [24,33,34]. Thus CNT gives values of J from 16 to 70 orders of magnitude lower than those obtained in MD simulation by the mean lifetime method. The higher the temperature, the larger the deviations. If the elastic energy in the work of formation of a critical nucleus is neglected (6), the discrepancy in J will be from 14 to 27 orders of magnitude. A visualization of the process of a crystal phase decay, as well as the decay thermodynamic and kinetic parameters (R∗, Z∗, D∗, n∗) points to the fact that at high superheatings (stretchings) it proceeds by the mechanism of nucleation and growth of liquid-phase nuclei. With a satisfactory agreement between CNT and the results of MD simulation in the preexponential factor there in a considerable discrepancy in the value of the activation barrier W∗. This cannot be connected with elastic stresses because even their absence gives no way of reducing the height of the activation barrier to values that the system can overcome in the process of MD simulation. Another possibility of reducing the height of W∗ is taking into account the size dependence of the surface free energy of a critical nucleus. If it is assumed that the rupture of a crystal phase proceeds by the nucleation mechanism and is described by CNT, then MD simulation data make it possible to estimate from Eqs. (1), (6) and (10) the effective crystal–critical nucleus surface free energy γe. For the case e = 0 the results of calculations of the dependence of γe on the surface curvature are presented in Fig. 6. The value of γe at e = 0 is 30–35% smaller than at a plane interface. Taking into account the elastic stresses increases this discrepancy. At T = 0.85 and 1.0 the deviations are already equal to 50%, and the character of the R–dependence of γe changes. The value of γe(R) becomes an increasing rather than a decreasing function of the curvature. The small excess of the values of γe at a temperature T = 0.55 over the values corresponding to T = 0.7 may be connected with the above-mentioned character of the crystal phase rupture close to the endpoint of the melting curve.

Fig. 6. Effective surface free energy of a liquid-phase nucleus as a function of the curvature of the dividing surface: 1 – T = 0.55, 2–0.7, 3–0.85, and 4–1.0. Arrows show the values of the surface free energy at a plane crystal–liquid interface.

coefficient, Zeldovich nonequilibrium factor, and critical nucleus size) have been made along four isotherms. Two isotherms had temperatures exceeding that of the triple point. The temperatures of the other two isotherms were close to those of the triple point and the endpoint of the melting line. At p = 0 and J = 10−7 the limiting superheating of a crystal phase was 1.2Tm, which is in good agreement with the data of Lu and Li [7]. To the spinodal crystal state on the zero isobar corresponds a superheating temperature Tsp = 1.29Tm. A pressure decrease and increase lead to insignificant changes in the relative value of the superheating of a crystal phase (Tn /Tm = 1.19 − 1.28). At the approach to the pressure of the endpoint of the melting line the boundaries of the limiting superheating and the spinodal approach each other [5,12,35]. The results of MD simulation point the fact that at high superheatings (stretchings) the decay of a crystal proceeds by the mechanism of nucleation and growth of competitive-phase nuclei. CNT gives a correct evaluation of the rate of the transition of nuclei through the critical size, but leads to considerable differences in the evaluation of the average number of critical nuclei in a unit of volume from data of MD simulation. In solids, as distinct from gases and liquids, the formation of a new phase is connected with the appearance of the elastic energy of internal stresses. We have calculated the contribution of the elastic energy to the work of formation of a critical (liquid) nucleus assuming it to be spherical. Taking into account the elastic energy in CNT formulas leads to an increase in the discrepancies between theory and data of MD simulation in the nucleation rate (up to 70 orders of magnitude). So considerable discrepancies are not observed for crystallization of supercooled liquid [30] and condensation of supersaturated vapor [28,34]. At cavitation in a stretched LJ liquid the discrepancies between data of CNT and MD simulation in J reach 20 orders of magnitude (T = 0.35) and are equal to 9–10 orders close to the temperature of the triple point [33,36]. As shown in Refs. [33, 36], this discrepancy is connected with the neglect in CNT of the dependence of the surface tension on the nucleus size. Assuming that at melting the discrepancy between CNT and MD simulation is also a result of the size effect, we have calculated the effective surface free energy of liquid phase critical nuclei γe in a superheated crystal. It has been found that at e = 0 the differences of γe from the surface free energy of a crystal–liquid plane interface are from 30 to 35%. Thus the investigation conducted has shown that in the absence of free surfaces the spontaneous melting of a crystal phase, at least qualitatively, may be described by CNT. A more detailed comparison of

5. Conclusion We have investigated the kinetics of a phase decay of a superheated LJ FCC crystal in MD simulation. Calculations of the main characteristics of the phase decay process (nucleation rate, nucleus diffusion 5

Journal of Non-Crystalline Solids xxx (xxxx) xxx–xxx

V.G. Baidakov, A.O. Tipeev

CNT and the results of MD simulation requires a direct calculation of the height of the activation barrier and the size dependence of the crystal–liquid drop surface free energy.

transitions in solids, Sov. Phys. Solid State 29 (1987) 732. [17] E.A. Brener, S.V. Iordanskii, V.I. Marchenko, Elastic effects on the kinetics of a phase transition, Phys. Rev. Lett. 82 (1999) 1506–1509, https://doi.org/10.1103/ PhysRevLett.82.1506. [18] A.L. Roitburd, D.E. Temkin, Plastic deformation and thermodynamic hysteresis at phase changes in solids, Fiz. Tverdogo Tela. 28 (1986) 775–784. [19] V.G. Baidakov, A.O. Tipeev, Nucleation of liquid droplets and voids in a stretched Lennard-Jones fcc crystal, J. Chem. Phys. 143 (2015) 124501, https://doi.org/10. 1063/1.4931108. [20] S. Glesston, K. Leidler, G. Eiring, Theory of Rate Processes, McGraw-Hill, New York, 1941. [21] A.V. Prokhorov, Equation of nucleation in phase space. Boiling-up of non-volatile liquids, Dokl. Akad. Nauk SSSR 239 (1978) 1323–1326. [22] V.P. Skripov, Metastable Liquids, Wiley, New York, 1974. [23] V.P. Skripov, et al., Thermophysical Properties of Liquids in the Metastable (Superheated) State, Gordon and Breach Science Publishers, New York, 1988. [24] V.G. Baidakov, S.P. Protsenko, Computer simulation of nucleation in a liquid under tension, Dokl. Phys. 49 (2004) 69–72, https://doi.org/10.1134/1.1686871. [25] D. Frenkel, B. Smit, Understanding Molecular Simulations: From Algorithms to Applications, Academic Press, London, 1996. [26] L. Verlet, Computer “experiments” of classical fluids. I. Thermodynamical properties of Lennard-Jones molecules, Phys. Rev. 159 (1967) 98–103, https://doi.org/10. 1103/PhysRev.159.98. [27] S. Auer, D. Frenkel, Numerical prediction of absolute crystallization rates in hardsphere colloids, J. Chem. Phys. 120 (2004) 3015–3029, https://doi.org/10.1063/1. 1638740. [28] J. Wedekind, R. Strey, D. Requera, New method to analyze simulations of activated processes, J. Chem. Phys. 126 (2007) 134103, https://doi.org/10.1063/1.2713401. [29] W. Lechner, S. Dellago, Accurate determination of crystal structures based on averaged local bond order parameters, J. Chem. Phys. 129 (2008) 114707, https:// doi.org/10.1063/1.2977970. [30] V.G. Baidakov, A.O. Tipeev, Crystal nucleation and the solid–liquid interfacial free energy, J. Chem. Phys. 136 (2012) 174510, https://doi.org/10.1063/1.3678214. [31] V.G. Baidakov, S.P. Protsenko, A.O. Tipeev, Temperature dependence of the crystal–liquid interfacial free energy and the endpoint of the melting line, J. Chem. Phys. 139 (2013) 224703, https://doi.org/10.1063/1.4837695. [32] V.G. Baidakov, Z.R. Kozlova, S.P. Protsenko, Thermal and caloric equations of state for stable and metastable Lennard-Jones fluids: I. Molecular dynamics simulations, Fluid Phase Equilib. 263 (2008) 55–63, https://doi.org/10.1016/j.fluid.2007.09. 019. [33] V.G. Baidakov, Spontaneous cavitation in a Lennard-Jones liquid: molecular dynamics simulation and the van der Waals-Cahn-Hilliard gradient theory, J. Chem. Phys. 144 (2016) 074502, https://doi.org/10.1063/1.4941689. [34] B.J. Block, et al., Curvature dependence of surface free energy of liquid drops and bubbles: a simulation study, J. Chem. Phys. 133 (2010) 154702, https://doi.org/10. 1063/1.3493464. [35] A.Yu. Kuksin, G.E. Norman, V.V. Stegailov, Molecular dynamics modelling of lifetime and decay of metastable crystals under superheating of stretching, Computational Physics: Proceedings of the Joint Conference of ICCP6 and CCP2003, Beijing Rinton Press Inc., 2004, p. 126. [36] V.G. Baidakov, K.S. Bobrov, Spontaneous cavitation in a Lennard-Jones liquid at negative pressures, J. Chem. Phys. 140 (2014) 184506, https://doi.org/10.1063/1. 4874644.

Acknowledgements The investigation has been conducted at the expense of a grant of the Russian Science Foundation (project № 18-19-00276). No potential conflict of interest was reported by the authors. References [1] V.P. Skripov, V.P. Koverda, Spontaneous Crystallization of Supercooled Liquids, Nauka, Moscow, 1984. [2] K.F. Kelton, A.L. Greer, Nucleation in Condensed Matter: Applications in Materials and Biology, Pergamon, Amsterdam, 2010. [3] S.E. Haykin, N.P. Bene, On the phenomenon of superheating of solids, Dokl. Akad. Nauk SSSR 23 (1939) 31–35. [4] G.I. Kanel, V.E. Fortov, S.V. Razorenov, Shock waves in condensed-state physics, Uspekhi Fiz. Nauk. 177 (2007) 809–830, https://doi.org/10.3367/ufnr.0177. 200708a.0809 (Phys. Uspekhi 50 771. 10.1070/pu2007v050n08abeh006327). [5] V.G. Baidakov, S.P. Protsenko, Singular point of a system of Lennard-Jones particles at negative pressure, Phys. Rev. Lett. 95 (2005) 015701, https://doi.org/10.1103/ PhysRevLett.95.015701. [6] V.G. Baidakov, A.O. Tipeev, Mechanical instability and nucleation in a LennardJones fcc crystal and limiting stretching, Chem. Phys. Lett. 643 (2016) 6–9, https:// doi.org/10.1016/j.cplett.2015.10.079. [7] K. Lu, Y. Li, Homogeneous nucleation catastrophe as a kinetic stability limit for superheated crystal, Phys. Rev. Lett. 80 (1998) 4474–4477, https://doi.org/10. 1103/PhysRevLett.80.4474. [8] S.N. Lou, et al., Maximum superheating and undercooling: Systematics, molecular dynamics simulations, and dynamic experiments, Phys. Rev. B 68 (2003) 134206, https://doi.org/10.1103/PhysRevB.68.134206. [9] V.G. Baidakov, A.E. Galashev, V.P. Skripov, Stability of superheated crystal in the molecular-dynamical model of argon, Fiz. Tverdogo Tela. 22 (1980) 2681. [10] S.-N. Lou, L. Zheng, A. Strachan, D.C. Swift, Melting dynamics of superheated argon: Nucleation and growth, J. Chem. Phys. 126 (2007) 034505, https://doi.org/ 10.1063/1.2424715. [11] J. Wang, J. Li, S. Yip, S. Phillot, D. Wolf, Mechanical instabilities of homogeneous crystals, Phys. Rev. B 52 (1995) 12627–12635, https://doi.org/10.1103/PhysRevB. 52.12627. [12] A.Yu. Kuksin, G.E. Norman, V.V. Stegailov, The phase diagram and spinodal decomposition of metastable states of Lennard-Jones system, High Temp. 45 (2007) 37–48 http://www.ihed.ras.ru/norman/eLibrary/HighTemperature_45_37.pdf. [13] V.G. Baidakov, A.O. Tipeev, Ideal and ultimate tensile strength of a solid body, High Temp. 56 (2018) 184–192, https://doi.org/10.1134/S0018151X18020013. [14] Ya. Zeldovich, Theory of formation of a new phase, Cavitation Zh. Exp. Teor. Fiz. 12 (1942) 525–538. [15] I.M. Lifshitz, L.S. Gulida, On local melting theory, Dokl. Akad. Nauk SSSR 87 (1952) 377–380. [16] V.I. Motorin, On the absence of elastic thermodynamic hysteresis at phase

6