Nanowire melting modes during the solid–liquid phase transition: theory and molecular dynamics simulations

Molecular dynamics simulations have shown that after initial surface melting, nanowires can melt via two mechanisms: an interface front moves towards the wire centre; the growth of instabilities at the interface can cause the solid to pinch-off and breakup. By perturbing a capillary fluctuation model describing the interface kinetics, we show when each mechanism is preferred and compare the results to molecular dynamics simulation. A Plateau-Rayleigh-type of instability is found and suggests longer nanowires will melt via an instability mechanism, whereas in shorter nanowires the melting front will move closer to the centre before the solid pinch-off can initiate. Simulations support this theory; preferred modes that destabilise the interface are proportional to the wire length, with longer nanowires preferring to pinch-off and melt; shorter wires have a more stable interface close to their melting temperature, and prefer to melt via an interface front that moves towards the wire centre.

Nanostructured objects have lower stability with respect to their molten phase due to large surface area to volume ratios [1][2][3] . In the case of nanowires, their stability has been studied at elevated temperatures both experimentally [4][5][6] and theoretically 7,8 indicating the presence of Plateau-Rayleigh (PR) type of instabilities which can cause a nanowire to neck and breakup into a chain of nanospheres. In fact, PR like instabilities have been used as a means of self-assembly of chains of nanospheres for several different initial geometries ranging from rings 9 , wires 10 , and thin films 11,12 . PR theory generally predicts that the wavelength of the perturbations which cause a liquid wire to become unstable are proportional to the initial wire circumference (i.e. wires become unstable when c > 2πR 0 ). Moreover, linear stability analysis predicts a preferred wavelength that will drive a liquid wire to breakup. Much work has been done in regards to nanocluster stability during solid-liquid coexistence 13 , and the stability of liquid nanojets 14,15 , nanocylinders 16 , and even in cylindrical metal alloys 17 , relatively few studies address the stability of the solid core of nanowires close to their melting point.
It was found that for finite-sized cylinders during phase coexistence, differences in curvature and fluctuations would lead to the formation of random breaches at the material interface, causing the growth of instabilities which lead to the melting of the solid 18 . For finite-sized boxes, different crystal geometries are realised by overcoming nucleation barriers, where a crystal nucleus surrounded by its own fluid could change from a slab geometry to a cylinder, and then to a solid droplet. It suggests that the solid prefers metastable forms as the box approaches the freezing (or melting) density 19 . This indicates that stable geometries of a given phase depend on the volume fraction of said phase (or medium) to the system volume [20][21][22] . Recent work has studied the breakage of gold nanofilaments connecting two nanoparticles where the filaments connecting the two nanoparticles would break apart by Joule heating 23 . Moreover, it was observed that the temperature at the breakage point had a strong dependence on the filament width, and had a dependence on the length in some, but not all cases 23 .
The thermally induced breaking of nanowires becomes important when considering the role they play in devices that utilise nanowire networks. Heat can be generated in nanowire networks via current passing through the network, and as such can influence the morphology and breakup of the nanowires making up the network 24,25 . This could be a hindrance to device stability, where it is important to understand the limitations of interconnecting materials like nanowires.
In this paper, we investigate the stability of metal nanowires as they approach their melting temperature for copper nanowires of varying lengths and radii. To describe the nanowire stability, we perturb a capillary fluctuation model that describes the kinetics of the solid-liquid interface. The model is then tested against molecular dynamics (MD) simulations, where it is found that longer nanowires are more unstable with respect to the melt.

Results
Capillary fluctuation model. Melting at the nanoscale is thought to initiate at the surface, and move from the outside inwards, with the interface consuming the solid as it melts. However, observations in nanowires show the solid will begin to neck and breakup as the nanowire approaches its melting point T m 8 . Figure 1a shows a topdown view of a nanowire at a temperature T that sits between its surface melting temperature T s and bulk melting temperature T m (where T m is the bulk melting temperature of finite-sized nanowire). Figure 1b shows that as T → T m the solid is consumed radially as the interface moves towards the wire centre. Figure. 1c, d shows the melting instability mechanism. Figure 1c shows a side-on view of a nanowire where T s < T < T m . However, as T → T m , rather than the interface moving towards the centre, a portion of the solid begins to thin out and neck, as seen in Fig. 1d, initiating the solid breakup and causing the remaining solid to be consumed.
We first consider the Gibbs free energy difference per unit length in an infinitely long cylinder surrounded by its own melt close to T m (see 18

supplementary theories)
The value of r that minimizes this equation represents the equilibrium solid radius for an infinitely long nanowire in the vicinity of T m Here, γ sl represents the solid-liquid interfacial energy, L v is the bulk latent heat of melting per unit volume, T c melting temperature of the bulk materials and T is the undercooling. Now we look at the interface velocity for a cylindrical nucleus 18 . For an infinite flat interface there is zero undercooling, and if T < T c , the solid-liquid interface will propagate towards the liquid phase with a velocity V [26][27][28][29] where V 0 represents a maximum velocity that depends on temperature, Q is defined to be a thermodynamic driving force, and k b is Boltzmann's constant. Q is defined as the difference between the solid and liquid phases per atom, so in a flat interface limit, it can be approximated as Q ≃ L v �T/N T c , where N is the number density. Taking Eq. (3), substituting for Q and taking T = T c − T , and expanding in the small undercooling limit the interface velocity can be linearised as where ζ is a kinetic coefficient. This gives the planar interface velocity in the small undercooling limit. We now consider the interface kinetics by looking at the dynamic behaviour of a cylindrical interface with a profile r(z, t) 18,30,31 where Ŵ = (γ sl + γ ′′ sl )T c /L v , with γ sl + γ ′′ sl the interfacial stiffness of the solid-liquid interface 32 . Assuming an isotropic solid-liquid interface and small anisotropy, this can be approximated as where C is a constant, and the delta functions suggest the noise is uncorrelated in space and time. We perturb the solid-liquid interface by a small parameter ε , and express it as the surface r(z, t) = r * + εe ikz+ωt , with k and ω being the wavenumber and instability growth rate respectively (see in Fig. 2). The term (r * + εe ikz+ωt ) −1 can be approximated as By substituting r into Eq. (5), using the expression in Eq. (6), solving to O(ε) , and using the definition of r * in Eq.
(2) we recover A PR type instability can be found by observing that ω > 0 when kr * < 1 , bringing us to the familiar solution Combining Eqs. (8) and (2) and taking T = T c − T m , the interface will remain stable when ω < 0 and leads to the moving interface front seen in Fig. 1a to b giving the condition If ω > 0 then T ω > T m , and perturbations will grow to destabilise the solid-liquid interface, giving the scenario in Fig. 1c-d. If * ∝ L then T ω becomes larger than T m quickly, giving a criteria describing when each melting mode is preferred.
Finally, we look at the equation for the bulk melting temperature of a cylindrically symmetric nanowire of radius R 0 from previous work 8 , where κ = 1−�γ /γ sl 1+�γ /γ sl ( �γ is the spreading parameter that determines the wettability of a material 8 ), ξ represents the correlation length, I 0 and I 1 represent modified Bessel functions of the first kind respectively. Equation (10) is found by solving a two-parabola Landau-type model for the free energy of a cylindrically symmetric nanowire 8 . Combining Eqs. (2) and (10), a equation for r * in terms of R 0 and interfacial energies κ can be found Methods MD simulations were performed using LAMMPS 33 using an embedded-atom-model (EAM) potential for copper 34 , with a bulk melting temperature of 1320 K. This potential was chosen because yielded good results for melting temperatures and dynamics 8 . Periodic boundary conditions in all directions were used, with the periodicity along the wire axis effectively simulating an infinitely long wire, and additionally suppress long-wavelength instabilities that may otherwise cause the wire to break apart prior to completely melting. As such, the maximum wavelength allowed by the system will be equal to the box size along the axis of the nanowire.
The equations of motion were integrated using a Verlet method using a timestep of 2.5 fs. The temperature was controlled with a Langevin thermostat, which effectively simulates Brownian motion, and is a appropriate choice of thermostat since To control the temperature a Langevin thermostat with a damping parameter of 1.0 ps −1 . This was to ensure a quick equilibration at each timestep of the simulation. The simulations were initialized at an initial temperature T i for 1.0 ns. Afterwards, a production phase for each wire was run from a temperature T i to a temperature T i + 1 . Then an equilibration phase around the temperature T = T i + 1 was run. Each production phase was 0.40 ns, with an equilibration phase of 0.60 ns, creating an effective heating rate of around 1K/ ns. This ensured us that at each temperature the wires were sufficiently close to equilibrium. See supplementary materials S2 and S3 for more information on computational details and methods used in this study. We can study the stability of the solid-liquid interface by examining when the solid core begins to pinch-off for two wires of the same radii, but different lengths, as shown in Fig. 3. As T → T m , the interface will either move towards the wire centre Fig. 3a or begin to pinch-off Fig. 3b. In Fig. 3a the size of the liquid nucleus is large compared to the much longer wire. The presence of solid atoms close to the liquid surface in Fig. 3b can be seen by the presence of a 'noisy' interface (red solid line). We will see in this section, wires with lengths that satisfy L > 2πR 0 will pinch-off and melt at a temperature consistently lower than wires with lengths L ≤ 2πR 0 .
To study the stability of the solid-liquid interface, the Fourier transform of the solid is taken to extract modes that destabilise the interface (see supplementary details S3). Figure 4 represents a stability diagram in terms of the fastest growing modes k sol r * against wire aspect ratio L/R 0 , where each k sol r * is the averaged value of kr * that represents the maximum Fourier transform amplitude (see supplementary details Fig. S3). The modes k sol r * for each wire aspect ratio are similar, which indicates destabilising modes depend strongly on the wire aspect ratio. As nanowires get shorter, modes that destabilise the interface approach unity, in violation of classical PR theory (see red dashed line). Also seen are two regimes that identify the preferred melting mode, as seen in Fig. 1. To the left (light pinkish region) we see the regime where T ω < T m which indicates the solid-liquid interface must move closer to the centre before the pinch-off can initiate. On the right (dark bluish region) we see the regime where T ω > T m and identify when the pinch-off melting mechanism is favoured. Observations from MD simulation agree with the theory developed, represented by Eqs. (7), (8) and (9). Longer wires will be more thermodynamically unstable since wire lengths will generally be greater than the circumference of the coexisting solid, as indicated by Eq. (8). Included in Fig. 4 is the scaling relation k sol r * ∝ 2πR 0 L , which follows the trend observed in MD simulations, as well as predictions by classical PR theory which states k max R 0 ≃ 0.697 . For wires with L < 2πR 0 , k sol r * approaches and exceeds unity, as predicted by the scaling relation, violating classical PR theory and indicating regions of interface stability. Moreover, we see that for higher wire aspect ratios, PR theory overpredicts the fastest growing modes in the solid. This has been previously reported when examining the stability of liquid nanojets, implying that at small length scales classical PR theory is not wholly sufficient to predict interface stability 16 . This too is the case in liquid nanowires (see supplementary details). www.nature.com/scientificreports/ We now examine how the wire length influences the stability of the solid-liquid interface by looking at how r * obtained from simulation behaves close to the melting point for wires of different radii and lengths. From Eq. (2) we can see that r * ∝ T −1 . As observed in Fig. 5, MD results consistently show that longer wires melt at a lower temperature since they are prone to the growth of instabilities that initiate the pinch-off. Shorter wires not only melt at a slightly higher temperature, but the value of r * at the point when the pinch-off initiates is consistently smaller, in agreement with the theory developed. It supports the idea that there are two melting mechanisms that depend on the wire aspect ratio, as evident in Fig. 4. The simulated results point to the bulk melting temperature of the potential used being between 1220 and 1230 K, rather than the 1320 K stated.
In Fig. 6 we see that r * depends not only on the initial wire radii R 0 , but also on the wire aspect ratio as well. This evidence supports the theory and previous claims that as the wire aspect ratio gets smaller, the solid core radius prior to the pinch-off decreases. The ratio of r * /R 0 appears to converge in the limit of large L/R 0 . Each Figure 4. This figure shows how the wire aspect ratio affects the fastest growing modes. A thick red line plots 2π/(L/R 0 ) , which we assume is proportional to k sol r * . The light shaded area shows the region where the interface is stable at temperatures close to T m (radial mode), whereas the darker shaded region shows where the interface is expected to be unstable (instability mode).

Figure 5.
Values of r * obtained from simulation against T −1 . Yellow circles are for r * when L ≤ 2πR 0 , and red diamond points represent r * and T for when L > 2πR 0 . Yellow circles indicate when the radial mode is the preferred melting mechanism, and red diamonds indicate where the instability mode is preferred. The first two points (bottom left) are for wires where R 0 ≈ 22 Å the middle two are for R 0 ≈ 30 Å, and the last two are for R 0 ≈ 38 Å. We assume yellow and red points represent the bulk and instability melting temperatures, T m and T ω , respectively. www.nature.com/scientificreports/ r * /R 0 for the aspect ratios studied are similar when L > 2πR 0 . Once L ≤ 2πR 0 the difference in r * /R 0 becomes appreciable. This could be indicative that quantities like the interfacial energies and correlation length ( γ sl and ξ respectively) become important quantities for small wires, implying size and curvature play key roles in observations for small aspect ratios. (See supplementary details for the mean and standard deviations for T m , r * and k sol r * ).

Discussion
By perturbing the interface of a surface melted metal nanowire, we can describe the existence of two mechanisms for nanowire melting; an instability mode or radial mode. The model showed the fastest growing modes that destabilise the interface are inversely proportional to the wire length, where a PR type instability for the solid in a surface melted nanowire is recovered. By using classical nucleation theory and exploring the nanowire stability in the vicinity of T m , we were able to define the condition that determines the preferred melting mechanism. Moreover, we recovered an expression for the equilibrium solid radius in terms of the initial wire radius and the interplay of the interfacial energies of the nanowire. Simulations show the fastest growing modes are inversely proportional to the wire length, and in fact that k sol r * ∝ 2πR 0 L . Additionally, we observe longer nanowires consistently melt at a lower temperature than shorter wires, in agreement with our developed theory and other recent observations 23 . The implication is that shorter nanowires have a more stable interface when close to their melting temperature. For nanowires where the instability mechanism is the preferred melting mode, once the pinch-off has initiated the remaining solid will be consumed. This is because the solid core tries to stabilise itself by forming into a sphere, minimising its surface energy. In some cases it was observed that for the longest, thinnest nanowires, the liquid-vapour interface would begin to neck, being driven by surface diffusion, which in turn influenced the breakup of the solid. For longer heating rates or overdamped Brownian dynamics, this feature would become more pronounced. However, due to the quick equilibration at each timestep, this was not an issue.
Evidence can also be seen that r * depends on the wire aspect ratio and not just the initial radius. This has been reported in previous work, where for nickel and aluminium nanowires of a single length but increasing radii, the solid core remained stable down to smaller radii 8 . Studies have explored the size-dependence of interfacial energies 35 . Our study, however, shows surface area of metal nanowires becomes an important factor in interfacial energies for small wire lengths. Curvature too plays a role in the interfacial energy, where for spherical clusters the solid-liquid interface energy is linear with inverse radius 36,37 . This size and curvature dependence explains why values of r * /R 0 begin to deviate away at low aspect ratios. The ratio of atoms at the surface compared to the bulk becomes far more appreciable for the smallest wires, giving the curvature a greater role in the solid-liquid interface dynamics. Given the fact that the fastest growing modes are inversely proportional to the nanowire length, it would be of no surprise that interfacial energies will depend on this too since their surface area will scale with radius and length. The theory and simulations show that long nanowires are thermodynamically unstable at high temperatures since the nanowire length will almost always be much greater than its equilibrium solid radius. This has ramifications when considering device stability that utilise nanowires subjected to heating. We observed that for long, thin nanowires, the liquid-vapour interface can begin to destabilise even before the solid begins to neck. This implies ultra-long, thin nanowires will be particularly unstable at elevated temperatures and should be considered when constructing nanowire devices. Figure 6. The ratio of equilibrium solid radius to initial wire radius against the wire aspect ratio. Values of r * /R 0 for each aspect ratio are clustered closely together, where they begin to deviate markedly when the initial wire length L < 2πR 0 .

Conclusion
We studied the stability of the solid in copper nanowires as they approach their melting temperature by perturbing a model describing interface kinetics and compared the results to MD simulations. The model found a stability criterion that dictates the preferred melting mode a nanowire will take. We found that longer nanowires are thermodynamically unstable, and will preferentially pinch-off and melt, indicating a melting mechanism driven by a PR type of instability. In shorter nanowires, the interface front moved radially towards the nanowire centre before the solid would breakup, indicating higher interface stability, with MD results in agreement with our model. Moreover, we proposed modes that destabilise the solid-liquid interface are dominated by the nanowire length, in contrast to PR theory which states they are proportional to the circumference. Additionally, it was observed from the MD simulations that longer nanowires consistently have a melting temperature a few degrees below shorter nanowires, indicating the nanowire aspect ratio influences the preferred melting mode and solid-liquid interfacial stability.

Data availability
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.