Self-instability of finite sized solid-liquid interfaces

In solid-liquid systems, macroscopic solids lose their equilibrium and melt in a manner that results in overall movement of the solid-liquid interface. This phenomenon occurs when they are subjected to temperature gradients or external stress, for example. However, many experiments suggest that the melting of nano- and micro-sized metallic nuclei follows a different process not described by traditional melting theory. In this paper, we demonstrate through simulation that the melting of solid nuclei of these sizes occurs via random breaches at the interfaces. Moreover, this breaching process occurs at the exact solid-liquid equilibrium temperature and in the absence of any external disturbance, which suggests the name “self-instability” for this melting process. We attribute this spontaneous instability to the curvature of the samples; based on the relationship between the sample’s instability and its curvature, we propose a destabilizing model for small systems. This model fits well with experimental results and leads to new insights into the instability behavior of small-sized systems; these insights have broad implications for research topics ranging from dendrite self-fragmentation to nanoparticle instability.

To make sure the length of nuclei long enough, they are selected as 20r 0 .

Supplementary Theories
The critical nuclei radius of cylindrical nuclei. According to the classical nucleation theory (CNT) 1~2 , to form a small solid sphere of radius r in an undercooled liquid, the change of the Gibbs free energy can be expressed as: where γ is the solid-liquid interfacial energy and ΔG V is the Gibbs free energy difference per unit volume between solid and liquid phases at the same temperature. In CNT, the approximated ΔG V can be expressed as: where L V is the latent heat of fusion per unit volume at the equilibrium meting point, T m is the equilibrium melting point, and ΔT(=T m -T) is the undercooling. The critical nucleus radius is obtained from equation (1): Considering an infinite long cylindrical nucleus, the difference in the Gibbs free energy and critical nuclei radius are amended as: For a given undercooling, there exists a critical radius r 0 , as shown in Supplementary Fig. S1.
If r < r 0 , the nucleus shrinks in order to lower the free energy of the system; if r > r 0 , it grows to lower the free energy; if r = r 0 , the nucleus can be in equilibrium (but metastable) with its surrounding liquid.
In our work, the cylindrical nuclei of Al from r=5 Å to r=50 Å were calculated using MD simulations. The potential in the simulations is Al1 3 , a kind of EAM.fs potential from Mendelev, that sets the melting point 3

The kinetic expression of interface fluctuation
A. Infinite flat interface. According to equation (5), ΔT is 0 for sample of infinite r 0 (flat interface), and the system keeps stable under this condition. Then if T < T m , the interface will move toward the liquid phase with a velocity V 5~9 : where Q is the thermodynamic driving force, defined as the difference of free energy between solid and liquid phases per atom, V 0 is a temperature-dependent factor representing the maximum velocity, and k B is the Boltzmann constant. For a flat interface, Q only includes the change of volume Gibbs free energy, which can be described as: where ρ is density of atoms.
In small undercooling limit, equation (6) can be linearized and yields: where the term multiplying the undercooling ΔT can be identified as the kinetic coefficient μ.
Under zero undercooling, the interface is stable, but the temperature fluctuation η introduces local movement of the interface with a velocity of μη. Therefore, the kinetic 10 of the interface can be expressed as: where h(x,t) is the profile of the interface, x is the direction along the interface, and Γ is related to the interfacial energy as follows: where γ+γ " is called interfacial stiffness. As a result of equation (9), concave and convex will form at different positions on this overall stable interface.
B. Finite size effect on interface dynamics. In the finite case, free energy Q not only includes the volume free energy but also the interfacial energy. For instance, for the growth of radius from r to r+dr, the difference in free energy is: where r is the radius of an infinite long cylindrical nucleus. Therefore the velocity expression in equation (8) is amended as: Therefore, the nucleus r = r 0 is stable under certain ΔT. The interface will move as a whole once r ≠ r 0 with a velocity V. As a consequence of equation (12), the kinetic expression of equation (9) should be amended as: where interfacial profile is described by r instead of h, which is the radius of the nuclei.
What must be pointed out that the radius r in equation (13) is treated just as a function of time t and the distance along the centerline of the nucleus x. However, it is also depends on θ, the angle around the centerline of the nucleus, which will be discussed in the later section (see Temperature profile of the interface atoms. According to statistical physics, the kinetic energy of atoms follows the Maxwell distribution: where the kinetic energy has relation to temperature by E i =3k B T i /2, T is the temperature of atom i. In our work, MD simulations of nuclei were simulated to obtain the statistic data of kinetic energy of interfacial atoms, defined by Morris 11 (see Supplemental Methods). In Supplementary Fig. S3 If there are n atoms exist on a finite element, Δx, along the centerline of the nucleus. The mean kinetic energy of these atoms follows the Gauss distribution according to the central limit theorem: Therefore, the temperature of this element Δx should be: where T-T x is thermal fluctuation ( η in equation (13) ) at x, represents the temperature difference between the local atoms and the system. The number of interfacial atoms, n, can be obtained by: where ρ is the number density of interfacial atoms. MD simulations of nuclei r=30 Å, 35 Å, 40 Å were used to obtained the value of ρ, were 0.16020 Å -2 , 0.15558 Å -2 , 0.15364 Å -2 respectively. In FDM simulations, a value of ρ=0.15 Å -2 was used for simplicity. Supplementary Fig. S4 shows the interface temperature in MD simulation of nuclei r=40 Å, where 18600 interfacial finite elements, Δx=1 Å, were counted (the black dots), which was consisted with the Gauss distribution (the red line).
It should be noted that, the thermal fluctuation also exist along the circumference. In this case, a finite element, ΔxΔθ, was considered. Here Δθ is a difference element along the circumference. Similarly, as the derivation above, the temperature of this finite element follows the Gauss distribution, and the equation (17) Therefore, the atoms whose parameter fall into [0.2, 0.4] were defined as interfacial atoms 11 . In situations of the shrinkage and growth, the interface moved simultaneously, the radius r(t) of the nuclei was calculated by averaging the radius of interfacial atoms.
In the research of the instability, the radius varied greatly on the breach of the nuclei. On this occasion, the interfaces were divided into series finite element, Δx, along the centerline of the nuclei, and the interfacial atoms belong to every finite element were picked out. Then radii r(x,t) were determined as the average of the interfacial atoms' radius. We consider this method to obtain the interfacial profile is rational, the error analysis was provided in Supplemental Discussion.

Profile of the perturbed interface.
According to the equipartition theorem 12 and capillary fluctuation method 13 , the relation between the interface profile and interfacial stiffness is: substitute r by h (h=r 0 -r), we get: The Gauss solution of equation (29) is: where C 0 and C 1 are constants only related to the initial interface morphology.

Dynamics of the interface.
To clarify the effect of the thermal fluctuation on the interfacial instability, FDM simulations were performed for nuclei r 0 =100 Å under four sets of conditions. The parameters used in FDM are shown in Supplementary Table S2.

Supplementary Discussion
The error analysis of the interfacial profile. In the kinetic expression of finite interface, equation (13), the dependence of the radius on θ, angle along circumference was neglected.
Considered this factor, the kinetic expression is amended as Because r(x,θ,t) is determined by η(x,θ,t), it follows to Gauss distribution, too. In simulations, the nuclei with r = 28 Å, 30 Å, 32 Å were relaxed at 892.975 K for 150 ps, the interfacial atoms were picked out for statistics. As showed in Supplementary Fig. S6, the radii tend to be Gauss distribution. From statistics, the 95% confidence interval of r vs. t was calculated, as shows in Supplementary Fig. S7. For the nucleus r = 30 Å, most radii fell into [25 Å, 35 Å], appeared stable in 150 ps; for the nucleus r = 28 Å, who shrank to 20 Å in 150ps, the error of the radii away from the average was no more than 5 Å; the same situation for the nucleus r = 32 Å, who shrank to about 40 Å in 150ps, most interfacial atoms weren't far from their average. Therefore, the method to obtain the radius r is reasonable.
In the research of the instability, the nuclei broke into two parts, gradually. The radius r cannot treated as uniform along the centerline, therefore, expressed as r(x,t). Then the thermal fluctuation along circumference led to deviation of r vs. θ. However, the distribution of interfacial atoms along circumference still followed to Gauss distribution. Supplementary Fig. S8 shows the interfaces' evolution of the nuclei r = 40 Å at 901.9 K. Three nuclei with different initial interface profiles were relaxed until to rupture. The error bar of 95 % confidence intervals of r(x,t) were given, the rationality of the method to extract r(x,t) was proved again.
The role of fluctuation in the process of instability. The processes of the instability are shown in Supplementary Fig. S9. From Fig. S9a and S9b, the breach emerged approximately at 5 ns for the straight interface with thermal fluctuation, the position of the rupture appeared randomly. For the perturbed interface in Fig. S9c and S9d, the original concave deviation on the interface played the role of the breach, leading to the rupture of the nuclei in 3.12 ns, no apparent difference was observed in spite of thermal fluctuation was introduced. We thus concluded that the thermal fluctuation was important for the formation of the breach and for the determination of the position of the rupture, but did not play a role after the breach formed.
The rapture time of interface. In previous section, the morphology of the rupture was obtained to be an approximate Gauss solution. The evolution of the breach tip is shown in Supplementary   Fig. S10, it kept bigger than 70 Å even at 2.5 ns. After that the nucleus broke into two parts promptly in 3.12 ns. When breaking starts, the velocity of the tip can be written as follows: C 0 /� 1 is the initial depth of the breaches on the perturbed interfaces, which is only related to the interfacial energy but not to the radius of the nuclei. The depth is around 4 Å for nuclei of r=10 Å to r=10 μm, see Supplementary Fig. S11. In this case, the term C 0 /� 1 in equation (35) can be treated as a constant.
With all the material's parameters were brought in, equation (35) can be obtained as: The result in the fitting is very close to the estimation of equation (36), which can be used to predict the rupture time of nuclei with different radius.