The one-dimensional Stefan problem with non-Fourier heat conduction

TL;DR

This paper models one-dimensional solidification using non-Fourier heat conduction theories, revealing that heat flux measurements better indicate non-classical heat transport during phase change.

Contribution

It introduces a reduced ODE system for non-Fourier solidification, linking heat flux to non-classical heat conduction effects during phase change.

Findings

Good agreement between reduced model and simulations

Small differences in silicon solidification kinetics

Significant deviations in heat flux evolution

Abstract

We investigate the one-dimensional growth of a solid into a liquid bath, starting from a small crystal, using the Guyer-Krumhansl and Maxwell-Cattaneo models of heat conduction. By breaking the solidification process into the relevant time regimes we are able to reduce the problem to a system of two coupled ordinary differential equations describing the evolution of the solid-liquid interface and the heat flux. The reduced formulation is in good agreement with numerical simulations. In the case of silicon, differences between classical and non-classical solidification kinetics are relatively small, but larger deviations can be observed in the evolution in time of the heat flux through the growing solid. From this study we conclude that the heat flux provides more information about the presence of non-classical modes of heat transport during phase-change processes.

Tables

Table 1. Table 1 : Values of the thermophysical properties of silicon used in this study [ 51 , 56 , 53 , 54 ] .

Property	Symbol	Value	SI Units
Thermal conductivity	$k^{\ast }$	22.1	W/m$\cdot$K
Specific heat capacity	$c^{\ast }$	1032	J/kg$\cdot$K
Density	${\rho }^{\ast }$	2296	kg/m${}^{\text{3}}$
Thermal diffusivity	${\alpha }^{\ast }$	9.33$\times$10${}^{\text{-6}}$	m${}^{\text{2}}$/s
Phase change temperature	$T_{\text{f}}^{\ast }$	1687	K
Latent heat	$L_{\text{m}}^{\ast }$	1787$\times$10${}^{\text{3}}$	J/kg
Relaxation time	${\tau }_R^{\ast }$	32.16	ps
Mean free path	${\ell }^{\ast }$	8.79	nm

Table 2. Table 2 : Values of the time and flux scales and the dimensionless parameters γ 𝛾 \gamma and η 2 superscript 𝜂 2 \eta^{2} depending on the size of the seed crystal. The Stefan number has been fixed to 10.

$s_c^{\ast }$ [nm]	${\tau }_D^{\ast }$ [ps]	$Q_0^{\ast }$ [W/mm${}^{\text{2}}$]	$\gamma$	$\eta$
2	0.43	1.9	75	7.7
5	2.7	0.77	12	3.0
10	11	0.38	3.0	1.5

Equations

$$c^{\ast }{\rho }^{\ast }\frac{\partial T^{\ast }}{\partial t^{\ast }}+\frac{\partial Q^{\ast }}{\partial x^{\ast }}=0,$$

$$Q^{\ast }=-k^{\ast }\frac{\partial T^{\ast }}{\partial x^{\ast }},$$

$${\tau }_R^{\ast }\frac{\partial Q^{\ast }}{\partial t^{\ast }}+Q^{\ast }=-k^{\ast }\frac{\partial T^{\ast }}{\partial x^{\ast }}+3{\ell }^{\ast 2}\frac{{\partial }^2{Q}^{\ast }}{\partial x^{\ast 2}},$$

$$T^{\ast }=T_{\text{e}}^{\ast },\text{at }{x}^{\ast }=0^{\ast }.$$

$$T^{\ast }=T_{\text{f}}^{\ast },\text{at }{x}^{\ast }=s^{\ast },$$

$${\rho }^{\ast }{L}_m^{\ast }\frac{\mathrm{d}{s}^{\ast }}{\mathrm{d}{t}^{\ast }}=-Q^{\ast },\text{at }{x}^{\ast }=s^{\ast },$$

$$s=s_c^{\ast },T^{\ast }=T_{\text{f}}^{\ast },Q^{\ast }=0,\text{at }{t}^{\ast }=0.$$

$${\tau }_R^{\ast }\frac{\partial Q^{\ast }}{\partial t^{\ast }}+Q^{\ast }=-k^{\ast }\left(\frac{\partial T^{\ast }}{\partial x^{\ast }}+\frac{3{\ell }^{\ast 2}}{{\alpha }^{\ast }}\frac{{\partial }^2{T}^{\ast }}{\partial x^{\ast }\partial t^{\ast }}\right),$$

$${\tau }_R^{\ast }\frac{{\partial }^2{T}^{\ast }}{\partial {t^{\ast }}^2}+\frac{\partial T^{\ast }}{\partial t^{\ast }}={\alpha }^{\ast }\frac{{\partial }^2{T}^{\ast }}{\partial {x^{\ast }}^2}+3{\ell }^{\ast 2}\frac{{\partial }^3{T}^{\ast }}{\partial x^{\ast 2}\partial t^{\ast }}.$$

$$x=\frac{x^{\ast }}{s_c^{\ast }},t=\frac{t^{\ast }}{{\tau }_D^{\ast }},T=\frac{T^{\ast }-T_{\text{f}}^{\ast }}{\Delta T},Q=\frac{Q^{\ast }}{Q_0^{\ast }},s=\frac{s^{\ast }}{s_c^{\ast }}.$$

$$T_t+Q_x$$

$$\gamma Q_t+Q$$

$$T=-1,$$

$$T=0,$$

$$\beta s_t=-Q,$$

$$T=Q=0,s=1,$$

$$\gamma Q_t+Q=-T_x-{\eta }^2{T}_{xt},$$

$$\gamma T_{tt}+T_t=T_{xx}+{\eta }^2{T}_{xxt}.$$

$$\frac{\partial }{\partial t}\mapsto \frac{\partial }{\partial t}-\xi \frac{s_t}{s}\frac{\partial }{\partial \xi },\frac{\partial }{\partial x}\mapsto \frac{1}{s}\frac{\partial }{\partial \xi },\frac{{\partial }^2}{\partial x^2}\mapsto \frac{1}{s^2}\frac{{\partial }^2}{\partial {\xi }^2}.$$

$$\gamma v_t-\gamma \frac{s_t}{s}\xi v_{\xi }+v-{\eta }^2\frac{1}{s^2}{v}_{\xi \xi }-\frac{1}{s^2}{u}_{\xi \xi }$$

$$v-u_t+\frac{s_t}{s}\xi u_{\xi }$$

$$u=-1,$$

$$u=0,$$

$$v=0,$$

$$v+\frac{s_t}{s}{u}_{\xi }=0,$$

$$\gamma w_t-\gamma \frac{s_t}{s}\xi w_{\xi }+w-{\eta }^2\frac{1}{s^2}{w}_{\xi \xi }=-\frac{1}{s}{u}_{\xi }.$$

$$w_{\xi }$$

$$w$$

$$\beta s_t=-w,\text{at }\xi =1.$$

$$u\approx -1+\text{erf}\left(\frac{\xi }{2\sqrt{\zeta t_0}}\right),$$

$$v\approx -\frac{\xi }{2\sqrt{\pi \zeta t_0^3}}\exp \left(-\frac{{\xi }^2}{4\zeta t_0}\right),$$

$$w\approx -\sqrt{\frac{\zeta }{\pi t_0}}\exp \left(-\frac{{\xi }^2}{4\zeta t_0}\right).$$

$$T(x,t)=\frac{x}{s}-1.$$

$$Q_t=-{\gamma }^{-1}\left[\frac{1}{s}+\left(\frac{{\eta }^2}{\beta s^2}+1\right)Q\right].$$

$$Q(0)=-{\eta }^2/\gamma .$$

$$⟨Q⟩=\frac{1}{s}{\int }_0^{s}Qdx.$$

Full text

1 Introduction

Heat conduction at small length or time scales is dominated by effects which are not well described by the classical equations [1]. For instance, it has been observed that the thermal conductivity, which was introduced as a material-intrinsic property and therefore independent of size or shape, exhibits size-dependent behaviour at small length scales [2]. The understanding and correct description of the underlying thermal physics is of vital importance for future applications of nanostructures. For instance, an incorrect description of the thermal response of nanoscale devices can lead to melting and eventually device failure [3]. Novel applications such as laser melting involve phase change at very short time scales [4, 5, 6].

When the characteristic length or time scales of a physical system are comparable to the characteristic length or time scales of the heat carriers, non-classical effects are expected to dominate the heat-transfer mechanisms [1, 7, 8, 9, 10]. In semiconductors such as silicon or germanium, these energy carriers are predominantly phonons, which may be understood as lattice vibrations. These quasi-particles have two associated important physical quantities: the mean free path and the thermal relaxation time. From a microscopic point of view, these represent the mean distance and time that phonons travel between collisions. We can then distinguish two thermal relaxation times, normal and resistive, depending on whether or not phonon momentum is conserved after the collision.

A number of different models extending the classical Fourier law can be found in the literature. In this paper we focus our attention to the Maxwell-Cattaneo equation [11, 12] and the Guyer-Krumhansl equation [13, 14]. Other popular continuum models are the dual phase lag model [15, 16] or the thermomass model [17, 18, 19]. The first incorporates a time lag into all the terms appearing in Fourier’s law, although the physical interpretation of the new parameters remains unclear. The latter uses the Einstein energy-mass relation to describe heat conduction as the motion of a gas formed by particles with non-zero mass called thermons. In this case, a non-linear term emerges in the governing equations, which complicates the mathematical analysis. In contrast to these models, the Guyer-Krumhansl and Maxwell-Cattaneo models involve linear equations where the parameters have a clear physical meaning.

The Maxwell-Cattaneo equation was proposed in the 1950s to remove the infinite speed of heat propagation implicit to Fourier’s law. Mathematically, this model extends Fourier’s law by incorporating a new term which accounts for the time lag between the imposition of a temperature gradient and the creation of a heat flux. This new term leads to the flux depending on the history of the temperature gradient and is thus referred to as a memory term. In the 1960s Guyer and Krumhansl [13, 14] proposed an extension to Fourier’s law derived from the Boltzmann transport equation (BTE) [20] which accounts for memory and non-local effects. Whereas memory effects are important for short times or high frequencies, non-local effects become relevant at small length scales. Although this formalism was initially thought to be applicable only in low-temperature situations ($\sim 2$ K), recent experimental observations show that it is also suitable for describing heat transfer at room temperature or above [8, 21, 22]. Jou et al. [23] later showed that both the Maxwell-Cattaneo and the Guyer-Krumhansl equations can be derived using the framework of extended irreversible thermodynamics (EIT).

The classical mathematical description of a phase-change process, known as the Stefan problem, is based on Fourier’s law [24]. Some authors have incorporated the Maxwell-Cattaneo equation into the formulation of the Stefan problem to model rapid solidification processes [25], cryopreservation [26, 27], cryosurgery for lung cancer treatments [28], and nanoparticle melting [29]. Other authors have studied the Stefan problem with Maxwell-Cattaneo conduction from a point of view of mathematical analysis [30, 31, 32]. Recently, Hennessy et al. [33] performed a detailed asymptotic analysis of the one-dimensional Stefan problem with Guyer-Krumhansl conduction, where it is shown that non-classical effects can lead to important differences in the solidification kinetics with respect to Fourier’s law. Sobolev [34, 35] solved different non-Fourier formulations under the assumption of constant interface velocity. To account for non-local effects, some authors [36, 37, 38] recently proposed a formulation of the Stefan problem based on Fourier’s law with an effective thermal conductivity dependent on the size of the solid. From a microscopic point of view, Font and Bresme [39] published a work where non-classical effects in a melting process were observed using molecular dynamics.

In this paper, we aim to provide tools for detecting non-classical heat transport mechanisms in phase-change processes. For this we study the solidification kinetics in a one-dimensional geometry using the Maxwell-Cattaneo and Guyer-Krumhansl models of heat conduction. The solid is assumed to grow from a small seed crystal, which splits the process into two major regimes, with the first capturing the heat conduction through the seed crystal and the second describing the solidification kinetics. Our asymptotic analysis shows that these two time regimes can be split into multiple sub-regimes, with each of them describing different physical phenomena. We show that, although both conduction models describe the initial transport of heat through the seed crystal in different ways, the solidification process can be captured by reducing the full problem to the same system of ordinary differential equations, which provides excellent agreement with numerical simulations. The results are then applied to the solidification of silicon, where we conclude that the motion of the interface does not necessarily give information about non-Fourier heat transport. Conversely, the evolution of the heat flux is able to provide valuable insights about non-classical modes of heat transfer.

2 Mathematical formulation

To elucidate the roles of non-Fourier conduction mechanisms during solidification, we consider a simple one-dimensional geometry consisting of a liquid bath, which is initially at the phase change temperature $T^{*}_{\text{f}}$, where the ∗ notation denotes dimensional quantities. Due to a cold temperature $T^{*}_{\text{e}}$ in the environment in contact with the bath at $x^{*}=0$, the liquid starts solidifying and a solid phase grows into the bath, occupying the space $[0,s^{*}(t^{*})]$. A schematic of the physical scenario is depicted in Fig. 1. By assuming that the liquid phase is at the phase-change temperature, any decrease in temperature will cause solidification to occur. In this way we can focus on the heat transfer in the solid, which drives the solidification process and is where the non-classical phenomena are expected to occur.

2.1 Governing equations

Conservation of energy in the solid is described by

[TABLE]

where $Q^{*}$, $T^{*}$, $c^{*}$ and $\rho^{*}$ are respectively the heat flux, the temperature, the specific heat capacity and the density of the solid. In the classical formulation, the mathematical description is completed by assuming that the heat flux is determined by Fourier’s law, which in one dimension is

[TABLE]

where $k^{*}$ is the bulk thermal conductivity of the solid. In the present study this will be replaced by the Maxwell-Cattaneo equation (MCE) or the Guyer-Krumhansl (GKE). In one spatial dimension, the GKE takes the form

[TABLE]

where $\tau^{*}_{R}$ is the resistive thermal relaxation time, which is the mean time between consecutive phonon collisions which do not conserve momentum, and $\ell^{*}$ represents the phonon mean free path, i.e. the mean distance between two collisions. The factor 3 in the last term appears due to different second-order terms collapsing to a single contribution in one dimension. The first term of the left-hand side and the second term of the right-hand side of (3) represent respectively memory and non-local effects. The prior are important at small time scales and therefore critical for the initial stages of the solidification process. Non-local effects appear due to two reasons. Firstly, they are caused by the non-negligible rate of normal scattering, i.e. phonon collisions are more likely to conserve their momentum and relaxation effects are delayed. Secondly, non-local effects are important when the characteristic size of the system is comparable to the mean free path.

The GKE can be reduced to the MCE by neglecting the term corresponding to non-local effects, i.e. by setting $\ell^{*}=0$ in Eq. (3). The two approaches, although similar, are based on different frameworks: the GKE is derived from the BTE whereas the MCE is a phenomenological approach to account for finite heat carrier speed. However, it was later shown that both can be derived within the EIT framework [23, 40].

2.2 Boundary and initial conditions

At the boundary $x^{*}=0$ we impose a fixed-temperature condition given by

[TABLE]

A boundary condition of the form of Eq. (4) predicts an initially infinite rate of phase change [24, 41, 42, 43] when the growing phase initially has zero size. In practice, other forms of boundary conditions should be applied, such as the Newton cooling condition. This has been used in the modelling of melting of nanoparticles [44] and nanowires [45], and leads to a significant increase of the melting times. The cooling condition has also been used in a recent study by Hennessy et al. [33], where a thorough asymptotic analysis of the Guyer-Krumhansl-Stefan problem is performed. The same authors also considered a cooling condition accounting for memory effects when studying the melting behaviour of nanoparticles using the MCE [29]. In this problem, we will assume that solidification begins from a seed crystal of finite size and thus the condition (4) can be used without inducing an unphysical infinite rate of phase change.

At the moving interface we impose that the temperature is equal to the phase-change temperature,

[TABLE]

which implies continuity in temperature across the interface. Other authors claim that a jump condition for the temperaure should be utilized at the interface [30, 31, 46]. Based on a diffuse interface model, Hennessy et al. [29] show that, when the MCE is used to model the melting of nanoparticles, the jump condition avoids the phenomenon of supersonic melting [35, 25, 47], which may appear when (5) is used. In the case of nanoparticle melting, the onset of supersonic phase change is due to the spherical geometry of system and the fact that less energy is required to melt the surface of smaller solid cores. Thus, the rate of melting increases and eventually becomes unbounded as the radius of the nanoparticle decreases to zero. As there is no mechanism for supersonic solidification in this one-dimensional Cartesian problem, it is sufficient to consider the case of temperature continuity, which is derived from the jump condition in the limit of slow phase change.

In its simplest form, the Stefan condition is given by

[TABLE]

where $L^{*}_{m}$ is the latent heat. We assume that solidification starts from an initial seed crystal of size $s_{c}^{*}$, formed, for instance, by heterogeneous nucleation, where the temperature is equal to the temperature of the liquid. Since we are using continuum models to describe the growth of the solid, additional assumptions on the size of the seed crystal need to be made. We will assume that its thickness is of the order of nanometres, large enough for continuum theory to hold [48]. The initial conditions are therefore

[TABLE]

2.3 Alternative formulations

Conservation of energy (1) and the GKE (3) can be combined to give

[TABLE]

with $\alpha^{*}=k^{*}/(\rho^{*}c^{*})$ being the bulk thermal diffusivity. This form shows that non-local effects are related to instantaneous changes of the temperature gradient, whereas memory effects are described temporal changes of the flux (which can be related to the time integral of the temperature gradient). In the special case where $\tau^{*}_{R}=3\ell^{*2}/\alpha^{*}$, it is easy to see that Fourier’s law is recovered from (8). This phenomenon has been termed Fourier resonance [49, 50]. Alternatively, it is possible to eliminate the flux from the equation to obtain the Guyer-Krumhansl heat equation (GKHE)

[TABLE]

In the limit $\ell^{*}\to 0$ we obtain the heat equation for the MCE, which in the literature is termed the hyperbolic heat equation (HHE).

2.4 Non-dimensional formulation

Since the classical solidification kinetics are expected to be recovered for sufficiently large times, we will employ the typical scales used to non-dimensionalise the classical equations. In this way, the importance of non-classical contributions will be described by new dimensionless parameters.

The natural length scale of the problem is the initial size of the solid, $x^{*},s^{*}=O(s^{*}_{c})$. Time is then scaled as $t^{*}=O(\tau^{*}_{D})$ with $\tau^{*}_{D}=s^{*2}_{c}/\alpha^{*}$, which corresponds to the time scale for diffusion through the initial seed crystal. The temperature scale is defined by the difference between the freezing and environment temperatures, $\Delta T=T^{*}_{\text{f}}-T^{*}_{\text{e}}$. Balancing terms in (1) then requires $Q^{*}=O(Q^{*}_{0})$, where $Q^{*}_{0}=k^{*}\Delta T/s^{*}_{c}$. Based on these typical scales we introduce the following non-dimensional variables:

[TABLE]

In these variables, conservation of energy (1) and the GKE (3) become

[TABLE]

where we have introduced the subscript notation to denote derivatives. The new parameters in (11b) are $\gamma=\tau^{*}_{R}/\tau^{*}_{D}$ and $\eta=3^{1/2}\ell^{*}/s_{c}^{*}$, which represent the dimensionless forms of the relaxation time and mean free path. In the literature, $\gamma$ is often called the Cattaneo number, whereas $\eta$ is proportional to the ratio of the mean free path to the initial size of the system and termed the Knudsen number.

The non-dimensional boundary, Stefan, and initial conditions become

[TABLE]

where $\beta=L_{m}^{*}/(c^{*}\Delta T)$ is the Stefan number and represents the ratio of latent heat to sensible heat.

In the dimensionless formulation, (8) becomes

[TABLE]

and Fourier resonance appears if the non-classical dimensionless parameters coincide, $\gamma=\eta^{2}$. Finally, the dimensionless GKHE (9) then takes the form

[TABLE]

2.5 Parameter estimation

Silicon is a semiconductor of high interest for many nanoscale applications. The phase-change temperature and latent heat of silicon are $T^{*}_{\text{f}}=1687$ K and $L^{*}_{\text{m}}=1787\times 10^{3}$ J/kg [51]. Since the coefficient of thermal expansion is of the order of $10^{-6}$ K*-1* [52], we can assume a constant value for the density. Based on Mills and Courtney [51], we take $\rho^{*}=2296$ kg/m3. Near the phase-change temperature, the thermal conductivity and specific heat capacity of silicon are $k^{*}=22.1$ W/m$\cdot$K and $c^{*}=1032$ J/kg$\cdot$K [53, 54]. The thermal diffusivity is therefore $\alpha^{*}=9.33\times 10^{-6}$ m${}^{\text{2}}$/s.

Determining the relaxation time $\tau^{*}_{R}$ and the mean free path $\ell^{*}$ is a more difficult task. Here we use the kinetic collective model (KCM) [55], which has been demonstrated to be a powerful tool for predicting the thermal properties of various materials [56, 8]. Since the calculation time increases and the accuracy decreases as we increase the temperature, in Fig. 2 we show results of the KCM up to 1000 K and then fit these values with an equation of the form $a\cdot(T^{*})^{-b}$ to find approximate values of the non-classical parameters near $T^{*}_{\text{f}}$. We find $a\approx 1.16\times 10^{-7}$ and $b\approx 1.11$ for the relaxation time, whereas for the mean free path the best fit is obtained taking $a\approx 4.58\times 10^{-5}$ and $b\approx 1.16$. From these expressions we find $\tau^{*}_{R}\approx 32.16$ ps and $\ell^{*}\approx 8.79$ nm at $T^{*}=1600$ K.

Independent of the underlying model, the dynamics of the solidification process are affected by the Stefan number $\beta$, which can be parametrized in terms of the temperature scale as $\beta=\mathcal{T}^{*}/\Delta T$. For silicon we find $\mathcal{T}^{*}\approx 1731.6$ K, hence even by setting the environmental temperature to [math] K (thus $\Delta T=1687$ K) we find $\beta>1$. Since our study involves an external temperature that is far above absolute zero, we therefore assume $\beta\gg 1$, in accordance with previous studies.

The Cattaneo number $\gamma$ represents the ratio of the time scales of the ballistic and diffusive heat transport regimes. For $\gamma\ll 1$, heat transfer is dominated by diffusion, whereas it becomes ballistic for $\gamma\gg 1$. The Cattaneo number may be parametrized as $\gamma=(L^{*}/s_{c}^{*})^{2}$, where $L^{*}=\sqrt{\alpha^{*}\tau^{*}_{R}}$ is an effective length where memory effects become non-negligible. Upon using the values given in Table 1, we find $L^{*}\approx 17.32$ nm. The Knudsen number $\eta$ corresponds to the ratio of the size of the initial crystal seed to the mean free path and therefore characterises the impact of non-local effects on the initial transport of heat. Using the values in Table 1 we obtain $\eta\approx\ell^{*}_{\text{eff}}/s_{c}^{*}$, where $\ell^{*}_{\text{eff}}=3^{1/2}\ell^{*}\approx 15.23$ nm is the effective mean free path. Since we assume the size of the seed crystal to be of the order of nanometres, we expect $\gamma,\eta\gg 1$ and thus heat conduction through the solid will initially be dominated by non-classical effects.

Finally, the typical scale of the flux may be parametrized in terms of the Stefan number and the initial size of the seed crystal, $Q_{0}^{*}=(k^{*}L^{*}_{m}/c^{*})\times(\beta s_{c}^{*})^{-1}$. For silicon, $k^{*}L_{m}^{*}/c^{*}\approx 3.8\times 10^{-4}$ W/m.

3 Numerical solution

To solve the problem numerically, we transform the growing domain $[0,s(t)]$ to the unit interval $[0,1]$ and then use finite differences to reformulate the problem as a set of algebraic equations.

3.1 Formulation in a fixed domain

Upon introducing the alternative space coordinate $\xi=x/s(t)$, which transforms the growing domain into the unit interval, the derivatives become

[TABLE]

To avoid performing this transformation on the second-order time derivative in (14), which becomes a particularly complicated expression, we introduce an auxiliary variable $v=T_{t}$ to write the problem solely in terms of first-order time derivatives. Hence, we let $T(x,t)=u(\xi,t)$, $T_{t}(x,t)=v(\xi,t)$ and $Q(x,t)=w(\xi,t)$. In these variables, $u$ and $v$ are determined by

[TABLE]

The transformed flux $w$ can be determined by the transformed GKE

[TABLE]

Note, in principle no boundary conditions for $w$ are needed since we could use (13) to write (17a) in terms of $w$, $w_{t}$, and derivatives of $u$ and $v$, the latter of which are already known from (16).

Finally, the transformed Stefan condition is

[TABLE]

From this formulation we obtain the numerical solution by solving (16), (17) and (18) at each time step using finite differences. Further details about the numerical scheme are provided in the Supplementary Material.

3.2 Small-time approximations

Numerical experimentation has shown that replacing the initial conditions for $u$, $v$, and $w$ with approximate solutions that are valid for arbitrarily small times $t\ll 1$ greatly improves the stability of the computational scheme when solving the GKE (for $\eta\neq 0$). By carrying out a small-time analysis of the model, which is described in the Supplementary Material, we find that the solution at time $t_{0}\ll 1$ can be approximated by

[TABLE]

where $\zeta=\eta^{2}/\gamma$ represents an effective thermal conductivity. When solving the MCE, the original initial conditions $u(\xi,0)=v(\xi,0)=w(\xi,0)=0$ can be used.

4 Reduction of the equations

As shown in the asymptotic analysis in the Supplementary Material, by the time that solidification occurs, the heat flux becomes approximately uniform in space: $Q_{x}\simeq 0$. This leads to the so-called quasi-steady regime, since differentiating the GKE (11b) with respect to $x$ then leads to $T_{xx}\simeq 0$ and hence time is present only due to the boundary condition (12b). This quasi-steady regime also appears in the classical formulation under the assumption $\beta\gg 1$, allowing the system to be reduced to a simple ordinary differential equation (ODE) for the position of the interface [57]. Our asymptotic analysis reveals that a similar procedure can be considered here even in the case of strongly non-classical effects. Applying the boundary conditions (12a) and (12b) implies that the temperature is given by

[TABLE]

Upon using (20) in (13), and then (12c) to eliminate the arising term with $s_{t}$, yields

[TABLE]

The asymptotic analysis, described in the Supplementary Material, shows that the appropiate matching condition is given by

[TABLE]

Note, this initial condition is valid for both the Guyer-Krumhansl and Maxwell-Cattaneo models, since the limit $\eta\to 0$ gives the initial condition $Q(0)=0$ which corresponds to the small time behaviour of the MC-Stefan problem. Hence, based on the asymptotic analysis we are able to reduce a problem involving two partial differential equations in a growing domain into a pair of coupled ordinary differential equations (12c) and (21).

5 Results and discussion

We are now in position to compare the different solidification kinetics in the limits of dominant non-classical effects and to assess the validity of the reduced formulation. Afterwards we apply the results to the solidification of silicon using the parameter values discussed in Sec. 2.5.

To compare the classical and non-classical heat conduction models, it will be insightful to consider the mean flux $\langle Q\rangle$, defined by

[TABLE]

Note, this quantity is uniform in space and, moreover, it coincides with the heat flux $Q$ when the temperature enters in the quasi-steady regime.

5.1 MC conduction with enhanced memory

For $\gamma,\beta\gg 1$ with $\gamma=O(\beta)$ or larger, the asymptotic analysis in the Supplementary Material shows that we can split the solidification process into three time regimes. The first time regime occurs when $t=O(1)$ and captures the initial propagation of a heat wave from the cooled boundary into the seed crystal. The second time regime, $t=O(\gamma^{1/2})$, describes the transport of heat across the entirety of the seed crystal. This regime also captures the reflections of thermal waves that occur when they collide with the boundaries of the crystal. We find that the temperature waves destructively interfere, thus keeping the temperature $O(1)$ in size. However, the flux waves undergo constructive interference and thus the magnitude of the flux scales like $O(t)$. Due to the large Stefan number, no appreciable solidification occurs during the first and second time regimes. In the third time regime, $t=O(\gamma^{1/2}\beta^{1/2})$, the temperature settles into its quasi-steady (linear) profile and solidification begins. Memory effects are relevant during the initial stages of solidification but diminish with time until Fourier’s law is recovered. If $\gamma\gg 1$ but $\gamma\ll\beta$, then memory effects diminish before solidification begins.

Figure 3 shows the solidification kinetics and evolution in time of the mean flux for different values of $\gamma$. Solid lines correspond to numerical simulations, whereas circles refer to the solutions of the reduced formulation. The results are shown for $\beta=10$. Since the asymptotic analysis reveals that for $\gamma\ll\beta$ non-classical effects become negligible by the time the solid starts growing, the case $\gamma=0$ is not shown in Fig. 3(a), as it is indistinguishable from the case $\gamma=1$. However, the mean flux for the $\gamma=0$ case is shown in Fig. 3(b), where it can be clearly observed how the limit $\gamma\to 0$, i.e. the Fourier case, represents a singular limit of the problem. In all cases, the agreement between the numerical simulations and the solutions of the reduced model is excellent. Simulations for larger values of $\beta$ (not shown here) also provide excellent agreement between both methods.

In Fig. 3(a) we can see how the start of the solidification process is delayed as $\gamma$ increases. This is caused by the fact that, for large values of $\gamma$, heat propagates in the form of a thermal waves with finite speed $\gamma^{-1/2}$ for a longer period and solidification begins after the wave hits the right extreme of the seed crystal when $t=O(\gamma^{1/2})$. The asymptotic analysis reveals that solidification begins when $t=O(\gamma^{1/2}\beta^{1/2})$. Hence, for $\gamma=10$ and $\gamma=100$, solidification begins when $t\approx 10$ and $t\approx 32$ respectively, which is confirmed in Fig. 3(a). The case $\gamma=100$ shows an interesting evolution of the interface. On one hand, the start of the interface motion is delayed due to the slower speed of the thermal wave through the crystal. On the other hand, the long-term rate of solidification is higher than for smaller values of $\gamma$ due to the delayed increase of the flux, which reaches its largest value when $t=O(\gamma)$, where its value is larger than the flux corresponding to smaller relaxation times. In this sense, solidification begins later as we increase $\gamma$, but the long-term solidification rate is larger until the classical solidification kinetics, where $s\sim t^{1/2}$, are recovered.

In Fig. 3(b) we observe that, as $\gamma$ decreases, the flux converges earlier to the profile described by the classical formulation. As $\gamma$ increases, the initial MC flux decreases, in accordance with the fact that $Q=O(\gamma)$ for early times. The decrease in the mean flux from its initial value of zero is due to the repeating sequence of constructive wave reflections that occurs at the boundaries of the solid. As $\gamma$ decreases and the thermal-wave speed increases, these reflections occur more frequently, thereby leading to a more rapid decrease in the mean flux. The Fourier case can therefore be interpreted as a limit whereby these reflections occur infinitely frequently, leading to an instantaneous jump in the mean flux from its initial value of zero. The asymptotic analysis also reveals an increase of the flux for $t\gg\gamma^{1/2}$, with the order of magnitude being proportional to $(\beta/\gamma)^{1/2}$. This can be observed in Fig. 3(b) and agrees with the asymptotic analysis.

5.2 GK conduction with enhanced memory

We now assume that $\gamma,\beta\gg 1$ with $\eta=O(1)$. The asymptotic analysis in the Supplementary Material reveals there are four time regimes. When $t\ll 1$, Fourier’s law is recovered with a small effective thermal conductivity (ETC) defined by $\zeta=\eta^{2}/\gamma$. Changes in the temperature and flux from their initial values only occur in a thin layer near $x=0$. When $t=O(1)$, heat propagates towards the liquid as a wave. The form of this wave is analogous to a viscoelastic pressure wave. Contrary to the MC case, here the wavefront (with position given by $x_{f}(t)$) is not sharp but diffuse with a (dimensionless) width that scales like $w=O(\zeta^{1/2}t^{1/2})$. When $t=O(\gamma^{1/2})$ we find that the thermal and flux waves undergo a sequence of reflections at the boundaries of the solid that is analogous to the MC case, where now the width of the diffuse wavefront continuously increases. The choice of the fourth time regime depends on the relative size of $\gamma$ to $\beta$. In all cases, non-local effects become negligible, solidification begins, and the dynamics become analogous to the third time regime of the MC case.

The evolution in time of the solid-liquid interface and the mean flux is shown in Fig. 4 for different values of $\gamma$. The mean flux corresponding to the Fourier model is also shown (dashed, gray line). The remaining dimensionless parameters have been set to $\eta=1$ and $\beta=10$. In addition, we have plotted the case $\gamma=1$, which reproduces the classical solidification kinetics due to the Fourier resonance phenomenon, as shown by the perfect match of the corresponding fluxes in Fig. 4(b). Similarly to the MC case, we note that the agreement between the numerical simulations and the solutions of the reduced model is remarkable.

By comparing Figs. 3(a) and 4(a), it can be observed that the solidification kinetics in the GK case are very similar to those in the MC case. This was already anticipated by the asymptotic analysis, which predicts that non-local effects become negligible by the time solidification begins. In fact, the key differences between the MC and GK formulations with large relaxation times do not appear in the solidification process, since it is essentially captured by the same equations (see Secs. SM3.3 and SM4.1.4), but in the initial conduction through the crystal, as can be observed by comparing Figs. 3(b) and 4(b). However, notice that these initial differences between the MC and GK formulations disappear for $t=O(\gamma^{1/2}\beta^{1/2})$ and hence both fluxes coincide by the time solidification begins. The different behaviour for smaller times is caused by the different underlying physics in each model. The MCE predicts a wave-like heat propagation through the crystal, whereas the GKE initially recovers Fourier’s law with an effective thermal conductivity $\zeta$. In the case of Fourier resonance, the ETC becomes $\zeta=1$ and it coincides with the bulk thermal conductivity.

As revealed by the asymptotic analysis, the flux in this specific limit of the GK model is very similar to the MC model and differences appear only for small times. The evolution of the flux profiles for the MC and GK cases is shown in Fig. 5. To clearly illustrate the wave-like propagation of heat in the MC and GK models, we increase $\gamma$ to 100 in both cases, using $\eta=1$ for the GK case, and fix $\beta=10$ as before. We can observe that differences in both formulations indeed appear only for small times and in the way the moving wave front $x_{f}(t)$ evolves in time. The right panels in Fig. 5 show the time regimes when heat propagates in the form of thermal waves. A sharp wavefront moving at a speed $\gamma^{1/2}$ can be clearly observed in the MC case, in contrast to the GK case, which has a diffuse interface.

5.3 GK conduction with enhanced non-local effects

In this case we assume $\eta,\beta\gg 1$ and $\gamma=O(1)$. From the asymptotic analysis we know that there are three time regimes to consider. In the first regime, $t=O(\eta^{-2})$, Fourier’s law is recovered with a large value of the ETC, $\zeta\gg 1$, and no solidification occurs. Due to the high value of the ETC, this regime describes the propagation of a large heat flux into the seed crystal. The second time regime depends on the relative size of $\eta^{2}$ to $\beta$, but in any case the system enters a quasi-steady regime. In the distinguished limit $\beta/\eta^{2}=O(1)$, the second time regime is given by $t=O(1)$, and the flux is driven by memory and non-local effects, i.e. the history and rate of change of the temperature gradient, rather than the temperature gradient itself. As time increases, the strength of these non-classical transport mechanisms diminishes and solidification decelerates, causing the position of the solid-liquid interface to reach a temporary steady state. However, as the third time regime, $t=O(\beta)$, is entered, the contributions to the flux from the temperature gradient become large enough to balance and then dominate the non-classical mechanisms, re-initiating the process of solidification. For large times, the classical solidification kinetics, $s\sim t^{1/2}$, are recovered. Similar time regimes were observed in Hennessy et al. [33], although in their case the magnitude of the flux for small times is not as large as here due to their use of a Newton condition at $x=0$ rather than a fixed-temperature condition.

The evolution of the interface and the mean flux is shown in Fig. 6, where lines and circles again represent solutions computed numerically and via the reduced formulation. Due to the large variation of the magnitude of $\langle Q\rangle$, the results in Fig. 6(b) are shown using logarithmic scales on both axes. As in the case of enhanced memory effects, we also show the results for $\eta=1$, which reduces to the classical model due to the Fourier resonance phenomenon. In general, there is good agreement between the solution to the reduced formulation and numerical simulations, with some small discrepancies in Fig. 6(a) that may be attributed to the terms being neglected in the asymptotic analysis. Larger discrepancies may be observed in Fig. 6(b) for $t\ll 1$, where the reduced formulation predicts smaller values (in magnitude) of the heat flux than the numerical simulations. This error is due to the reduced model not capturing these small-time regimes.

As predicted by the asymptotic analysis, the initial solidification rate increases as the non-local parameter $\eta$ is increased. This is attributed to the increase in the ETC, which enhances the transport of thermal energy to the solidification front. The rapid initial growth of the solid is solely driven by non-Fourier conduction mechanisms; the contribution to the flux from the instantaneous temperature gradient is negligible for small times. As time increases and the solid grows in size, memory and non-local effects weaken and the solidification rate process undergoes a substantial decrease, as seen in Fig. 6 when $\eta^{2}=100$. The mean flux relaxes to roughly the Fourier contribution, $\langle Q\rangle\simeq-1/s$; however, due to the rapid initial growth of the solid ($s\gg 1$), this contribution is strongly diminished compared to the case of pure Fourier conduction. As a result, when solidification is re-initiated on $O(\beta)$ time scales, it proceeds at a greatly reduced rate.

In summary, the dynamics that occur in the case of GK conduction with strong non-local effects are the opposite of those that occur in the limit of strong memory effects or for MC conduction. In the former case $\eta\gg 1$ with $\gamma=O(1)$, the mean thermal flux is initially much greater and then much smaller than the Fourier prediction, which leads to solidification kinetics that are much faster and then much slower than the Fourier case. However, if $\gamma\gg 1$ with $\eta=O(1)$ (or $\eta=0$), then there is a prolonged period where the mean thermal flux is initially much higher and then slightly less than the Fourier prediction, resulting in relatively slower and then faster solidification kinetics.

5.4 Application to the solidification of Si

Silicon is a material of high interest for nanoscale applications and an obvious choice for the validation of theoretical models. In Fig. 7 we plot the evolution of the solid-liquid interface and the mean flux according to the classical, MC and GK formulations for different sizes of the seed crystal. The corresponding non-classical parameters are given in Table 2. The results have been calculated numerically for $\beta=10$, which corresponds to an initial temperature jump $\Delta T^{*}=173.2$ K. The case $\beta=100$ ($\Delta T^{*}=17.3$ K) leads to all curves in Figs. 7(a), 7(c) and 7(e) collapsing and thus no differences between classical and non-classical solidification kinetics can be observed. To give the results in dimensional form, we have written $\langle Q^{*}\rangle=Q^{*}_{0}\langle Q\rangle$ and $t^{*}=\tau^{*}_{D}t$, with $Q^{*}_{0}$ and $\tau^{*}_{D}$ given in Table 2 for each of the seed crystal sizes.

Figures 7(a), 7(c) and 7(e), where we plot the evolution of the solid-liquid interface, show that the presence of non-classical effects does not cause any major differences in the solidification kinetics. In fact, the Fourier and GK models lead to virtually identical behaviours in the evolution of the solid-liquid interface for the three seed crystal sizes considered. In the case of the MC model, there is a delay in the start of the solidification, which from a physical point of view can be understood from the fact that heat is initially propagated in the form of thermal waves with finite speed, contrary to the Fourier and GK models. Furthermore, the asymptotic analysis of the MC model reveals that in the cases $s_{c}^{*}=2$ nm and $s_{c}^{*}=5$ nm, where $\gamma$ is relatively large, solidification starts for $t^{*}\approx 11.78$ ps and $t^{*}\approx 29.36$ ps respectively, which is in good agreement with Fig. 7(a).

Although differences in the plots corresponding to $s^{*}$ are small, Figs. 7(b), 7(d) and 7(f), which show the evolution in time of the mean flux, suggest that the choice of the underlying model does have a crucial impact on how heat will propagate through the crystal for small times ($<100$ ns). It can be observed how the qualitative behaviour of the Fourier and GK fluxes is similar whereas they have large differences with respect to the MC flux. These differences are a direct consequence of the different behaviours for small times, since in the latter heat is initially propagated in the form of heat waves whereas for the GK and classical models the initial heat propagation corresponds to a diffusive process. Therefore, the results shown in Fig. 7 suggest that non-classical mechanisms can be detected by measuring the total heat flux through the solid.

The similarities between the solidification kinetics of the Fourier and GK models could exist for two reasons. Firstly, the asymptotic analysis shows non-classical effects become less important for the solidification process as the Stefan number increases. Secondly, the ratio $\zeta=\eta^{2}/\gamma=3\ell^{*2}/(\alpha^{*}\tau^{*}_{R})$ does not depend on the initial size of the seed crystal. Using the values provided in Table 1, for silicon we find $\zeta\approx 0.77$, which implies that the similarities may be caused by the system being close to the case of Fourier resonance ($\zeta=1$).

In Fig. 8(a) we show the solidification kinetics for three different seed crystal sizes and for a larger temperature jump of $\Delta T^{*}=500$ K, which reduces the Stefan number to $\beta\approx 3.46$. The results show that the similarities between the Fourier and GK models persist, which indicates that the largeness of the Stefan number was not causing the similarities in the solidification kinetics observed in Figs. 7(a), 7(c) and 7(e).

Thus, we now explore whether these similarities are caused by Fourier resonsance phenomenoma by varying the size of the ETC $\zeta$. In Fig. 8(b) we show how the solidification kinetics vary when changing the order of magnitude of $\tau_{R}^{*}$ while keeping the value of $\ell^{*}$ fixed, which leads to variations in $\zeta$ across two orders of magnitude. In the cases where $\zeta\ll 1$ or $\zeta\gg 1$, the non-Fourier signatures predicted by the asymptotic analysis can be observed. This shows that the similarities between the solidification kinetics in the Fourier and GK models are indeed caused by the fact that $\zeta$ is close to unity. Hence, evidence of non-classical heat transport can be found in measurements of the solidification front only in the cases $\zeta\gg 1$ or $\zeta\ll 1$.

From the results of this study we can conclude that non-classical conduction mechanisms do not lead to significant changes in the solidification kinetics in situations near Fourier resonance, despite the fact that the dimensionless parameters characterising these effects are large. Hence, in this case non-classical heat transport mechanisms are unlikely to be detected from measurements of the solidification front. However, the analyses predict that large deviations can occur in the thermal flux, which is a better quantity to use in the detection of non-classical transport mechanisms and which can be measured both experimentally and in molecular dynamics simulations [36].

6 Conclusion

In this paper we have studied non-classical heat transport in phase-change processes as described by the Maxwell-Cattaneo (MC) and Guyer-Krumhansl (GK) equations via the solidification kinetics of a one-dimensional liquid bath. Based on a detailed asymptotic analysis we proposed a simplified formulation which reduces the moving boundary problem to a set of two ordinary differential equations for the position of the interface and the mean flux through the solid. Comparisons of the results provided by this reduced formulation against numerical simulations have shown excellent agreement. When the non-classical terms become dominant, the study shows that we can observe large deviations from the classical solidification kinetics described by Fourier’s law. However, in GK conduction near Fourier resonance, for example with silicon, the non-classical effects almost cancel each other out. This leads to similar evolution of the phase-change front between the GK and Fourier models. However, the propagation of heat through the solid is affected by the choice of the constitutive equation to a greater extent. We conclude that a possible strategy to determine the presence of non-classical effects is studying the heat conduction through the growing solid rather than the position of the interface, since the prior seems to give more information even in cases that are close to Fourier resonance. Moreover, due to the different small-time behaviours, measuring the flux could also indicate which of the considered models is more realistic, since the MC model does not have Fourier resonance. Although we focused on MC and GK conduction laws, the reduction framework developed here can be applied to alternative non-Fourier models based on continuum theory and used to study other modes of phase change such as melting as well.

Acknowledgements

M.C.S. acknowledges that the research leading to these results received funding from “La Caixa” Foundation and the Ferran Sunyer i Balaguer Foundation. M.G.H and T.G.M received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 707658. T.G.M. acknowledges the support of a Ministerio de Ciencia e Innovación grant MTM2017-82317-P. The authors have been partially funded by the CERCA programme of the Generalitat de Catalunya. The authors thank F. X. Alvarez and L. Sendra for providing the KCM data.

SM1 Introduction

In this Supplementary Material we focus on the numerical solution and asymptotic analysis of the Maxwell-Cattaneo-Stefan and the Guyer-Krumhansl-Stefan problems. The relevant, dimensionless equations are

[TABLE]

For $\eta=0$ the GKE reduces to the MCE and hence the Guyer-Krumhansl-Stefan problem reduces to the Maxwell-Cattaneo-Stefan problem.

SM2 Numerical Scheme

In this section we consider the problem in the alternative variables $T(x,t)=u(\xi,t)$, $T_{t}(x,t)=v(\xi,t)$ and $q(x,t)=w(\xi,t)$, where $\xi=x/s$. The relevant transformed equations have been presented in the main text and take the form

[TABLE]

and

[TABLE]

The system (S2) determines $u$ and $v$, whereas (S3) is solved to determine $w$. Finally, the solid-liquid interface is determined by

[TABLE]

subject to $s(0)=1$.

We discretize the unit interval into $N+1$ equidistant points of the form $\xi_{i}=i\Delta\xi$, where $i=0,...,N$ and $\Delta\xi=1/N$. Similarly, we consider time instants of the form $t_{n}=t_{0}+n\Delta t$, where $t_{0}$ is the starting time and $n\geq 0$.

At each time step we first solve (S2), then use this information to update the heat flux by solving (S3), and finally we update the position of the solid-liquid interface using (S4). For this we discretize the derivatives of $u,v,w$ implicitly, i.e.

[TABLE]

whereas at $\xi=1$ we use

[TABLE]

to retain second-order accuracy in space. The boundary conditions for $w$ are discretized analogously. Upon writing $s\approx s^{n}$ and $s_{t}\approx s_{t}^{n}$, the discretized form of (S2) is a set of $2N+2$ linear equations that can be expressed as

[TABLE]

for certain matrices $A^{n}_{1},A^{n}_{2}\in\mathbb{R}^{2(N+1)\times 2(N+1)}$ and where $\mathbf{e}_{1}=(1,0,....,0)^{T}\in\mathbb{R}^{2(N+1)}$ and $\mathbf{u}^{n+1}=(u^{n+1}_{0},...,u^{n+1}_{N})^{T}$, $\mathbf{v}^{n+1}=(v^{n+1}_{0},...,v^{n+1}_{N})^{T}$. To find $\mathbf{w}^{n+1}$ we discretize (S3) at each grid point, using central differences in interior points and forwards/backwards differences at the boundaries. This gives a linear system for $\mathbf{w}^{n+1}$ of the form

[TABLE]

for certain matrices $B^{n}_{1},B^{n}_{2},B^{n}_{3}\in\mathbb{R}^{(N+1)\times(N+1)}$. Finally, we obtain the updated value of the position and speed of the interface by

[TABLE]

SM3 Asymptotic analysis of the Maxwell-Cattaneo-Stefan problem

In this section we consider (S1) with $\eta=0$ and perform an asymptotic analysis of the problem assuming $\gamma,\beta\gg 1$, which corresponds to a limit where memory effects are strong. In the case when $\gamma=O(1)$ and $\beta\gg 1$, we find that the leading-order solidification kinetics are identical to those predicted by Fourier’s law.

We can combine Eqs. (S1a) and (S1b) to eliminate $Q$, which yields the hyperbolic heat equation (HHE)

[TABLE]

Mathematically, this equation describes the propagation of a temperature wave with a finite speed $\gamma^{-1/2}$. Damping is initially negligible, but appears gradually as the wave travels through the seed crystal. From a microscopic point of view, damping appears as the number of resistive collisions increases.

Heat conduction through the crystal can be split into three time regimes. The first time regime captures the propagation of the disturbance at $x=0$ and the second time regime captures the collision of the wave against the solid-liquid interface. In the third time regime we enter a pseudo-steady state where the flux becomes a function of time and the solidification process eventually begins.

SM3.1 First time regime

Due to the fact that $T=-1$ at $x=0$ and $s(0)=1$ we have $T,s=O(1)$. Additionally, the initial wave-like behaviour of the temperature implies $x=O(\gamma^{-1/2}t)$ and damping effects are negligible for $t\ll\gamma$. Finally, balancing terms in Eq. (S1b) requires $Q=O(\gamma^{-1/2})$. By introducing the scaled variables defined by $t=\hat{t}$, $x=\gamma^{-1/2}\hat{x}$, $Q=\gamma^{-1/2}\hat{Q}$, $T=\hat{T}$ and $s=\hat{s}$, the Stefan condition becomes $\hat{s}_{\hat{t}}=O(\gamma^{-1/2}\beta^{-1})$, hence $\hat{s}=1$ at leading order, which reduces the problem to

[TABLE]

The leading-order problem can be solved using the method of characteristics, giving $\hat{T}=-\mathcal{H}(\hat{t}-\hat{x})$ and $\hat{Q}=-\mathcal{H}(\hat{t}-\hat{x})$, where $\mathcal{H}$ is the Heaviside function. The solution therefore corresponds to a travelling wave that propagates with constant velocity. The position of the wavefront is given by $\hat{x}_{f}(\hat{t})=\hat{t}$. An examination of the large-time behaviour shows that $\hat{x}_{f}=O(\hat{t})$ and $\hat{Q},\hat{T}=O(1)$ as $\hat{t}\to\infty$, which determine the scales for the second time regime.

SM3.2 Second time regime

By looking at the scaled system (S11), a new balance arises for $\hat{t}=O(\gamma)$, which implies $x=O(\gamma^{1/2})$ and hence $x\gg 1$. This suggests that one should first consider the time regime where the heat wave hits the boundary, i.e., we force $x_{f}=O(1)$ and hence $t=O(\gamma^{1/2})$. Additionally, we have $T,s=O(1)$ and $Q=O(\gamma^{-1/2})$ in order to match to the solution in the previous time regime. In the corresponding scaled variables $\tilde{x},\tilde{s},\tilde{t},\tilde{Q}$ and $\tilde{T}$, the Stefan condition becomes ${\tilde{s}}_{\tilde{t}}=O(\beta^{-1})$, hence we still do not have any interface motion at leading order. In these variables, the problem is similar to (S11), with the only difference in the spatial domain,

[TABLE]

with $\tilde{T}\sim-\mathcal{H}(\tilde{t}-\tilde{x})$ and $\tilde{Q}\sim-\mathcal{H}(\tilde{t}-\tilde{x})$ as $\tilde{t}\sim 0$ as matching conditions.

We first solve the leading-order problem for the temperature, which satisfies the wave equation $\tilde{T}_{\tilde{t}\tilde{t}}=\tilde{T}_{\tilde{x}\tilde{x}}$ in the bulk. The matching condition corresponds to a right-moving travelling wave and satisfies the bulk equation. Thus, the solution for $0<\tilde{t}<1$ is given by $\tilde{T}=-\mathcal{H}(\tilde{t}-\tilde{x})$. When $\tilde{t}=1$, this wave collides with the boundary at $\tilde{x}=1$ and reflects. To determine the dynamics of this reflection, we shift the origin using the transformation $\tilde{x}=\tilde{x}^{\prime}+1$ and write the solution for the temperature in terms of an incoming wave $\tilde{T}_{I}$ and a reflected wave $\tilde{T}_{R}$ as $\tilde{T}=\tilde{T}_{I}(\tilde{t}-\tilde{x}^{\prime})+T_{R}(\tilde{t}+\tilde{x}^{\prime})$. Imposing the condition $\tilde{T}=0$ when $\tilde{x}^{\prime}=0$ gives that $T_{R}=-T_{I}$, implying that the reflected wave destructively interferes with the incoming wave. The resulting solution for the temperature is a left-moving travelling wave given by $\tilde{T}=-\mathcal{H}(\tilde{t}+\tilde{x}-2)$ valid for $1<\tilde{t}<2$. When $\tilde{t}=2$, the wave collides with the left boundary at $\tilde{x}=0$, again creating a reflected wave. By writing the temperature as $\tilde{T}=-1+\tilde{T}^{\prime}$, the same analysis can be used to show that the reflected wave destructively interferes with the incoming wave, resulting in a right-moving travelling wave given by $\tilde{T}=-\mathcal{H}(\tilde{t}-\tilde{x}-2)$ for $2<\tilde{t}<3$. As this solution is identical (up to a translation in time) to that for $0<\tilde{t}<1$, the sequence of destructive reflections repeats itself. Importantly, we find that $\tilde{T}=O(1)$ for all times in this regime.

The dynamics for the flux are much more interesting and can be determined by solving (S12a) using the known solutions for the temperature. For $0<\tilde{t}<1$, the flux is also a right-moving travelling wave given by $\tilde{Q}=-\mathcal{H}(\tilde{t}-\tilde{x})$. However, the collision with the boundary at $\tilde{x}=1$ results in constructive interference. This can be seen by writing the solution for the temperature for $1<\tilde{t}<2$ as $\tilde{T}=\tilde{T}_{I}(\tilde{t}-\tilde{x}^{\prime})-\tilde{T}_{I}(\tilde{t}+\tilde{x}^{\prime})$. The corresponding solution for the flux is $\tilde{Q}=\tilde{T}_{I}(\tilde{t}-\tilde{x}^{\prime})+\tilde{T}_{I}(\tilde{t}+\tilde{x}^{\prime})$; thus, the right-moving incoming wave and left-moving reflected wave add in a constructive manner. We find that $\tilde{Q}=-1-\mathcal{H}(\tilde{t}+\tilde{x}-2)$ for $1<\tilde{t}<2$, corresponding to a left-moving travelling wave which takes on values between $-2$ and $-1$. Repeating the analysis shows that the collision with the boundary at $\tilde{x}=0$ is also constructive, resulting in a solution for the flux given by $\tilde{Q}=-2-\mathcal{H}(\tilde{t}-\tilde{x}-2)$ for $2<\tilde{t}<3$, which, again, is a right-moving travelling wave. In this case, the flux ranges from $-3$ to $-2$ across $\tilde{x}$. The sequence of constructive reflections continues, each time reducing the flux by one unit. Thus, the asymptotic behaviour of the flux is given by $\tilde{Q}=O(\tilde{t})$ for $\tilde{t}\gg 1$.

Having determined the behaviour of the temperature and flux for large times, we are now in a position to determine the relevant scales for the next time regime. It turns out there are two main choices depending on the relative sizes of $\beta$ and $\gamma$. One one hand, if $\gamma\ll\beta$, a new balance in the governing equations is obtained by setting $\tilde{t}=O(\gamma^{1/2})$ and $\tilde{Q}=O(\gamma^{1/2})$,

[TABLE]

whereas the Stefan condition becomes

[TABLE]

This scaling implies that the travelling waves start to become damped, which is associated with a relaxation to the classical Fourier profiles; however, solidification has yet to begin. On the other hand, if $\gamma\gg\beta$, then we can balance both sides of the Stefan condition by setting $\tilde{t}=O(\beta^{1/2})$ and $\tilde{Q}=O(\beta^{1/2})$. In this case, the rescaled equations are

[TABLE]

implying that solidification begins before the thermal waves become damped. Both of these scalings result in a quasi-steady scenario whereby the flux is uniform ($\tilde{Q}_{x}=0$) and the temperature profile is linear (which is seen by differentiating the MCEs with respect to $\tilde{x}$).

SM3.3 Third time regime

As previously discussed, the choice of the time scale depends on the relative size of $\gamma$ to $\beta$. For simplicity, we can consider the distinguished limit where $\gamma^{-1/2}\beta^{1/2}=\alpha=O(1)$. Based on the large-time analysis of the previous time regime we take $t=O(\gamma^{1/2}\beta^{1/2})$, $Q=O(\gamma^{-1/2}\beta^{1/2})$ and $x,s,T=O(1)$. In the corresponding scaled variables $\bar{t},\bar{x},\bar{s},\bar{Q}$ and $\bar{T}$, the problem becomes

[TABLE]

with matching conditions given by the previous time regime. Differentiating Eq. (S19a) with respect to $x$ gives that $T$ is linear in space and hence, applying Eqs. (S19c) and (S19d), we find

[TABLE]

Using this expression in (S19a) yields

[TABLE]

subject to the initial condition $\bar{Q}(0)=0$ due to the large time behaviour $\tilde{Q}\sim-\tilde{t}$ in the previous time regime. Note, for $\bar{t}\gg 1$ we find $\bar{Q}\sim-1/(\alpha\bar{s})$ and hence, using (S19e), the classical behaviour $\bar{s}\sim\bar{t}^{1/2}$ is obtained. The reduced model given by (S19e) and (S21) can also be used when $\alpha\ll 1$, corresponding to the case when $\gamma\gg\beta$. From an asymptotic point of view, this case can be studied using matched asymptotic expansions in terms of $\alpha$ through the introduction of additional time and length scales.

SM4 Asymptotic analysis of the Guyer-Krumhansl-Stefan problem

As we show now, the presence of non-local effects leads to a completely different small-time behaviour of the problem. This is due to the fact that governing equations are mathematically very different in the cases $\eta=0$ and $\eta>0$. In the first case we obtain the HHE, which is a hyperbolic equation and describes the propagation of a heat wave, contrary to the GKHE, which is a parabolic equation. In this section we will perform an asymptotic analysis of (S1) assuming that either memory or non-local effects are strong.

SM4.1 Enhanced memory effects

Let us assume that $\gamma\gg\eta=O(1)$, i.e. that heat transport in the seed crystal is initially dominated by memory effects. Additionally we assume $\beta\gg 1$.

The initial heat conduction through the solid can be divided into four time regimes. In the first time regime, heat conduction is described by an effective Fourier law. In the second time regime, both memory and non-local effects enter the leading-order problem, giving rise to a wave-like propagation of heat that is analogous to one-dimensional viscoelastic wave propagation. The third time regime captures the refections of the thermal waves with the boundaries. In the fourth regime, the system tends towards a quasi-steady state and solidification begins. Non-local effects are negligible but memory effects persist.

SM4.1.1 First time regime

At the beginning of the process we have $T,s=O(1)$ due to the boundary condition at $x=0$ and the initial condition for $s$. Assume now $t=O(\epsilon)$, where $\epsilon$ is an arbitrarily small parameter. Balancing terms in (S1b) requires $xQ^{-1}=O(\epsilon)$. Since $\epsilon\ll 1$ and $\gamma\gg 1$, the terms $Q_{t}$ and $Q$ in (S1a) cannot balance. The same happens for the terms $T_{x}$ and $Q_{xx}$ on the r.h.s. of (S1a). Hence, $Q_{t}$ must balance either $T_{x}$ or $Q_{xx}$. The first case yields $x=O(\gamma^{-1}\epsilon^{2})$ and $Q=O(\gamma^{-1}\epsilon)$, which would eliminate $Q_{t}$ from the leading order equation. Hence, the only sensible balance is given by the two non-classical contributions, which gives $x=O(\gamma^{-1/2}\epsilon^{1/2})$ and $Q=O(\gamma^{-1/2}\epsilon^{-1/2})$. In terms of the new variables defined by $T=\hat{T}$, $x=\gamma^{-1/2}\epsilon^{1/2}\hat{x}$, $s=\hat{s}$, $t=\epsilon\hat{t}$ and $Q=\gamma^{-1/2}\epsilon^{-1/2}\hat{Q}$, the Stefan condition becomes ${\hat{s}}_{\hat{t}}=O(\gamma^{-1/2}\beta^{-1}\epsilon^{1/2})$, thus $\hat{s}\approx 1$, and the G-K equation becomes

[TABLE]

Since $\epsilon$ is arbitrarily small, the leading order problem becomes

[TABLE]

Equations (S23a) and (S23b) can be combined to give

[TABLE]

from where we recover Fourier’s law after integrating with respect to time and applying Eq. (S23e). The problem is therefore equivalent to a classical heat conduction problem in a semi-infinite domain with an effective thermal conductivity $\eta^{2}$. As suggested by the fact that the length scale has remained undetermined, this problem admits a solution in terms of a similarity variable $\hat{y}=\hat{x}/\sqrt{\hat{t}}$. In terms of the original non-dimensional variables, the solution takes the form

[TABLE]

where $\zeta=\eta^{2}/\gamma$.

SM4.1.2 Second time regime

Equation (S22) indicates that the previous time regime breaks down for $t=O(1)$, since the term $T_{x}$ enters the leading order problem. From the previous time regime we find $x,Q=O(\gamma^{-1/2})$, whereas $s,T=O(1)$. Upon defining the new variables $T=\tilde{T}$, $x=\gamma^{-1/2}\tilde{x}$, $s=\tilde{s}$, $t=\tilde{t}$ $q=\gamma^{-1/2}\tilde{q}$, we obtain again $\tilde{s}_{\tilde{t}}=O(\gamma^{-1/2}\beta^{-1})$, hence $\tilde{s}\approx 1$, and

[TABLE]

By eliminating the flux $\tilde{Q}$, we find that the temperature satisfies a linear viscoelastic equation given by

[TABLE]

with the additional condition that $\tilde{T}_{\tilde{t}}(\tilde{x},0)=0$. The large-time limit of (S27) implies that the solution for the temperature is given by a right-moving travelling wave that has an interior, moving boundary layer located at $\tilde{x}=\tilde{t}$ that captures weak non-local (“viscous”) effects. The solution for the temperature for $\tilde{x}=O(\tilde{t})$ with $\tilde{t}\gg 1$ can be written as

[TABLE]

where the upstream and downstream solutions $\tilde{T}_{\text{up}}$ and $\tilde{T}_{\text{down}}$ are, in principle, obtained by matching to the solution for $O(1)$ times and the inner solution in the boundary layer is

[TABLE]

This solution reveals that the temperature will remain $O(1)$ in size for large times. The large-time approximation to the flux can be obtained via (S26a) and shown to remain $O(1)$ in size for large times as well. The boundary layer at $\tilde{x}=\tilde{t}$ acts as the diffuse wavefront of the propagating thermal wave.

SM4.1.3 Third time regime

This time regime is analogous to that analysed in Sec. SM3.2 for the MCE, in which thermal waves propagate and collide with the boundaries of the crystal. The first collision occurs at the right boundary when $t=\gamma^{1/2}$. Therefore, in this time regime we choose a length scale of $x=O(1)$ and a time scale of $t=O(\gamma^{1/2})$. The heat flux has not changed in magnitude, therefore we choose a flux scale of $Q=O(\gamma^{-1/2})$ as in the previous time regime. Additionally, we have $s,T=O(1)$. In terms of the corresponding scaled variables $\bar{x},\bar{s},\bar{t},\bar{Q}$ and $\bar{T}$, the Stefan condition takes the form $\bar{s}_{\bar{t}}=O(\beta^{-1})$ and hence $\bar{s}\approx 1$, which leads to the problem

[TABLE]

where the flux has been eliminated from the GKE to produce (S30a). Naively taking $\gamma\to\infty$ in the GKHE (S30a) shows that the temperature satisfies the same leading-order problem as in Sec. SM3.2 and thus we can expect similar wave propagation and destructive interference to occur. The main difference is that non-local effects will influence the temperature and flux profiles near the wavefront. Remarkably, they do not influence the reflection dynamics to leading order.

We first consider the dynamics that occur before the thermal wave collides with the right boundary ($0<\bar{t}<1$). The leading-order solutions for the temperature and flux can be constructed using matched asymptotic expansions and are found to be given by $\bar{T}(\bar{x},\bar{t})=\mathcal{F}(\bar{x},\bar{t})$ and $\bar{Q}(\bar{x},\bar{t})=\mathcal{F}(\bar{x},\bar{t})$ where

[TABLE]

Equation (S31) is analogous to the Heaviside function that was found in Sec. SM3.2 but it now accounts for the diffuse nature of the wave front which has a finite width of $O(\gamma^{-1/4}\eta t^{1/2})$.

To understand the dynamics of the reflection that occurs at $\bar{x}=1$ when $\bar{t}=1$, we write $\bar{x}=1+\eta^{2}\gamma^{-1/2}\bar{x}^{\prime}$ and $\bar{t}=1+\eta^{2}\gamma^{-1/2}\bar{t}^{\prime}$ so that the temperature now satisfies $\bar{T}_{\bar{t}^{\prime}\bar{t^{\prime}}}=\bar{T}_{\bar{x}^{\prime}\bar{x}^{\prime}}+\bar{T}_{\bar{x}^{\prime}\bar{x}^{\prime}\bar{t^{\prime}}}$. The bulk equation admits wave solutions of the form $\bar{T}\sim\mathrm{e}^{i(k\bar{x}^{\prime}-\omega\bar{t}^{\prime})}$, where $\omega=\pm(ik/2)[k+(k^{2}-4)^{1/2}]$. Since $\omega$ is complex for wavenumbers that satisfy $k<2$, these short-wave modes are attenuated during wave propagation. By considering a single wave mode, writing the temperature as an incoming and reflected wave, $\bar{T}=\bar{T}_{I}\mathrm{e}^{i(k\bar{x}^{\prime}-\omega\bar{t}^{\prime})}+\bar{T}_{R}\mathrm{e}^{i(-kx-\omega\bar{t}^{\prime})}$, and applying the boundary condition $\bar{T}=0$ when $\bar{x}^{\prime}=0$, we find that $\bar{T}_{R}=-\bar{T}_{I}$. That is, the reflected wave has the same amplitude as the incoming wave but opposite sign, leading to destructive interference again. Moreover, non-local effects to not alter the amplitude of the reflected wave. By considering a single-mode wave solution for the flux, we find from conservation of energy that $\bar{Q}=(\omega/k)[\bar{T}_{I}\mathrm{e}^{i(k\bar{x}^{\prime}-\omega\bar{t}^{\prime})}+\bar{T}_{I}\mathrm{e}^{i(-k\bar{x}^{\prime}-\omega\bar{t}^{\prime})}]$. Therefore, the reflection of the flux wave leads to the same constructive interference as before.

Using the superposition principle, the solutions for the temperature and flux after the first reflection, but before the second reflection with the left boundary, can be written as

[TABLE]

for $1<\bar{t}<2$. In fact, since the second terms of (S32) are exponentially small for $\bar{t}<1$, these expressions are valid for $0<\bar{t}<2$. By further exploiting the asympototic limits of the function $\mathcal{F}$, the solution for times given by $0<\bar{t}<2N$, where $N$ is an integer, is given by

[TABLE]

Physically, the expression in (S33) can be interpreted as a collection of right- and left- moving travelling waves, each of which has a diffuse wavefront that broaden in time according to a $O(t^{1/2})$ scaling law. From these solutions we can deduce that $\bar{T}=O(1)$ and $\bar{Q}=O(\bar{t})$ as $\bar{t}\to\infty$.

We are now in a position to determine the next time regime using the asymptotic behaviour of the temperature and flux for large $\bar{t}$. Using the current scaling, the GKE can be written as $\bar{Q}_{\bar{t}}+\gamma^{-1/2}\bar{Q}=-\bar{T}_{\bar{x}}+\eta^{2}\gamma^{-1/2}\bar{Q}_{\bar{x}\bar{x}}$, indicating that the second term on the left-hand side becomes relevant when $\bar{t}=O(\gamma^{1/2})$. Alternatively, balancing both terms in the Stefan condition implies that solidification will occur when $\bar{t}=O(\beta^{1/2})$. The choice of the time scale will therefore depend on the relative size $\gamma$ to $\beta$.

SM4.1.4 Fourth time regime

We consider the distinguished limit where $\gamma^{-1/2}\beta^{1/2}=\alpha=O(1)$. We choose $t,Q=O(\gamma^{1/2}\beta^{1/2})$, whereas $x,s,T=O(1)$. In the new variables $\check{x},\check{s},\check{t},\check{Q}$ and $\check{T}$, the leading order problem takes the form

[TABLE]

with matching conditions provided by the previous time regime. This system is analogous to (S19) and hence the GKE has effectively reduced to the MCE, since non-local effects have become negligible. Interestingly, this has happened before the solidification process has begun.

SM4.2 Enhanced non-local effects

Here we assume $\beta,\eta^{2}\gg\gamma=O(1)$, therefore non-local effects are the dominant non-classical feature of heat transport. There are three time regimes to consider. The first time regime captures the initial heat flux through the seed crystal described by an effective form of Fourier’s law. In the second time regime, both memory and non-local terms enter the leading-order problem; however, unlike the case of enhanced memory effects where the temperature satisfied a viscoelastic wave equation, now the problem is quasi-steady and the temperature profile is linear. Solidification begins in the second time regime but eventually terminates due to the diminishing influence of the non-classical conduction mechanisms that drive phase change. In the third time regime, the classical contribution to the flux from temperature gradients becomes dominant, allowing solidification to restart.

SM4.2.1 First time regime

We consider heat flow across the entire seed crystal and take $x=O(1)$. Since $T=O(1)$, conservation of energy implies that $Q=O(t)$. The only sensible balance in Eq. (S1a) is obtained by balancing the non-classical contributions, which gives $t=O(\eta^{-2})$. Defining $T=\hat{T}$, $x=\hat{x}$, $s=\hat{s}$, $t=\eta^{-2}\hat{t}$ $Q=\eta^{2}\hat{Q}$ we find $\hat{s}_{\hat{t}}=O(\beta^{-1})$ and thus solidification is negligible at leading order. Upon neglecting terms of order $\eta^{-4}$, the remaining equations become

[TABLE]

We can combine (S35b) and (S35a) to obtain (after neglecting the $O(\eta^{-2})$ term)

[TABLE]

from where we recover Fourier’s law. The problem is therefore equivalent to a classical heat conduction problem in the unit interval and with an effective thermal conductivity $\gamma^{-1}$. The large-time solution for the temperature is $\hat{T}\sim\hat{x}-1$, whereas $\hat{Q}\sim-\gamma^{-1}$.

There are two possible choices for the next time regime. On one hand, for $\hat{t}=O(\beta)$ the terms in the Stefan condition balance and the solidification enters the leading order problem. On the other hand, a new balance in (S35a) is obtained for $\hat{t}=O(\eta^{2})$. The choice of the next time scale therefore depends on the relative size of $\eta^{2}$ to $\beta$. However, notice that (S35b) reduces to $\hat{Q}_{\hat{x}}\sim 0$ for either choice and hence the pseudo-steady state is reached before solidification begins.

SM4.2.2 Second time regime

As discussed previously, the choice of the time scale depends on the relative size of $\eta^{2}$ to $\beta$. Let us consider the distinguished limit $\beta\eta^{-2}=:\alpha=O(1)$ and choose $t=O(\alpha)$. From the previous time regime we have $Q=O(\eta^{2})$ and $x,s,T=O(1)$. In the corresponding scaled variables $T=\tilde{T}$, $x=\tilde{x}$, $s=\tilde{s}$, $t=\alpha\tilde{t}$ and $Q=\eta^{2}\tilde{Q}$, the Stefan condition becomes

[TABLE]

and, upon neglecting terms of order $\eta^{-2}$, the remaining equations take the form

[TABLE]

Taking derivatives with respect to $\tilde{x}$ in (S38a) yields ${\tilde{T}}_{\tilde{x}\tilde{x}\tilde{t}}=0$. Additionally, from the previous time regime we have the matching condition ${\tilde{T}}_{\tilde{x}\tilde{x}}\sim 0$ for $\tilde{t}=O(\beta^{-1})$, therefore ${\tilde{T}}_{\tilde{x}\tilde{x}}=0$ and hence, applying the boundary conditions,

[TABLE]

which reduces (S38a) to

[TABLE]

Using (S37) and the fact that ${\tilde{Q}}_{\tilde{t}}=(\partial\tilde{Q}/\partial\tilde{s}){\tilde{s}}_{\tilde{t}}$ we find that $\tilde{Q}$ is given by

[TABLE]

where we have applied the matching condition $\tilde{Q}\sim-\gamma^{-1}$ and $\tilde{s}\sim 1$ for $\tilde{t}=O(\beta^{-1})$. Replacing $\tilde{Q}$ in (S37) by this expression yields

[TABLE]

In particular, we observe that the initial solidification kinetics are similar to those predicted by the classical formulation, since $\tilde{s}\sim 1$ for $\tilde{t}\ll 1$ and hence $\tilde{s}_{\tilde{t}}\sim 1/(\gamma\tilde{s})$. In the original dimensionless variables we have $s_{t}=O(\eta^{2})$, which indicates that non-local effects accelerate the solidification process. Interestingly, as the solid grows into the liquid, we find $\tilde{s}_{\tilde{t}}\to 0$ as $\tilde{s}\to(1+\sqrt{1+4\alpha^{-1}})/2=\tilde{s}_{0}$ and therefore $\tilde{Q}\to 0$ as the interface approaches this value. Hence, the problem becomes stationary at the end of this regime and the interface tends to the value $\tilde{s}_{0}=O(1)$. However, as the flux decreases and the time increases, the term $T_{x}$ enters the leading order problem and solidification starts again. The only possible balance for $\tilde{Q}\ll 1$ and $\tilde{t}\gg 1$ is when $\tilde{Q}=O(\eta^{-2})$ and $\tilde{t}=O(\eta^{2})$, which also keeps the terms in the Stefan condition balanced, allowing solidification to start again as we enter the next regime.

SM4.2.3 Third time regime

From the previous time regime we have $t=O(\beta)$ and $s,x,T,Q=O(1)$. Introducing the scaled variables $\bar{t},\bar{x},\bar{s},\bar{T},\bar{Q}$, the equations become

[TABLE]

Upon differentiating (S43a) with respect to $\bar{x}$ gives $\alpha\bar{T}_{\bar{x}\bar{x}}+\bar{T}_{\bar{x}\bar{x}\bar{t}}=0$. Since the temperature profile $\bar{T}$ is linear in the previous time regime we find $\bar{T}_{\bar{x}\bar{x}}=0$ and hence

[TABLE]

Thus, the flux is given by

[TABLE]

and hence $\bar{s}$ is determined by

[TABLE]

In particular, we find $\bar{s}_{\bar{t}}\sim\bar{s}^{-1}$ for $\bar{t}\gg 1$, recovering so the classical solidification kinetics by the end of this time regime. The solution to (S46) can be given implicitly,

[TABLE]

or in terms of the LambertW function,

[TABLE]

The value of $C$ is found by imposing $\bar{s}\sim\tilde{s}_{0}$ for $\bar{t}=O(\eta^{-2})$ we find that $\bar{C}=\tilde{s}_{0}^{2}/2+\alpha^{-1}\log\tilde{s}_{0}$.

As suggested by the large-time behaviour of (S46), the classical dynamics are recovered in the fourth time regime, given by $\bar{x},\bar{s}\sim\bar{t}^{1/2}$, $\bar{Q}\sim\bar{t}^{-1/2}$ and $\bar{t}=O(\eta^{2})$.