Including nonequilibrium interface kinetics in a continuum model for melting nanoscaled particles.

The melting temperature of a nanoscaled particle is known to decrease as the curvature of the solid-melt interface increases. This relationship is most often modelled by a Gibbs-Thomson law, with the decrease in melting temperature proposed to be a product of the curvature of the solid-melt interface and the surface tension. Such a law must break down for sufficiently small particles, since the curvature becomes singular in the limit that the particle radius vanishes. Furthermore, the use of this law as a boundary condition for a Stefan-type continuum model is problematic because it leads to a physically unrealistic form of mathematical blow-up at a finite particle radius. By numerical simulation, we show that the inclusion of nonequilibrium interface kinetics in the Gibbs-Thomson law regularises the continuum model, so that the mathematical blow up is suppressed. As a result, the solution continues until complete melting, and the corresponding melting temperature remains finite for all time. The results of the adjusted model are consistent with experimental findings of abrupt melting of nanoscaled particles. This small-particle regime appears to be closely related to the problem of melting a superheated particle.

The size-dependent nature of melting nanoscaled metal particles has received attention from a number of experimental studies123. A well accepted model for the melting temperature is given by the Gibbs–Thomson formula where R* is the radius of the spherical particle and is the bulk melting temperature, which is the temperature at which the material would melt if the interface was flat. The behaviour has been verified experimentally for a number of different metals including silver1, gold4, lead5, tin67 and aluminium8 (see Mei and Lu9 for a more complete summary of experimental results). The physical constants ρ, L and σ* are, respectively, the density of the material in the solid phase, the latent heat of fusion, and a parameter proportional to surface energy effects acting on the solid-melt interface.

The effect of equation (1) is that for small R*, the melting temperature of the spherical particle is significantly reduced. This size dependence on melting temperature is a consequence of nanoparticles having a much larger surface-to-volume ratio than bulk materials, and occurs for both for round and facetted particles. Since both the solid and liquid molecules on a curved surface are more weakly bonded than their counterparts in the solid and liquid bulk, the difference between the binding of liquid and solid molecules on the surface is a driving factor for this reduced melting temperature. Thus, ultimately the size dependence on melting temperature of a nanoscaled spherical particle is due to the very high surface-to-volume ratio and lower interfacial energy of the liquid phase. We note that there is no observable reduction in melting temperature for macrosized particles; this is a small-scale phenomenon only, as the length scale ω in equation (1) is typically of the order of nanometres (see data in Table 1 for various metals).

Table 1

Approximate thermodynamic constants used to calculate the nondimensional parameters (16) including references, except for and ρ, which are well known. The constant ω has been calculated with equations (1)–(3), except when values for and could not be found. The column ωG reproduces the values of Guisbiers et al.10 for an alternative model to equation (2)

	L × 103 (Jkg−1)	ρs × 103 (kgm−3)	ρℓ × 103 (kgm−3)	(Jm−2)	(Jm−2)	ω (nm)	ωG10 (nm)	cs (Jkg−1K−1)	cℓ (Jkg−1K−1)	kℓ/ks	(K)	β
Ag	10256	10.5	9.3	1.2157	0.9057	0.42	0.94	233	29258	0.4811	1234	0.35
Au	62.74	19.3	17.31	1.384	1.144	0.26	0.92	12912	1694	0.3312	1336	0.36
Pb	22.913	11.34	10.66	0.6113	0.4813	0.85	1.45	12812	14812	0.4612	600	0.30
Sn	58.93	7.18	6.98	0.667	0.557	0.47	1.02	22759			505	0.51
Al	3968	2.7	2.385		0.9160		1.28	89759	108058	0.4258	933	0.47
Cu	20556	8.96	8				0.82	38559	48058	0.8558	1358	0.39
Ti	29659	4.5	4.11				0.89	52459	70058	0.9158	1940	0.29

The constant σ* in equation (2) is a measure of the surface energy effects, also referred to as surface tension. One long-standing model for σ* is where and are the surface energies of the particle in the solid and liquid states, respectively1415. The Gibbs–Thomson relation (1)–(3) has been shown to be consistent with the experimental data for a variety of metals12416. To take one example, in Fig. 1 we show a comparison of this model against melting temperatures of gold particles, measured by Dick et al.2. Note that gold particles with a radius of 10 nm have a melting temperature of approximately 1250 K, which is noticeably lower than the bulk melting temperature ( for gold). For smaller particles we see the agreement in Fig. 1 is still quite good, even down to a radius of less than 1 nm.

Figure 1

The size dependent nature of the melting temperature of gold nanoparticles demonstrated by experiment data (○) from Dick et al.2.

The measurements are in good agreement with the Gibbs–Thomson relation equations (1)–(3), plotted here (dashed) with the thermodynamic constants for gold24.

There are several other models for σ* (including Kofman et al.5, Semenchenko17 and Wronski18) with some agreeing better with the experimental data than others13. These different models can lead to varying values of σ*. For example, equation (3) for lead gives σ* = 0.11 Jm−2 (using the values provided in Table 1), while other models1920 give σ* = 0.031 Jm−2 or 0.046 Jm−2 < σ* < 0.145 Jm−2. In fact, there is a whole series of theoretical models for ω in equation (1) (see the summary by Guisbiers21.), some of which do not follow equation (2). To take one example, there is a model due to Guisbiers, Wautelet and coworkers102223 which predicts higher values of ω, as demonstrated in Table 1. We revisit the uncertainty in these measures in the Discussion section.

In a continuum model for the melting of a free-standing nanoscaled particle, the Gibbs–Thomson rule (1) acts as a boundary condition on the moving solid-melt interface r* = R*(t*)12242526272829. This approach brings with it a number of issues that need addressing. First, as it is, the melting temperature (1) is singular as R* → 0, and indeed implies that for a sufficiently small radius (R* < ω), the melting temperature must drop below 0 K. Thus any use of equation (1) must be applied with the understanding that there is a limiting radius below which equation (1) can never hold. Second, one must question the reliability of a continuum model that involves classical notions of melting and thermodynamics for very small radii, especially when the particle size is of the order of a couple of nanometres. Lastly, the continuum model for melting a spherical particle is known to exhibit a form of mathematical blow-up, which predicts that the speed of the solid-melt interface increases without bound as the radius reaches a critical value 26. According to the model, the flux of heat into the interface also blows up in this limit. Mathematical solutions of this form may be interpreted as modelling a type of abrupt melting that has been observed when melting nanoscaled particles9. Here, when the particle's radius reaches a critical value, the melting process for the remaining core appears to occur instantly. This phenomenon is known to occur for particles of lead913 and gold30, for example. Despite this natural phenomenon, a more satisfactory output from a continuum model should be that solutions are mathematically well behaved for all time.

The Gibbs–Thomson effect (1) can be derived by considering the temperature at which a solid sphere of radius R* is in equilibrium with its melt3132, but this derivation requires the system to be in equilibrium; that is, a key assumption behind equation (1) is that the interface is not moving. In practice this implies that equation (1) is valid when the departure from equilibrium is small, while in reality, equation (1) should include a correction for nonequilibrium kinetic effects. The most simple form of this correction is a linear additive term so that the adjusted Gibbs–Thomson rule becomes where the length scale ω is given in equation (2) (or calculated using one of the many other models21) and we may refer to the parameter as the kinetic coefficient3233. Note that a more general correction term which is not linear in the interface speed is discussed in a number of studies33343536.

Given the material properties of a particular metal, we are unable to simply plot the melting temperature against the particle radius R* with equation (4) without information about the speed of the solid-melt interface dR*/dt*, which itself will have some complicated dependence on the size of the particle and other experimental conditions. For moderate interface speeds we may expect the kinetic parameter to be sufficiently small that the magnitude of the second term on the right-hand side of equation (4) is much smaller than the first, so that in practice the melting temperature will be well approximated by equation (1). However, as noted above, solutions to the continuum model with equation (1) have interface speeds that continue to grow without bound26, and on the experimental side there are observations of particularly high interface speeds for very small particles913. Thus, regardless of how small is, there is a small-particle regime in which the full condition (4) appears more appropriate than equation (1).

In this paper, we consider a continuum model for melting a sphere that imposes equation (4) as the relevant boundary condition on the solid-melt interface, and find that the resulting solutions exhibit features that shed light upon the issues listed above. By numerical simulation, we show that including a nonzero kinetic term ( in equation (4)) in the model acts to suppress the mathematical blow-up, so that the solution naturally continues for all time until the particle is completely melted. By treating the regime in which the nondimensional kinetic parameter is much less than one, we show that solutions are consistent with the notion of abrupt melting. Further, with numerical solutions to the full continuum model, we are able to plot equation (4) and show that, even with a small kinetic parameter, the melting temperature is well behaved for all particle radii; in particular, with we no longer have the unphysical prediction that as R* → 0. Thus when we include the correction term for nonequilibrium interface kinetics, we no longer need to cut off the model at an arbitrary particle radius.

The present study is similar in nature to the recent work of the authors24, which also considers the effect of interface kinetics on the melting of a sphere. Two minor differences are that the problem in Back et al.24 is formulated on a finite domain, which is slightly artificial from an experimental perspective, and the nondimensionalisation is not the same. However, the key difference is that the mathematical model in Back et al.24 is a simplified one-phase approximation that arises under the assumption that heat conduction occurs much more easily in the liquid phase than the solid phase. This assumption is not appropriate for melting nanoparticles (especially metals such as gold, leads, silver, etc.). While the results in that paper are instructive as they deal with regularising mathematical blow-up via a kinetic term, it is important to treat the full two-phase problem formulated without unrealistic simplifications.

In the Results section, we outline the continuum model for melting a spherical particle. This model includes the linear heat conduction equation in both the liquid and solid phases, and a dimensionless version of equation (4) imposed on the moving boundary r* = R*(t*). We then document numerical results for the cases in which the kinetic term is zero and nonzero. These results demonstrate how the inclusion of nonequilibrium kinetic effects in the Gibbs–Thomson rule acts to prevent the mathematical blow-up, mentioned above, from occurring. In the Discussion section, we summarise the implications of our results on the use of the Gibbs–Thomson rule and, more generally, in terms of mathematical modelling of nanoscaled particles. Finally, we discuss how abrupt melting is driven by a counter-intuitive heat flux from the particle into the solid-melt interface. We argue that this process is closely related to melting a superheated particle, even without the effects of surface tension or interface kinetics.

Results

Continuum model

We consider a single radially symmetric solid nanoparticle of radius , initially at a temperature . This particle is suspended in an infinite liquid with the temperature at infinity held at a constant for all time. The liquid has an initial temperature with the property as r* → ∞. If , the particle will melt. In this model, we assume that the inner solid core and the outer semi-infinite liquid region are separated by a distinct boundary, the solid-melt interface r* = R*(t*), which moves inwards towards the centre of the sphere as the particle melts. We seek to solve for the temperatures in the solid and liquid regions, and , respectively, as well as for the position of the interface R*(t*).

A dimensional continuum model for this idealised scenario is as follows. Assuming the heat diffuses through each phase via conduction only, we have where, despite the use of equation (3) to compute σ*, we now assume the density ρ is constant during the phase change process. On the solid-melt interface r* = R* we have the boundary condition referred to as the Stefan condition31, which models how the solid at the interface absorbs latent heat to phase change into liquid. This coupling of the temperature gradient on the interface and the speed of the moving front shall be important later when we discuss finite-time blow-up and abrupt melting. The thermodynamic constants in each phase are the thermal conductivity k,, specific heat c, and latent heat of fusion L. A second boundary condition on the moving interface is that the temperature is continuous between the two phases, and that the temperature here is equal to the melting temperature. That is, on r* = R*, where comes from the Gibbs–Thomson law (4). Lastly, we have a no-flux boundary condition at the centre of the sphere and the condition at infinity, as well as initial condition for the position of the moving front, which we write as at t* = 0.

It proves particularly insightful to scale the problem using the following dimensionless variables: With these new variables, equations (5)–(6) become subject to the Stefan condition and the remaining relevant initial and boundary conditions

Here, the dimensionless parameters are the ratio of thermal conductivity, the ratio of specific heat, and the Stefan number, which is a measure of the latent heat absorbed by the solid phase during melting. By nondimensionalising our continuum model according to equation (8), in which the representative length, time and temperature scales are unique for each material, the resulting dimensionless parameters in equation (16) can be determined using available data such as that listed in Table 1 and do not depend on the boundary or initial conditions. This property will be important later, when we discuss universal features of the model.

The final parameter in the model (9)–(13) is the dimensionless kinetic parameter A nonzero value of in the model is required to implement the adjusted Gibbs–Thomson rule (4), which accounts for nonequilibrium kinetic effects. While, in principle, values of the parameter may be determined for a given material, we shall argue in what follows that, regardless of its value, solutions to the continuum model that include nonequilibrium kinetic effects have a number of appealing features. A further important remark on the scalings (8) is that, since the standard Gibbs–Thomson law (1) can never hold for R* < ω, our dimensionless continuum model with cannot ever be used for the dimensionless regime R < 1 (see equation (12)). On the other hand, with , the boundary condition (12) does not necessarily lead to the same restriction. Indeed, provided the magnitude of dR/dt is large enough, mathematical solutions for can continue until the particle is completely melted.

The model governed by equations (9)–(15) is commonly known as a two-phase Stefan problem with surface tension and a kinetic term (referred to as ‘kinetic undercooling' in the Stefan problem literature). Without the effects of interface kinetics (setting but keeping σ > 0), variations of this spherically-symmetric melting problem have been studied mathematically by a number of authors using a variety of numerical and analytical techniques1225262728293738. Neglecting both the kinetic term and surface tension leads to the classical problem of melting a macroscaled sphere (see McCue et al.39 and the references therein). One-dimensional Stefan problems with kinetic undercooling () but without surface tension (σ = 0) have also received attention in the applied mathematics literature, both to model solidification or melting processes35364041 and in the context of diffusion through glassy polymers4243. In numerical studies of two- and three-dimensional models for instabilities and pattern formation in crystal growth phenomena, it is common to include both surface tension and nonequilibrium interface kinetics in the Gibbs–Thomson law444546. As far as we are aware, with the exception of the recent work of Myers and coworkers1247, there has been no in-depth discussion in the physics literature on the use of two-phase models like equations (9)–(15) for melting nanoscaled metal particles. While we acknowledge that metal particles are not always spherical (nanoscaled particles can be polyhedral in shape4849), our radially-symmetric model is significantly easier to handle than its three-dimensional counterpart (cf. the detailed asymptotic study of melting a three-dimensional body50), and we expect the key results will carry over.

Numerical results

The dimensionless problem (9)–(15) is solved computationally using the numerical scheme detailed in the methods section. To illustrate our numerical results, we have chosen to use the (dimensionless) parameter values β = 0.30, k = 0.46 and c = 1.16, which corresponds to melting a pure lead nanoparticle. Temperature profiles are shown in Fig. 2, while the dependence of the solid-melt interface on time is illustrated in Fig. 3.

Figure 2

Temperature profiles for β = 0.30, k = 0.46, c = 1.16 and Rinit = 23.65.

The initial temperatures are shown by the (blue) dot-dashed lines. The lead particle is initially at T,init = 0.9, while the surrounding liquid melt is at T,init = T,∞ = 1.25, except near Rinit where the initial condition is chosen such that the temperature is continuous, as discussed in the methods section. The (green) solid lines are for while the (black) dashed lines are for , and the (red) thin solid lines and the (red) thin dashed lines are the melting temperatures calculated from equation (4) for each case. (top) Temperature profiles for Regime I, 0 < t < tI = 103.34, are shown (from right to left) for times corresponding to R = 22 and R(tI) = 18.26. The inset gives a magnified view of the temperature near the moving boundary, where ∂T/∂r → 0+ as . (bottom) Temperature profiles for Regime II, tI < t < tc, are shown (right to left) for times equivalent to R(tI), R = 8, 4, 2 and the critical radius Rc = 1.32. Finite-time blow-up occurs for the case at tc = 400.98. From the inset we see that finite-time blow-up is avoided for the case at the time tc, as the flux ∂T/∂r is finite.

Figure 3

The radius of the particle R(t) versus time t for the same parameter values as in Fig. 2.

The (green) solid line is for the case , while the (black) dashed line is . For the case, we have that the speed of the boundary blows up at Rc = 1.32, corresponding to tc = 400.98. The solution with a nonzero kinetic term follows that for very closely, except for times near the blow-up (inset). Here, the solution deviates away from the case, so that the moving boundary propagates inwards past the critical radius and through the blow-up regime. The speed of the boundary increases dramatically until the boundary reaches the centre, and complete melting is achieved.

As we argue below, there are universal features of the model (9)–(15) that are independent of the initial conditions and far-field condition. With this in mind, we note that for the results shown in Figs. 2 and 3, we have chosen to use T,init = 0.9 and T,init = T,∞ = 1.25, which correspond to the dimensional temperatures and . Further, we use Rinit = 23.65, which is chosen to correspond to if we use the length scale ω from equations (2)–(3). However, as discussed in the Introduction, using other models for ω would provide different dimensional lengths (for example, the model of Guisbiers and Wautelet102223 would lead to the larger initial radius ).

No kinetic term, . Before considering the effects of interface kinetics, it proves useful to summarise the results found by setting in equation (12) (see also McCue et al.26). These are represented in Figs. 2 and 3 by the (green) solid curves.

The solutions for are well described qualitatively by dividing the dimensionless time domain into two regimes; 0 < t < tI and tI < t < tc. Here tI is defined to be the time at which heat ceases to flow from the solid-melt interface into the solid, while tc is the critical time at which finite-time blow-up occurs. Regime I (0 < t < tI) can be thought of as ‘normal melting': heat flows from the outer liquid phase R(t) < r < ∞ into the solid-melt interface r = R(t), providing enough latent energy there to transform the solid into liquid. This behaviour is demonstrated by schematic (a) in Fig. 2 (top). The surplus of heat energy then flows from the solid-melt interface into the solid 0 < r < R(t), raising the temperature of the solid. Temperature profiles for Regime I are shown in Fig. 2 (top).

At early times, the radius of the particle is large enough such that the Gibbs–Thomson effect is relatively small and the size-dependent melting temperature is approximately equal to the bulk melting temperature (which is consistent with the experimental data in Fig. 1). As time increases and the particle size decreases, the Gibbs–Thomson effect lowers the melting temperature on the moving boundary, which in turn decreases the temperature gradient in the inner solid phase at the interface. This trend continues until ∂T/∂r(R, t) = 0 at the time we define to be tI. Thus, at the precise moment t = tI, there is no flux of heat from the interface into the solid (see the inset in Fig. 2 (top), as well as schematic (b)).

In Regime II (tI < t < tc), we have that the temperature gradient ∂T/∂r(R, t) < 0 (see Fig. 2 (bottom) and schematic (c)), which means that now heat is flowing from the solid phase into the solid-melt interface and that the temperature in the solid phase in the neighbourhood of the interface is higher than the melting temperature itself. In this sense the solid can be thought of as being locally superheated (although the temperature here is still less than the bulk melting temperature). As heat is still flowing into the interface from the liquid phase, the melting process accelerates quickly, leading to a dramatic increase in interface speed |dR/dt| via equation (11) (see Fig. 3). The reduced melting temperature results in an ever decreasing temperature gradient ∂T/∂r(R, t) < 0, until the solution blows up at the critical time tc. This form of mathematical blow-up is characterised by ∂T/∂r(R, t) → −∞, dR/dt → −∞ and as , where Rc > 0 is the critical radius (see schematic (d)). As a consequence, according to the continuum model (9)–(15) with , the particle does not melt completely.

An interesting and insightful feature of this mathematical model (for ) is that the behaviour of the solutions near blow-up appears to be largely independent of the initial and boundary conditions. As a means to test this idea, for a fixed set of parameter values for c, k and β, we ran many simulations for a variety of different initial temperature profiles T,init and T,init, including constant values and spatially-dependent functions. Further, we have used a number of different values for the initial radius Rinit. Performing this exercise leads to vastly different temperature profiles for sufficiently small times, as would be expected, since small time behaviour is heavily dependent on initial conditions. However, for times close to the blow-up time tc, we found that regardless of the initial or boundary conditions, the temperature profiles appear almost identical. Thus, for example, the temperature profiles for Fig. 2 (top) would not be replicated if the initial or boundary conditions were different. On the other hand, we found that the profiles in Fig. 2 (bottom), especially those near the blow-up time t = tc, were essentially reproduced (with some very small quantitative differences) for a large number of initial and boundary conditions. This apparently universal near blow-up behaviour for was not observed by McCue et al.26 because the problem there was not scaled according to equation (8) (or some equivalent scaling that involves only material parameters).

Nonzero interface kinetics, . We now consider what happens to our model when we include a kinetic term, so that both σ > 0 and in equation (12). To illustrate the main results, we choose the specific value , together with our parameter values for lead. This value of is not taken from experimental data, but is used to demonstrate representative behaviour for small values.

Figure 2 (top) shows that in Regime I, the effect of turning on a kinetic term with a small kinetic parameter is negligible, as the temperature profiles for ((black) dashed lines) are virtually indistinguishable from the solutions described in the subsection “No kinetic term, ”. As such, the melting process appears to enter Regime II at essentially the same particle radius.

For much of Regime II the effects of the kinetic term are still small. In Fig. 2 (bottom) we see the temperature profiles for when follow those for the case very closely. However, close to the blow-up regime, the speed of the interface becomes sufficiently large that the kinetic term changes the qualitative behaviour. The melting temperature (12) for is increased, and the blow-up is avoided. In particular, we see in the inset of Fig. 2 (bottom) that the profiles for have a finite slope, so that the flux of heat from the solid into the solid-melt interface does not increase without bound. While not shown in this figure, the solutions for continue and remain regular. These features are also visible in Fig. 3, where the inclusion of interface kinetics in the melting temperature (12) has prevented the interface speed from blowing up at a finite particle radius.

Now that finite-time blow-up is avoided for , the moving boundary r = R(t) continues towards the centre into a third extremely short regime tc < t < te, where te is the extinction time characterised by R(te) = 0. As we see in Fig. 3, the interface speed during this regime is large but finite. In summary, we see that the addition of a small kinetic parameter has regularised the singular behaviour, so that the solution continues through the blow-up regime until complete melting occurs.

Discussion

This work concerns the use of a standard Stefan-type continuum model for melting a sphere, but with a non-classical size-dependent melting temperature that is motivated by experiments on melting nanoscaled metal particles. This form for the melting temperature is often taken to be the Gibbs–Thomson law (1). One obvious limitation of equation (1) is that it is unphysical for R* < ω (predicting negative temperatures), thus results from any mathematical model with equation (1) must be invalid for some small-particle regime. Further, the continuum model itself will lose its validity for very small particles, as we discuss below.

It turns out that the two-phase Stefan problem with equation (1) undergoes a form of mathematical blow-up before melting is completed, so the issue of cutting off the model at some arbitrary R* does not necessarily arise. This form of blow-up is interesting as it appears to be universal in the sense that the behaviour of the mathematical solution at times just before blow-up is essentially independent of boundary or initial conditions. In addition, this blow-up behaviour involves the speed of the solid-melt interface increasing without bound, which appears to be consistent with experimental observations13 and results from molecular-dynamics simulations30 of nanoscaled metal particles melting abruptly9.

In order to appreciate the physics behind this form of blow-up, it is worth comparing with the idealised problem of melting a superheated solid. Using the language of the subsection “No kinetic term, ”, we see that once the melting process has entered Regime II, there is a counter-intuitive heat flux from the bulk of the particle into the solid-melt interface, which is acting like a heat sink. This is precisely what occurs when one melts a superheated solid, even on a macroscale without the effects of surface tension. A very simple model for melting a superheated solid is41515253 For this well-studied problem one can consider the quantities

The term Q1 describes the amount of dimensionless heat energy that needs to be removed from the superheated ball in order to reduce the temperature from T(r, t) to the melting temperature (which, for the idealised problem, is Tmelt = 0). The second term, Q2, is the dimensionless latent heat energy required to melt the solid ball. If Q = Q1(0) − Q2(0) > 0, then there is initially more heat energy in the solid than is required to melt it. There is nowhere for this heat to escape and, as such, the condition Q > 0 necessarily leads to finite-time blow-up51 with , , as .

Returning to our problem, equations (9)–(15) with , we can construct a similar argument, although the details are slightly more complicated. In this case we use Tmelt(r) = 1 − 1/r in our definition of Q1. From the time t = tI, the melting proceeds with T > 1 − 1/R(t) for 0 < r < R(t), and so have the lower bound Q1 > 2πR(t)2/3. Now, some straightforward algebra shows that Q2 < 2πR(t)2/3 for R(t) < 1/(2β), so we conclude that there is a finite particle radius (1/(2β)) below which it is certain that Q1 > Q2. This excess in heat energy is amplified by the additional heat that is conducting in from the liquid phase. Thus we see that, while it is likely that Q1 < Q2 initially, and even for much of the melting process, it is unavoidable that Q1 > Q2 at some point in time. Therefore, for the problem (9)–(15) with , blow-up is inevitable.

With all this in mind, our main result from the full model (9)–(15) with is that including interfacial kinetic effects in the Gibbs–Thomson equation regularises the singularity at R = Rc so that the solution does not blow-up, but is instead regular for all time up until complete melting.

In practice, is it likely to be difficult to measure the precise value of , however we can assume that it is small, and so its effects only come into play when the interface speed is very large. In terms of mathematical modelling, the key point is that the unphysical singular behaviour that occurs in the model for can be regularised by including some kinetic term, no matter how small. Thus by including some in the model, we find results that are physically acceptable for all time, but are still consistent with observations of abrupt melting.

This regularisation is reminiscent of what happens when equivalent kinetic effects are included in the simple model (18)–(22) described above41. From a mathematical perspective, we expect the formal asymptotics of the singular limit to follow the work of King & Evans41, although the details are likely to be even more cumbersome. We do not pursue these issues further here.

Much of the predictive capability of using the continuum model considered here relies on the accurate measurement of the quantity ω in the Gibbs–Thomson rule (1), which depends on the surface tension σ* if using the model (2). We use this quantity as our representative length scale, so each dimensionless unit of length in our results is equivalent to the dimensional length ω. Our time scale defined in equation (8) is proportional to ω2. As mentioned in the Introduction, there are a range of models for σ* that are used in the literature, not all of which agree with the others (for lead there are suggested1920 values for σ* between σ* = 0.031 Jm−2 and 0.145 Jm−2). Indeed, there are many other models for ω that do not follow equation (1)21. A welcome advance would be experimental procedures to refine our quantitative understanding of equation (1) (and ω in particular), especially for metal particles.

This point also goes to the issue of whether it is appropriate to apply our continuum model to nanoscaled phenomena. While it is obvious that a continuum model cannot hold when the metal particles contain only a few atoms, it is not clear what the minimum number of atoms or molecules needs to be. Arguments about which factors affect the limitations of continuum models are summarised by Font & Myers12, who themselves choose to truncate their results for gold particles at R* = 1 nm. Detailed discussion about the validity of thermodynamic arguments is provided in Guisbiers et al.48, while Aguado & Jarrold.54 contains a review of the physics of melting metal clusters on a very small scale (a couple of hundred of atoms or less per cluster). The precise relationship between our dimensionless results and real dimensional quantities depends heavily on the measured values of ω. More certainty about these measurements will inform the applicability of our predictions regarding abrupt melting, and so on.

We close by mentioning the effects of density on our model. Instead of including kinetic effects in the Gibbs–Thomson rule, we may attempt to regularise the singular behaviour by allowing the model to account for different densities in the liquid and solid phases. In this case, the governing equations are generalised to include an advective term in equation (9) and an additional term in the Stefan condition (11). The Gibbs–Thomson condition (1) remains the same. Indeed, a very recent study focuses on such a generalised model with different densities47, but does not elaborate on temperature profiles or consider the near blow-up regime. We may expect that the density difference would regularise the blow-up described in this study, however such a model with equation (1) (and no kinetic term) would still carry with it the same limiting features, such as predicting negative melting temperatures, and the need to cut off the model at an arbitrary particle radius. We leave the study of these issues for further research.

Methods

The dimensionless problem defined by equations (9)–(15) is solved numerically by first applying front-fixing transformations, and then using the method of lines with finite difference spatial discretisations. We approximate the semi-infinite liquid domain R < r < ∞ with the finite domain R < r < rout, making sure that rout is large enough such that all results are independent of it (typically, if Rinit ≈ 24 then we choose ), and the far field condition (14) is approximated with The initial conditions T,init in 0 ≤ r ≤ Rinit and T,init in Rinit ≤ r ≤ rout can lead to unusual effects, such as the development of a ‘chilled region' where the solid region first swells, before melting. While interesting, this phenomenon is not the focus of the current study, and we choose an initial condition such that this does not occur. The initial temperature is chosen to be as follows: in the inner and outer phases, the temperatures are predominantly at the constant temperatures T,init and T,init, respectively, except in the interior layer where we choose the temperature to smoothly transition to the melting temperature (12).

The Landau-type transformations with have the effect of fixing the shrinking solid domain 0 ≤ r ≤ R(t) to 0 ≤ ξ ≤ 1 and the expanding liquid domain R(t) ≤ r ≤ rout to 0 ≤ ζ ≤ 1. Next, a uniform mesh is introduced in each domain, and we replace the spatial derivatives in the new transformed system with second order finite difference approximations. We use one-sided differences to discretise the spatial derivatives at the centre of the solid core ξ = 0 and on the moving boundary ξ = 1, ζ = 0. The Dirichlet boundary condition at rout provides the temperature at ζ = 1.

The result is a semidiscrete system of coupled nonlinear ordinary differential equations in time, which is solved using the built in fully implicit Matlab solver ode15i for the temperature at the mesh points, as well as the position of the moving interface. Typically the temporal stepsize is approximately Δt = 10−3, but the solver allows for a variable stepsize. This is particularly useful for later times where the solver can decrease the stepsize to Δt = 10−11 if need be, such as near blow-up where large rates of change are encountered. While we are solving for the transformed temperatures φ and ψ in the fixed domains 0 ≤ ξ ≤ 1 and 0 ≤ ζ ≤ 1, recall that the original temperatures T and T are on the time varying domains 0 < r < R(t) (shrinking) and R(t) < r < rout (expanding), respectively. This means that as time progresses, the mesh in the outer domain becomes more coarse. Further, the resolution of the mesh in the inner shrinking domain will become more fine, which is important as it is this inner phase which encounters rapid changes in temperature and large spatial derivatives, particularly during the onset of finite-time blow-up. For these reasons, all results were generated with 10000 nodes in each phase, and the solutions were observed to have converged within visual accuracy for both the temperature and interface profiles.

The solver ode15i allows the user to specify the sparsity pattern of the Jacobian matrix, so that it may be formed efficiently using shifted evaluations and forward difference quotients55. With this option in place, accurate simulations with tens of thousands of mesh nodes can be performed with less than two minutes of runtime on a standard desktop machine.

Author Contributions

J.M.B. and S.W.M. conceived and designed the project. J.M.B. and T.J.M. developed the numerical algorithms. J.M.B. conducted the simulations. All authors analysed the results and wrote the manuscript.