## Abstract

Wide applications of nanofilms in electronics necessitate an in-depth understanding of nanoscale thermal transport, which significantly deviates from Fourier's law. Great efforts have focused on the effective thermal conductivity under temperature difference, while it is still ambiguous whether the diffusion equation with an effective thermal conductivity can accurately characterize the nanoscale thermal transport with internal heating. In this work, transient in-plane thermal transport in nanofilms with internal heating is studied via Monte Carlo (MC) simulations in comparison to the heat diffusion model and mechanism analyses using Fourier transform. Phonon-boundary scattering leads to larger temperature rise and slower thermal response rate when compared with the heat diffusion model based on Fourier's law. The MC simulations are also compared with the diffusion model with effective thermal conductivity. In the first case of continuous internal heating, the diffusion model with effective thermal conductivity under-predicts the temperature rise by the MC simulations at the initial heating stage, while the deviation between them gradually decreases and vanishes with time. By contrast, for the one-pulse internal heating case, the diffusion model with effective thermal conductivity under-predicts both the peak temperature rise and the cooling rate, so the deviation can always exist.

## 1. Introduction

Wide applications of semiconductor nanofilms in electronics, e.g. silicon on insulator device, necessitate an in-depth understanding of nanoscale thermal transport [1]. Electrons (or electron holes) scatter with lattice vibrations (phonons) and lose net energy to phonons, consequently heating up the structures via the mechanism known as Joule heating which can greatly lower the efficiency of electronics. Phonons predominate the heat conduction in semiconductors [2]. In nanofilms whose thickness is comparable with the phonon mean free path, owing to phonon-boundary scattering, heat conduction deviates from Fourier's law which neglects the boundary effect. The phonon Boltzmann transport equation is usually applied to characterize the nanoscale thermal transport [3],
*v*_{g} is the group velocity, *f* is the phonon distribution function, *f*_{0} is the equilibrium distribution function, *τ* is the relaxation time and

Various methods have been used to investigate the in-plane thermal transport in nanofilms with internal heating [4–8]. Through solving the phonon Boltzmann transport equation, Sverdrup *et al.* [4] found that the peak temperature rise in the nanofilms could be significantly larger than the prediction using the heat diffusion equation based on Fourier's law. Narumanchi *et al.* [5] discussed the influence of the range and duration of heat generation on the temperature distribution evolutions within the nanofilms. Pop *et al.* [6] analysed the heat generation and transport in nanoscale transistors using the Monte Carlo (MC) simulations. Wong *et al.* [7] studied the effect of ballistic phonon transport in the nanofilms with heat generation via the MC simulations. Besides, the molecular dynamics simulation was also used to characterize the transient thermal phenomena of heat conduction in the silicon nanofilms with internal heating [8]. In general, it has been found that the phonon-boundary scattering significantly suppresses the in-plane thermal transport in the nanofilms with internal heating, resulting in larger temperature increase and slower thermal response rate when compared with the heat diffusion equation based on Fourier's law.

Great efforts have focused on the reduction of in-plane effective thermal conductivity in the nanofilms [9–23] and tried to characterize the thermal transport with phonon-boundary scattering by the modification of effective thermal conductivity [24–26]. Analytical effective thermal conductivity models were derived based on the steady-state phonon Boltzmann transport equation [9,10], extended irreversible thermodynamics [11] and phonon hydrodynamics [12,13]. Besides, the molecular dynamics [14] and MC [15] simulations were also employed to investigate the effective thermal conductivity. Sellitto *et al.* [16] investigated the steady-state radical heat transport in silicon thin layers based on phonon hydrodynamics and found that two-dimensional radial situations cannot be reducible to Fourier law with the effective thermal conductivity. As for experiments, Liu & Asheghi [17,18] obtained the in-plane effective thermal conductivity of silicon films by introducing a steady-state uniform internal heat source and measuring the average temperature increase. In the work of Johnson *et al.* [19], a transient internal heat source was introduced via diffraction of a laser beam and the in-plane thermal conductivity of the silicon films was measured by fitting the temperature attenuation rate. Besides, the thermoreflectance and 3-omega techniques are other powerful tools to measure the effective thermal conductivity of nanostructures in the transient heating state [20–23]. Although researchers have concentrated on investigating the size-dependence of in-plane effective thermal conductivity, an essential question is still ambiguous, namely whether the combination of the heat diffusion model and effective thermal conductivity can accurately characterize the in-plane thermal transport in nanofilms with internal heating, especially for transient states. This question does matter, the applications of effective thermal conductivity in electronics thermal design and management.

In this study, the transient in-plane thermal transport in the semiconductor nanofilms with internal heating was studied using the MC simulations. Two different cases were simulated: continuous heating and one-pulse heating. The heat diffusion model with effective thermal conductivity has also been solved for comparing with the MC simulations. Furthermore, the heat diffusion model with effective thermal conductivity has been analysed based on Fourier transform.

## 2. Simulation details and problem formulation

### (a) Monte Carlo simulation

A MC technique [27–33] is applied to simulate the phonon transport in nanofilms. The MC technique simulates phonon transport processes by random number samplings, equivalent to directly solving the phonon Boltzmann transport equation. In [27–30], a large ensemble of phonon bundles is simulated at the same time, while the tracing processes in [31–33] are independent for each single phonon bundle. The second tracing algorithm is adopted in our simulations due to it being less computationally expensive [33]. The phonon tracing algorithm is as follows: (1) Draw the initial properties (such as position *r*_{0} and time *t*_{0}) of phonons. (2) Calculate the travelling length, Δ*r*, until the first scattering event (*l*_{0} is the mean free path and *R* is the uniformly random number ranging from 0 to 1), and renew the position and time of phonons, *r*_{new}=*r*_{0}+Δ*r*, *t*_{new}=*t*_{0}+Δ*r*/*v*_{g}. (3) If a phonon collides with a boundary at *r*_{b}, set *r*_{new}=*r*_{b} and *t*_{new}=*t*_{0}+∥*r*_{b}−*r*_{0}∥/∥*v*_{g}∥, according to the nature of the boundary, determine whether the phonon re-emits from the boundary [31]. (4) If a phonon does not collide with boundaries, the phonon should re-emit at *r*_{new}, and set *r*_{0}=*r*_{new}, *t*_{0}=*t*_{new}. (5) If *t*_{new}>*t*_{final}, the tracing process of this phonon is finished, then we proceed to (1) and begin the tracing of the next phonon.

The bulk phonon properties of silicon at room temperature are chosen in our simulations. For clarity, we adopt the grey media approximation, which assumes that the phonon properties are frequency-independent. Therefore, phonons travel with one group velocity and the scattering rate is described by single phonon mean free path. The phonon mean free path is calculated as *l*_{0}=3*k*_{bulk}/(*ρc*_{V}*v*_{g}), where *k*_{bulk} is the bulk thermal conductivity, *c*_{V} is the specific heat, *ρ* is the mass density and *v*_{g} is the average group velocity. As for silicon at room temperature, *k*_{bulk} is 150 W (m K)^{−1}, *c*_{V} is about 700 J (kg K)^{−1}, *ρ* is 2330 kg m^{−3} and *v*_{g} is 6400 m s^{−1} [34]. Thus, the phonon mean free path *l*_{0} is about 43 nm. Arguments still exist on the phonon mean free path of silicon at room temperature [34]. Based on a more detailed dispersion model and only considering acoustic phonons that carry most of heat, the phonon mean free path of silicon ranges from 200 to 300 nm [34,35]. However, as the longer mean free path is chosen, the corresponding specific heat and group velocity should also be chosen. Here, the mean free path is chosen as 43 nm with its corresponding specific heat and group velocity. In fact, because all the quantities in the MC simulations and the models are dimensionless, the choice of mean free path does not influence the comparisons between them.

### (b) Problem formulation

The in-plane thermal transport in nanofilms with internal heating can be characterized by the two-dimensional Boltzmann transport equation, *x*-directional thickness is much longer than the *y*-directional thickness, i.e. *L*_{x}≫*L*_{y}, the phonon transport is mainly influenced by the *y*-directional boundary constraints. The *x*-directional boundaries are considered as phonon blackbody, that is, phonons are completely absorbed at them, while the lateral boundaries are adiabatic and the phonon-boundary scattering at them is assumed to be completely diffusive. The internal heat source is implemented by setting the phonon source inside the nanofilms, that is, phonons with the definite energy originating within the nanofilms. Then, phonons are traced in the domains until they exit through the *x*-directional boundaries, and the intrinsic scattering processes, such as phonon-impurities and phonon–phonon scatterings, are treated in the relaxation-time approximation. The temperature profile can be obtained from the phonon distribution calculated by the MC simulations [33]. Moreover, it should be noted that when the system is comparable with the phonon mean free path, the local thermodynamic equilibrium assumption cannot be achieved, and the temperature will lose its conventional meaning of representing a thermal equilibrium state. In this case, the temperature calculated in the MC simulations is a representation of the average energy of all phonons around a local point [34]. The tracing number of phonon bundles is equal to 10^{8}. The unit control volume is Δ*x*=0.02*L*_{x}, and the time step is 40 *τ*, in which *τ*=*l*_{0}/*v*_{g}. Two heating schemes are studied: continuous internal heating and one-pulse internal heating. In the continuous heating case, the sudden internal heating is imposed in the nanofilms at the beginning and then lasts all the time. Thus, the temperature will increase and reach the steady state at last. By contrast, as for the one-pulse internal heating case, the internal heating just lasts for a limited duration. Therefore, the temperature rise will increase in the heating duration and then decrease and vanish with time. In our simulations, the heating durations, Δ*t*, are set as 40 *τ*, 400 *τ* and 4000 *τ*. It should be pointed out that the heat wave effect does not need to be taken into account in this work, because both the heating and observation times are far beyond the relaxation time, *τ* [36–42]. In fact, when the heating duration is comparable with the relaxation time, i.e. Δ*t*=*τ*, due to the heat wave effect, the diffusion equation is not sufficient to derive an effective thermal conductivity, and the Maxwell–Cattaneo (or Cattaneo–Vernotte) equation incorporated with a size-dependent effective thermal conductivity may be applied to characterize the boundary and heat wave effects simultaneously [11].

## 3. Results and discussions

### (a) In-plane effective thermal conductivity of nanofilms

In the modelling work about the in-plane effective thermal conductivity of the nanofilms, a nanofilm is generally assumed to be in contact with two phonon baths of different temperatures and the temperature difference induces the heat flow [9–13]. Then the effective thermal conductivity can be calculated using Fourier's law. Based on the temperature difference scheme stated above, the theoretical model was derived from the phonon Boltzmann transport equation [9],
*Kn*_{y} is the Knudsen number defined as *Kn*_{y}=*l*_{0}/*L*_{y}, where *L*_{y} is the thickness of nanofilm. Equation (3.1) can be approximately simplified in the form of Matthiessen's rule, *k*_{eff_TD}/*k*_{bulk}≈1/(1+3/8*Kn*_{y}) [43]. In the internal heating case, researchers have been concerned more about the temperature rise induced by the internal heat source [25,26]. The effective thermal conductivity in this case should be obtained by comparing the actual temperature rise with the analytical solution of the diffusive heat conduction equation based on Fourier's law. Li & Cao [44,45] studied the effective thermal conductivity of the nanostructures with internal heat source by the non-equilibrium molecular dynamics simulations, and found that the effective thermal conductivity with internal heat source was significantly lower than that with temperature difference. Based on the phonon Boltzmann transport equation and the MC simulations, Hua *et al.* [46] investigated the ballistic-diffusive heat conduction in the nanostructures with internal heat source, and the in-plane effective thermal conductivity model for the nanofilms with internal heat source was also given
*α*=0.65 which is larger than 3/8.

The effective thermal conductivity, *k*_{eff_IHS}, is derived from average temperature rise within the nanofilms with internal heating [46]. Therefore, when we discuss about the in-plane thermal transport with internal heating, *k*_{eff_IHS} may be the more appropriate choice. As a verification, the steady-state in-plane thermal transport with internal heating is simulated and compared with the diffusion model with *k*_{eff_TD} and *k*_{eff_IHS}, respectively. As shown in figure 2, the dimensionless temperature, *T**, is defined as *η*, is *x*/*L*_{x}. It is found that the temperature rise increases with the Knudsen number *Kn*_{y} increasing. As *Kn*_{y}=0.01, the phonon-boundary scattering has no significant influence on the temperature distributions and the MC simulation results can be consistent with the predictions based on Fourier's law, which confirms the validity of our simulations. The diffusion model with *k*_{eff_IHS} can predict the steady-state temperature distributions, while the model with *k*_{eff_TD} significantly under-predicts the MC simulation results since *k*_{eff_TD}>*k*_{eff_IHS}. Therefore, the effective thermal conductivity, *k*_{eff_IHS}, obtained in the internal heating case can be the better choice, and the heat diffusion model with *k*_{eff_IHS}

### (b) Transient thermal transport with continuous internal heating

Figure 3 illustrates the peak temperature varying with time in the continuous heating case. As *Kn*_{y}=0.01, the boundary effect can be nearly neglected, so that the simulation results are consistent with the prediction by Fourier's law. With the Knudsen number, *Kn*_{y}, increasing, the influence of the phonon-boundary scattering grows up, leading to the larger peak temperature and more time to reach steady state. Furthermore, it is also found that the diffusion model with the effective thermal conductivity under-predicts the peak temperature predicted by the MC simulations at the initial heating stage, and the deviation between the model and the MC simulation decreases and vanishes with time.

The temperature distributions within the nanofilms are shown in figure 4 for different Knudsen numbers, *Kn*_{y} (the ratio between the phonon mean free path and the nanofilm thickness). Figure 4*a* shows that when the Knudsen number, *Kn*_{y}, is equal to 0.01, the temperature distribution evolutions can agree with the predictions by the diffusion model based on Fourier's law all the time. As the nanofilm thickness is comparable with the phonon mean free path, i.e. *Kn*_{y}=1.0 (figure 4*b*) and 2.0 (figure 4*c*), the temperature distributions calculated by the heat diffusion model with the effective thermal conductivity are significantly lower than the simulation results at the initial stage, and the deviations decrease and vanish with time.

In order to further analyse the deviations between the model and the MC simulations, we can define the temperature deviation function as
*T*_{MC}(*x*) and *T*_{Model}(*x*) are the temperature distributions calculated by the MC simulation and the heat diffusion model with the effective thermal conductivity, respectively. As shown in figure 5, the deviation function, *δ*, is always positive, indicating that the model always under-predicts the simulation results in the continuous heating case. As *Kn*_{y}=0.01, the deviation is only about 0.05, and the deviation increases with the Knudsen number increasing at the same moment. Furthermore, the deviation function also decreases with time. In the MC simulations, phonons emit within the nanofilms uniformly and can scatter with the lateral boundaries from the beginning. Therefore, the temperature rise can be significantly larger when compared with the heat diffusion equation based on Fourier's law even at the initial heating stage [4,5]. However, the heat diffusion model with the effective thermal conductivity cannot reflect the boundary effect at the initial stage. This point will be further investigated in §3(d) through Fourier transform.

### (c) Transient thermal transport with one-pulse internal heating

Figure 6 shows the peak temperature rise varying with time in the one-pulse heating case with Δ*t*=40 *τ*. As *Kn*_{y}=0.01, the diffusion model based on Fourier's law predicts the results calculated by the MC simulations. With the Knudsen number increasing, the phonon-boundary scattering causes the larger peak temperature at the initial heating stage and the lower cooling rate of the peak temperature rise. When compared with the simulations, the heat diffusion model with the effective thermal conductivity under-predicts both the peak temperature rise at the initial stage and the cooling rate of the temperature rise. The under-prediction of peak temperature rise comes from the failure of the modification of the effective thermal conductivity to reflect the boundary effect at the initial stage. Furthermore, since the effective thermal conductivity, *k*_{eff_IHS}, was derived from the average temperature rise in the steady state [46], it does not predict the thermal response rate. Therefore, different from the continuous heating case, the deviation between the model and the simulation does not decrease and vanish with time. When the peak temperature rise occurs at the initial heating stage, the peak temperature rise predicted by the model does not even increase with the increasing Knudsen number, indicating that the modification of the effective thermal conductivity for the heat diffusion model fails to take the boundary effect into account at the initial stage. Moreover, due to the under-prediction of the cooling rate, the model will over-predict the temperature rise with time.

The temperature distributions within the nanofilms of different Knudsen numbers for the one-pulse heating case are illustrated in figure 7. As shown in figure 7*a*, when the Knudsen number is equal to 0.01, the diffusion model based on Fourier's law predicts the temperature distribution evolutions calculated by the MC simulations. By contrast, when *Kn*_{y} increases as shown in figures 7*b* (*Kn*_{y}=1.0) and 7*c* (*Kn*_{y}=2.0), the boundary effect becomes significant, and the diffusion model with the effective thermal conductivity cannot predict the simulation results. The model under-predicts the temperature rise at the initial stage, and it will over-predict the temperature rise with time due to the under-prediction of the cooling rate. Figure 8 shows the temperature deviations varying with time in the one-pulse heating case. As *Kn*_{y}=0.01, the deviation is only about 0.05. However, different from the continuous heating case, when the thickness is comparable with the phonon mean free path (*Kn*_{y}=1.0 and 2.0), the deviation function, *δ*, changes from positive to negative with time. Furthermore, the absolute value of the deviation decreases with the Knudsen number decreasing.

Furthermore, for the nanofilm with given thickness (*Kn*_{y}=2.0), we have changed the heating duration and investigated the influence of the heating duration on the in-plane thermal transport with the boundary effect. Figure 9 shows the peak temperature varying with time in the one-pulse heating case with different heating duration (Δ*t*=400 *τ* and 4000 *τ*). The peak temperature rise will increase with the heating duration increasing, while the cooling rate of the peak temperature rise will not change. Despite the different heating durations, the diffusion model with the effective thermal conductivity both under-predicts the temperature rise at the initial stage and the cooling rate calculated by the MC simulations. The temperature distributions for the one-pulse heating with Δ*t*=400 *τ* and Δ*t*=4000 *τ* are illustrated in figure 10*a* and *b*, respectively. Similar to the case with Δ*t*=40 *τ*, the diffusion model with the effective thermal conductivity still fails to predict the simulation temperature distribution evolutions due to the under-predictions of peak temperature rise and thermal response rate. Figure 11 compares the temperature deviations between the continuous heating and the one-plus heating cases with different heating durations. The temperature deviation for the continuous heating case can decrease and nearly vanish with time, while the temperature deviation for the one-pulse heating case will always exist and change from positive to negative with time. Besides, it is also found that when the Knudsen number is given, the absolute value of the deviation decreases with the heating duration increasing, as longer heating duration makes the thermal transport state closer to the continuous heating case.

### (d) Analyses based on Fourier transform

The heat diffusion model with the effective thermal conductivity is analysed through Fourier transform to further investigate the temperature deviations. The Fourier transform is performed to the spatial domain, ignoring the *x*-directional boundaries. The heat diffusion model with the effective thermal conductivity, equation (3.3), can be rewritten as,
*g*(*x*,*t*) ranging from 0 to 1.0 describes the spatial and temporal dependence of the internal heat source. The Fourier transforms for the functions *T*(*x*,*t*) and *g*(*x*,*t*) are expressed as,
*k*_{eff}, has been cancelled out.

Equation (3.10) indicates that the modification of the effective thermal conductivity for the heat diffusion model to reflect the boundary effect cannot work at the initial heating stage, which explains the reason why the heat diffusion model with the effective thermal conductivity significantly under-predicts the temperature rise at the initial stage. Moreover, according to equation (3.9), the cooling rate of the temperature peak is directly related to the effective thermal conductivity, and the reduction of the effective thermal conductivity causes the cooling rate reduced. However, the effective thermal conductivity, *k*_{eff_IHS}, which predicts the steady-state temperature rise (as shown in figure 2), significantly under-predicts the cooling rate calculated by the MC simulations. That is because the effective thermal conductivity, *k*_{eff_IHS}, is obtained by calculating the average temperature rise in the steady state but not the thermal response rate [46].

## 4. Conclusion

The transient in-plane thermal transport in the semiconductor nanofilms with internal heating has been investigated using the MC simulations. It is found that the phonon-boundary scattering significantly suppresses the in-plane thermal transport, leading to larger temperature rise and slower thermal response rate when compared with the heat diffusion equation based on Fourier's law.

The MC simulation results are compared with the predictions by the heat diffusion model with the effective thermal conductivity which takes the boundary effect into account. In the case of continuous internal heating, the diffusion model with the effective thermal conductivity under-predicts the temperature rise predicted by the MC simulations at the initial heating stage, while the deviation between them gradually decreases and vanishes with time. By contrast, for the case of one-pulse internal heating, the model under-predicts both the peak temperature rise and the cooling rate of the temperature rise.

Furthermore, on the basis of the Fourier transform, we have found that the modification of the effective thermal conductivity for the heat diffusion model to reflect the boundary effect cannot work at the initial heating stage. Therefore, the heat diffusion model with the appropriate effective thermal conductivity can characterize the in-plane thermal transport with internal heating in steady state and the state close to steadiness for the continuous heating case, but fails to predict the temperature rise at the initial heating stage and in the one-pulse heating case. A size- and frequency-dependent effective thermal conductivity [47,48] may be a proper choice to describe the transient in-plane thermal transport with boundary effects far from the steady state, but the frequency-dependent behaviour of effective thermal conductivity is still waiting for more in-depth investigation.

## Authors' contributions

All authors contributed equally.

## Competing interests

All authors have no competing interests.

## Funding

This work is financially supported by National Natural Science Foundation of China (nos. 51322603, 51136001, 51356001), and Science Fund for Creative Research Group (no. 51321002), the Tsinghua National Laboratory for Information Science and Technology of China (TNList).

## Acknowledgements

The authors appreciate the assistance from group members in this research.

- Received November 27, 2015.
- Accepted January 18, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.