Efficient simulation of multidimensional phonon transport using energy-based variance-reduced Monte Carlo formulations

We present a Monte Carlo method for obtaining solutions of the Boltzmann equation to describe phonon transport in micro- and nanoscale devices. The proposed method can resolve arbitrarily small signals (e.g., temperature differences) at small constant cost and thus represents a considerable improvement compared to traditional Monte Carlo methods, whose cost increases quadratically with decreasing signal. This is achieved via a control-variate variance-reduction formulation in which the stochastic particle description solves only for the deviation from a nearby equilibrium, while the latter is described analytically. We also show that simulation of an energy-based Boltzmann equation results in an algorithm that lends itself naturally to exact energy conservation, thereby considerably improving the simulation fidelity. Simulations using the proposed method are used to investigate the effect of porosity on the effective thermal conductivity of silicon. We also present simulations of a recently developed thermal conductivity spectroscopy process. The latter simulations demonstrate how the computational gains introduced by the proposed method enable the simulation of otherwise intractable multiscale phenomena.

Figure 1

In standard particle methods, the moments of the distribution are stochastically integrated.

Figure 2

In a control-variate formulation, the stochastic part is reduced to the calculation of the deviation from a known state, which is much smaller.

Figure 3

Example of a periodic nanostructure. Each periodic cell comprises two rectangular voids with diffusely reflecting walls. This nanostructure, and in particular the influence of the parameter $d$, is studied further in Sec. 6a.

Figure 4

Transient temperature profile in a one-dimensional ballistic system whose boundary temperatures undergo an impulsive change at $t=0$. Initially, the system is in equilibrium at temperature $T_0=300$ K. At $t=0^{+}$, the wall temperatures become $T_0\pm \Delta T$; here, $\Delta T=3$ K.

Figure 5

Heat conduction in a silicon slab due to an imposed temperature gradient in the $y$ direction. Slab is infinite in the $x$ and $y$ directions.

Figure 6

Simulation geometry. Boundaries at $z=0$ and $z=d$ are diffusely reflecting. Infinite domain in the $y$ direction is terminated by use of periodic boundaries $L=100$ nm apart.

Figure 7

Spatial variation of the axial (in the $y$ direction) heat flux in a thin film with a thickness $d=100$ nm, computed theoretically and compared to the result of the deviational simulation.

Figure 8

Theoretical values of the thin film thermal conductivity at $T_0=300$ K, computed by numerical integration of the theoretical expressions (40) and (41), and compared with the values obtained from the deviational simulation.

Figure 9

Comparison of relative statistical uncertainties for equilibrium systems at temperature $T_1$ with $\Delta T=T_1-T_0$ and $T_0=300$ K.

Figure 10

(a) Temperature field in a unit cell of a periodic nanoporous material. (b) Thermal conductivity as a function of parameter $d$.

Figure 11

(a) Thermal conductivity accumulation as a function of the mean free path. The bulk conductivity was computed by numerical integration of the thermal conductivity per unit frequency, $\tau v^2{C}_{\omega }/3$ (Refs.15,24, and1); here, $C_{\omega }$ is the heat capacity per frequency unit. (b) Normalized thermal conductivity accumulation, highlighting the influence of ballistic effects on the thermal conductivity.

Figure 12

System composed of a slab of aluminum on a semi-infinite silicon wafer, used for transient thermoreflectance experiments. At $t=0$, a laser pulse induces a temperature field $T(r,z,t)$. The temperature field evolution after the pulse is computed by assuming that the aluminum surface is adiabatic.

Figure 13

Variance-reduced temperature field in an aluminum slab and supporting silicon wafer after initial heating by a laser pulse. The picture shows the aluminum slab (100 nm thickness) and a portion of the silicon wafer (100 nm thickness).

Figure 14

Surface temperature at the hot spot (averaged over the region $0⩽r⩽2\mu \text{m}$, $0⩽z⩽5\text{nm}$) as a function of time after initial heating by laser pulse. The difference from the solution based on the Fourier model is a result of ballistic effects.

Figure 15

At a given point in space, the solid angle can be divided into three distinct regions in which the distribution of phonons is known; here, ${\theta }_l$ is given by $\cos \left({\theta }_l\right)=x/\left[V_g(\omega ,p)t\right]$, while ${\theta }_r$ is given by $\cos \left({\theta }_r\right)=(L-x)/\left[V_g(\omega ,p)t\right].$