A slip-based model for the size-dependent effective thermal conductivity of nanowires

The heat flux across a nanowire is computed based on the Guyer-Krumhansl equation. Slip conditions with a slip length depending on both temperature and nanowire radius are introduced at the outer boundary. An explicit expression for the effective thermal conductivity is derived and compared to existing models across a given temperature range, providing excellent agreement with experimental data for Si nanowires.


Introduction
Nanotechnology is currently the focus of extensive research due to its wide range of applications in fields such as industry and medicine [1][2][3][4][5][6][7]. Nanowires, in particular, are being used in technologies relating to solar cells [4], flexible screens [5], detection of cancerous cells [7], and energy storage [6]. A key issue facing the practical use of nanodevices is thermal management [8]. Inefficient regulation of heat can lead to large temperatures and melting, possibly resulting in device failure. Understanding and predicting heat flow on the nanoscale is therefore crucial for the manufacturing and operation of nanotechnologies.
It is widely known that many thermophysical material properties become size-dependent at the nanoscale [9][10][11][12][13][14][15][16]. Buffat and Borel [9] showed a dramatic decrease of the melting temperature of gold nanoparticles of almost 50% from the bulk value. For aluminium nanoparticles, a decrease in latent heat by a factor of four has been reported [10]. Experimental observations also demonstrate that the thermal conductivity in silicon nanowires is much lower than the theoretical value predicted by kinetic theory [14]. For instance, it is reported that, at room temperature, the thermal conductivity of Si nanowires with a diameter of 37 nm decreases by approximately 87% with respect to the bulk value. Furthermore, when the characteristic size of the system is much smaller than the phonon mean free path, the thermal conductivity shows an approximately linear dependence on size [17][18][19][20].
The size dependence of the thermal conductivity of nanosystems is attributed to the fact that, on the nanoscale, the transport of thermal energy is a ballistic process driven by infrequent collisions between thermal energy carriers known as phonons. This is in contrast to macroscopic heat transfer, which is a diffusive processes driven by frequent phonon collisions. As the size of a device becomes commensurate with the phonon mean free path, bulk phonons are more likely to collide with a boundary than with each other. The ability of a nanodevice to conduct thermal energy, therefore, becomes strongly influenced by the scattering dynamics at the boundary as well as the geometrical structure (e.g., size and shape) of this boundary.
Due to the fundamentally different manner in which heat is transported across nanometer length scales in comparison to heat flow at the macroscale, Fourier's law is unable to provide an accurate description of heat conduction in this regime [21]. Different approaches to modelling nanoscale heat flow have been developed in order to capture the ballistic nature of energy transport and size dependence of the effective thermal conductivity (ETC). These approaches can be classified into three main categories: microscopic, mesoscopic and macroscopic models. Microscopic approaches, such as molecular dynamics or Monte-Carlo methods [22], focus on the evolution of every single phonon while mesoscopic models group them together depending on their wavelength and wavevector. Micro and mesoscopic models are mainly based on the Boltzmann transport equation (BTE) and its solution under different approximations. A popular example is the equation of phonon radiative transfer (EPRT) [23], where an expression for the ETC similar to the classical expression from kinetic theory is derived, although here an effective mean free path is now considered. However this model is based on the gray approximation and thus considers a single phonon group velocity and lifetime. A more general model where these quantities are mode-dependent was presented by McGaughey et al. [24]. Starting again from the BTE, Alvarez et al. [18,25] extract a continued-fraction expression to describe the ETC in thin films. Other models, such as those of Callaway [26] and Holland [27], also consider phonon distributions rather than single phonons, but Mingo et al. [28] showed that they fail when predicting the ETC for Si or Ge nanowires.
Macroscopic models aim to describe global variables of the system, such as the temperature and the heat flux. A recent approach at the macroscopic scale is the thermomass model, where heat carriers are assumed to have a finite mass determined by Einstein's mass-energy relation [19,29,30]. Other approaches are based on the Guyer-Krumhansl (G-K) equation [31,32], which is derived from a linearized BTE in dielectric crystals. This equation has become popular since it is analogous to the extensively studied Navier-Stokes (N-S) equations and it is one of the simplest extensions to Fourier's law that includes memory and non-local effects. Models based on the G-K equation are included in the framework of phonon hydrodynamics [33]. For instance, Alvarez et al. [20] use the analogy between the G-K equation and the N-S equations to derive an expression for the ETC in circular nanowires, splitting the heat flux into two separate contributions. This work has been extended to elliptical and rectangular nanowires [34]. However, since the size of the phonon mean free path depends on temperature, the assumptions on which they base their reductions are only valid for low temperatures or very small sizes. Dong et al. [35] find an expression for the ETC by solving the full G-K equations at steady state, although their no-slip boundary condition leads to a quadratic dependence of the ETC on the characteristic size of the device for large Knudsen numbers instead of the known linear behaviour. This does not happen in the case of the expression derived by Ma [17], where a fixed flux is imposed on the outer boundary. However, Ma's solution for the nanowire shows a very poor match to data (their paper contains an error in that they plot the thin film solution in their figure for the nanowire; this actually shows reasonable agreement with the data). Further, most of the existing models are only validated at room temperature and therefore a deeper assessment of their accuracy is required.
In this paper, we introduce a new phenomenological slip boundary condition that, when used with the G-K equations, results in ETC predictions that are in excellent agreement with experimental data for Si nanowires over a range of radii and temperatures. The proposed Figure 1: A circular nanowire with radius R * and length L * is held at different temperatures T * 0 , T * 1 at the left and right ends respectively, which induces a heat flux Q * .
model is remarkably simple and only requires knowledge of the temperature dependence of the bulk thermal conductivity and phonon mean free path, both of which may be obtained experimentally or computationally, making it well suited for use in practical applications. A detailed comparison of the proposed and existing models is performed, the results of which show that the proposed model consistently yields the most accurate predictions of the ETC compared to existing models. Furthermore, this comparison establishes the validity of each model in terms of temperature and nanowire radius.

Mathematical modelling
We consider a circular nanowire (NW) of radius R * and length L * that is suspended in a vacuum; see Fig. 1. The radius of the NW is assumed to be much smaller than its length, i.e., R * /L * ≡ 1. A temperature gradient ∆T = T * 0 − T * 1 > 0 is imposed along the axial direction of the NW by fixing the temperature at its left and right ends to be T 0 and T 1 , respectively. The thermal flux that is driven by this temperature gradient is assumed to be axisymmetric. Therefore, it is sufficient to consider a two-dimensional model with radial and axial coordinates r * and x * , respectively. The mathematical model will consist of an equation representing conservation of thermal energy and the G-K equation describing the evolution of the thermal flux.

Conservation of energy requires
where u * (T * ) is the internal energy per unit mass and T * is the temperature. The thermal flux Q * = v * r + w * x is assumed to satisfy the G-K equation where v * and w * are the radial and axial components of the heat flux, τ * (T * ) is the total mean free time, k * (T * ) is the bulk thermal conductivity in the Kinetic Collective Model (KCM) framework [36], and * (T * ) is a non-local length related to the bulk phonon mean free path (MFP), i.e., the mean distance between phonon-phonon collisions. Zhu et al. [37] propose that * is the geometric mean of the bulk MFP and a local MFP, the latter of which decreases near a boundary. Here, we opt for simplicity and model the decrease in MFP near the boundary through the slip boundary condition. For convenience, we will not write the temperature dependence of the parameters explicitly unless it is required due to the context. However, we note that, for the temperature ranges considered here, both the bulk thermal conductivity and non-local length monotonically decrease with temperature. As will be shown in Sec. 4, the non-local length of silicon decreases from about 5 µm at 50 K to 55 nm at 300 K. The first term on the left-hand side of (2) captures memory effects and, in particular, the dependence of heat flux on the history of the temperature gradient. The second term on the right-hand side of (2) captures non-local effects, such as the interaction of phonons with the boundary of the NW. When the characteristic time and length scales are much larger than the resistive mean free time and non-local length, the G-K equation (2) reduces to Fourier's law.
For the remainder of the paper, we restrict our attention to the case of steady-state heat flow. This focus is motivated by the available experimental data. Under the steady-state assumption, ∂u * /∂t = ∂Q * /∂t = 0. Conservation of energy (1) and the G-K equation (2) then reduce to This system is clearly analogous to an incompressible, viscous flow with a source term proportional to the velocity. Under this analogy, the parameter µ * = * 2 /k * may be interpreted as a thermal viscosity. These observations allow us to use well-known techniques from viscous flow to analyse the problem.

Boundary conditions
The boundary conditions for the temperature at the endpoints of the NW are where T * 0 > T * 1 . Due to the symmetry of the problem the boundary conditions at r * = 0 are straightforward, We assume that no heat flows across the outer boundary, i.e., v * | r * =R * = 0.
Finally, on the edge of the NW we continue the analogy with fluid dynamics and employ a slip condition, where * s is the slip length. In accordance with previous authors (see, e.g., Refs. [20,37]), we will assume that the slip length is proportional to the non-local length, i.e., * s = C * , where C is a dimensionless parameter that encodes detailed information about phonon-boundary scattering and surface roughness. This form of boundary condition is discussed in more detail in Refs. [20,38]. The basic idea is that a slip condition can capture crucial contributions to the thermal flux from reflected phonons. A no-slip condition neglects such contributions, resulting in predictions of the thermal flux and hence ETC that decrease too rapidly with the radius of the NW.
Various forms for the parameter C appear in the literature. Alvarez et al. [20] write C in terms of the specular parameter p, which describes the precise nature of phonon scattering (i.e., specular or diffuse). A similar approach is presented in Zhu et al. [37], although the authors essentially treat p, and hence C, as a fitting parameter. Sellitto et al. [38] write C in terms of the surface geometry and then expand it as a power series in the temperature, resulting in a model with an excessive number of fitting parameters. Here, we take a more practical approach and write C as a simple exponential, C = exp(−R * / * ), which is similar in form to the expression derived for and used in the specific case of a two-dimensional rectangular nanolayer by Zhu et al. [37]. The proposed exponential form does not involve any fitting parameters and it inherits a temperature dependence through the bulk non-local length * (T * ). The motivation for this expression is as follows. If the non-local length is much less than the NW radius, * R * , then phonons are more likely to collide with each other than with the boundary and the influence of reflected phonons will be small, corresponding to a no-slip condition. Conversely, if the nonlocal length is much larger than the NW radius, * R * , then phonon-boundary scattering is more likely than phonon-phonon collisions. In this case, phonons reflections will strongly influence the flow, which is captured in a slip condition.

Reduction of the governing equations
We may reduce the governing equations by exploiting the separation between the axial and radial length scales. First, we introduce non-dimensional variables where v * 0 and w * 0 are (unknown) typical values of the axial and radial components of the heat flux and k * 0 is a reference value for the bulk thermal conductivity. A relationship between v * 0 and w * 0 is obtained by requiring both terms in (3a) to have the same magnitude, which ensures that energy is conserved in the leading-order problem. This yields v * 0 = w * 0 . Finally, the classical scale for the axial flux is chosen, w * 0 = k * 0 ∆T /L * , which ensures that Fourier's law is recovered in the classical limit, i.e., when the Knudsen number Kn = * /R * tends to zero.
After non-dimensionalising the model given by (3) The boundary conditions become

Calculation and interpretation of the effective thermal conductivity
The ETC is defined as the ratio of the heat flux per unit area to the temperature gradient driving this flux. In non-dimensional form, the ETC may be expressed as where the heat flux per unit area is Calculating the integral in (13) requires solving Eqn. (9b) for the axial component of the flux w, which is trivial because the temperature does not depend on the radial coordinate. Upon solving (9b) and applying the boundary conditions (11b), we find that where I ν is the modified Bessel function of the first kind of order ν. The heat flux per unit area may then be calculated as in (13), which yields According to its definition in (12), the ETC is finally given by In principle, the temperature can be calculated from (15) by first noting that, at steady state, the flux Q must be uniform in the axial direction: if a different amount of energy enters the wire to that leaving, then the temperature must vary with time. Due to the temperature dependence of the parameters appearing in (15), we have that T satisfies a nonlinear ODE of the form which is analogous to the Reynolds equation in fluid mechanics [39]. Solving the first-order ODE (17) and imposing the boundary condition T = 1 at x = 0 determines the temperature in terms of the flux Q. Subsequently applying the boundary condition T = 0 at x = 1 to the solution for the temperature enables the flux to be obtained. Due to the dependence of Kn = * (T * )/R * on temperature and the radius of the NW, Eqn. (16) can be reduced to simpler expressions in some limiting cases. For instance, at very low temperatures or for very small radii we have Kn 1. Using I 0 (ξ) = 1 + O(ξ 2 ) and I 1 (ξ) = ξ/2 + O(ξ 3 ) for ξ 1, as well as C(Kn) = 1 − Kn −1 + O(Kn −2 ) for Kn 1, we find that (16) can then be reduced to from which we deduce k * eff ∝ R * , in agreement with previous theoretical results [17][18][19][20]. In deriving the leading-order term in (18), only the leading contribution to C, given by C ∼ 1, is required. Thus, any alternative expression of C that has a large-Kn expansion of the form C = 1 + O(Kn −α ), with α > 0, will result in the asymptotic behaviour k * eff /k * ∼ (2Kn) −1 as Kn → ∞. For larger devices or at high temperatures, where Kn 1, the ETC should reduce to its classical value as non-local effects become negligible. To show this, we use the relation [40] hence the conductivity tends to the bulk value as Kn → 0. The asymptotic behaviour k * eff /k * ∼ 1 − 2Kn as Kn → 0 will be true for all functional forms of C with the limit C → 0 as Kn → 0.
To aid with the physical interpretation of the model solutions and their asymptotic limits, the radial dependence of the axial flux w is plotted for three Knudsen numbers (Kn = 0.05, 1, and 10) in Fig. 2 (a). The ETC, given by (16), is plotted as a function of the Knudsen number in Fig. 2 (b) for two choices of the function C. The figures make it clear that there are three distinct regimes to consider, corresponding to diffusive (Kn 1), ballistic (Kn 1), and mixed (Kn = O(1)) modes of thermal energy transport. The three values of the Knudsen number in Fig. 2 (a) are chosen to clearly illustrate how the axial flux varies across these three regimes.
For small Knudsen numbers, Kn 1, the phonon MFP is small compared to the radius of the NW. Phonons in the bulk are therefore more likely to collide with each other before reaching the boundary of the NW and scattering. This corresponds to the boundary condition w → 0 as Kn → 0. Due to the relatively small influence of non-local effects, thermal energy in the bulk is transported with little resistance across the NW, with an axial flux that approximates the classical flux w F predicted by Fourier's law. However, there is a thin boundary layer near the edge of the NW where the flux rapidly decreases to zero in order to satisfy the no-slip condition. The (dimensional) width of this boundary layer is O( * ), reflecting the fact that boundary effects become relevant on length scales that are commensurate with the phonon MFP. As the Knudsen number tends to zero, the width of this boundary layer does so as well. The restricted transport of thermal energy in the boundary layer leads to a slight reduction in the ETC compared to its bulk value, which is captured in the small-Kn limit of k eff given by (19) and shown in Fig. 2 (b) by the slow reduction in k * eff /k * as Kn increases from 10 −3 to 10 −1 .
In the limit of large Knudsen number, Kn 1, the phonon MFP greatly exceeds the radius of the NW. Phonons therefore collide more frequently with the boundary than with each other, resulting in a substantial decrease in the flow of thermal energy. The strong influence of nonlocal effects and large slip length lead to an axial flux that is uniform along the radial direction and, in particular, close to its value on the boundary, w ∼ Kn −1 w F /2. Therefore, as Kn → ∞, corresponding to an arbitrarily large MFP relative to the NW radius, the axial flux tends to zero. The restricted transport of energy in this regime results in a marked decrease in the ETC, the magnitude of which also scales linearly with Kn −1 , as demonstrated in Eqn. (18).
The proportionality between the flux w and the inverse Knudsen number Kn −1 can be understood by drawing on the analogy between nanoscale heat transfer and fluid dynamics. In the limit of large Kn (analogous to large viscosity), the net force exerted by the temperature gradient (analogous to pressure) on an arbitrary cross section of the NW must balance the net traction (i.e., shear stress) along its circumference, From Eqn. (20) it is seen that the rate of shear strain, ∂w/∂r, induced by the temperature gradient is small and of size O(Kn −2 ). The scaling for the flux (analogous to velocity) is determined by the slip condition (11b), yielding w ∼ Kn ∂w/∂r ∼ Kn −1 w F /2, where the rate of strain has been eliminated using (20) and C has been approximated by its leading-order contribution C ∼ 1 as Kn 1. Thus, the scaling with Kn  Figure  2 (b) shows the ETC when C = exp(−1/Kn) (solid line) and C = 1 (dashed line). The former corresponds to the situation where the slip length * s is smaller than the non-local length * (recall that * s = * C). Taking C = 1 implies that the slip length and MFP are equal. The exponential form of C results in a slight kink near Kn = 1. In physical terms, this kink is associated with a reduced ETC compared to the C = 1 case and is the consequence of a smaller slip length and hence axial flux. Using C = 1 smooths out this kink and yields higher values of the ETC due to greater slip.

Model validation and comparison
The prediction of the ETC is now compared against experimental data measured from Si nanowires with diameters D * = 2R * = 37 nm, 56 nm, and 115 nm at temperatures ranging from T * = 50 K up to 300 K [14]. Evaluating the ETC prediction requires knowledge of the temperature dependence of the non-local length * and thermal conductivity k * . These two quantities are obtained from first principles in the KCM framework [36] using an open-source code [41] and plotted as functions of temperature in Fig. 3. The bulk MFP ranges from 55 nm to 5 µm, corresponding to Knudsen numbers between 2.98 and 262 (D * = 37 nm), 1.97 and 173 (D * = 56 nm), and 0.958 and 84.2 (D * = 115 nm). Thus, the exponential form of C will be relevant when comparing the ETC prediction to the experimental data. Indeed, we find poor agreement when the the dependence of C on Kn is neglected and C is held at unity. Figure 4 shows the ETC predicted by (16) (solid lines) and experimentally measured values of the ETC (symbols) as functions of temperature for three NW diameters (Figs. 4 (a)-(c)) and as functions of the diameter for a fixed temperature (Fig. 4 (d)). The ETC k * eff is normalised against the bulk value k * to highlight the relative change that occurs as the diameter of the NW and non-local length decrease. In all cases, the theoretical predictions are in excellent agreement with the experimental data. The dashed, dotted, and dashed-dotted lines in Fig. 4 represent theoretical predictions of the ETC computed from existing models found in Refs. [17,18,35], respectively. It is clear that that the model proposed here consistently yields predictions that are more representative of the experimental data. The models of Ma [17] and Dong et al. [35], both of which are based on the G-K equation, are unable to quantitatively capture the experimental data associated with all three NW diameters. These inaccuracies can be attributed to the different boundary conditions that are used in their models. Dong et al. [35] employ a no-slip boundary condition, equivalent to C = 0 in our model. Imposing this boundary condition leads to particularly large errors at low temperatures (or large Knudsen numbers). Moreover, the calculated ETC becomes quadratically dependent on the NW radius, i.e., k * eff /k * ∝ R * 2 , which contradicts the experimentally measured linear dependence between these two quantities. Ma [17] instead considers a prescribed flux on the boundary, which leads to substantial inaccuracies across the whole temperature range. The expression for the ETC presented by Alvarez and Jou [18] gives the most accurate prediction of the three existing models, although it generally underestimates the experimentally measured values. A distinguishing feature of Alvarez and Jou's expression is that it is derived from the linearised Boltzmann equation rather than the G-K equation. This derivation does not explicitly account for the cylindrical geometry of the nanowire and thus no boundary conditions can be imposed along the outer edge, which may contribute to the errors seen in Fig. 4. In comparison to the models derived from the G-K equation, the model of Alvarez and Jou predicts that the ETC converges more rapidly to the macroscopic value with increasing NW diameter; see Fig. 4 (d). However, the lack of experimental data for microscale NW diameters prevents the models from being validated in this regime.

Conclusion
We have proposed a new phenomenological slip model that can be coupled with the Guyer-Krumhansl equation in order to obtain an analytical expression for the size-dependent ETC in thin nanowires. A key advantage of the proposed model is that it only requires knowledge    [17,18,35], respectively; and circles denote experimental data [14]. of the bulk thermal conductivity and MFP and does not involve any fitting parameters. The prediction of the ETC of Si nanowires is found to be in excellent agreement with experimental measurements, which is particularly remarkable given the simple nature of the underlying model. When compared against existing models for the ETC, we observed that the proposed model yielded the most accurate predictions in all cases.
Crucial to the success of our model is the condition imposed on axial flux w at the outer boundary of the nanowire. Similar models, but with different boundary conditions, were shown to be in poor agreement with much of the experimental data. It is surprising that the present choice, based purely on an analogy with viscous fluid flow, works so well without resorting to any parameter fitting. Consequently, in our future work, we intend to further investigate the precise form of the boundary condition and, ideally, produce an expression based on a detailed physical model.
Finally, although we focused on circular nanowires, the model could, in principle, be applied to arbitrary geometries. Analytical solutions could be sought in the case of simple geometries such as rectangular nanofilms. In the case of complex geometry, the model could be straightforwardly implemented in finite element software, enabling a broad range of scenarios to be explored.