Variance-reduced Monte Carlo simulation of phonon transport in silicon films
The variance-reduced Monte Carlo simulation of phonon transport solves the following equation:
\[ \frac{\partial e^{\text{d}}}{\partial t} + \boldsymbol{v}_{\mathrm{g}} \cdot \nabla e^{\text{d}} = -\frac{e^{\mathrm{d}} - (e^{\mathrm{eq}}(T)-e^{\mathrm{eq}}(T_{\mathrm{eq}}))}{\tau} \]
Here, relaxation time approximation is used, and \(\tau\) takes ab initio calculated value, see the previous post. In this equation, \(e^{\mathrm{d}} = \hbar\omega (f(\omega)-f_{\mathrm{BE}}(\omega, T_{\mathrm{eq}}))\) for each phonon mode, where \(f\) is the distribution function (not necessarily equilibrium). While \(e^{\mathrm{eq}}(T)=\hbar\omega (f_{\mathrm{BE}}(\omega, T)-f_{\mathrm{BE}}(\omega, T_{\mathrm{eq}}))\) follows the Bose-Einstein distribution. Note that here \(T_{\mathrm{eq}}\) is a reference temperature, and \(T\) is the local temperature that deviates from \(T_{\mathrm{eq}}\).
Method
To solve the above equation in a Monte Carlo way, we need to understand the meaning of each sampling step. The traditional Monte Carlo is easy to understand, which simulates the original phonon BTE:
\[ \frac{\partial f}{\partial t} + \boldsymbol{v}_{\mathrm{g}} \cdot \nabla f = -\frac{f - f_{\mathrm{eq}}(T)}{\tau} \]
and it is direct that each sample represents a fixed amount of phonons. Analogously, in the energy-based variance-reduced formulation, each sample represents a fixed amount of energy, and please keep that in mind, which is important for understanding the rules of initialization. Also, it is also worth noting that \(e^{\mathrm{d}}\) can be positive or negative, and that will give each sampled phonon a sign of energy.
To put things simple enough, here we just consider the cross-plane transport of a silicon film. In that case, only temperature boundaries are considered.
Phonon initialization
To draw a phonon for the start, we have to make sure that the boundary and initial conditions are met. For the temperature boundaries, one thing that we are sure is that, from \(t=0\) to \(t_{\mathrm{total}}\), the total energy input from the left boundary is:
\[ E_{\mathrm{left}} = t_{\text{total}}\sum_{\boldsymbol{q}s, \text{ s.t. } v_{gx}>0} \hbar\omega f(\omega, T_{\text{left}}) v_{gx} \]
thus the deviated energy is:
\[ E_{\mathrm{left}}^{\mathrm{d}} = t_{\text{total}}\sum_{\boldsymbol{q}s, \text{ s.t. } v_{gx}>0} \underbrace{\hbar\omega (f(\omega, T_{\text{left}})-f(\omega, T_{\text{eq}}))}_{\approx c_{\boldsymbol{q}s} \Delta T_{\text{left}}} v_{gx} \]
And the right boundary is the same:
\[ E_{\mathrm{right}}^{\mathrm{d}} = t_{\text{total}}\sum_{\boldsymbol{q}s, \text{ s.t. } v_{gx}<0} \underbrace{\hbar\omega (f(\omega, T_{\text{right}})-f(\omega, T_{\text{eq}}))}_{\approx c_{\boldsymbol{q}s} \Delta T_{\text{right}}} v_{gx} \]
For the internal area, the story is not the same, as the deviated energy does nothing with the time \(t_{\text{total}}\), and it is only related to the initial temperature:
\[ E_{\mathrm{internal}}^{\mathrm{d}} = \sum_{\boldsymbol{q}s} \hbar\omega (f(\omega, T_{\text{init}})-f(\omega, T_{\text{eq}})) L \approx \sum_{\boldsymbol{q}s} c_{\boldsymbol{q}s} \Delta T_{\text{init}} L \]
Recall that each sample represents a fixed amount of energy, we should thus draw each phonon mode according to their contribution to the total deviated energy. It makes the things worse as the deviated energy can be positive or negative, but the good thing is that we can add an abosolute value. That leads to the following strategy:
- If we want to study \(N\) phonons in the simulation, we should calculate the energy that each phonon we draw should have, which is \(\frac{|E_{\mathrm{left}}^{\mathrm{d}}| + |E_{\mathrm{internal}}^{\mathrm{d}}| + |E_{\mathrm{right}}^{\mathrm{d}}|}{N}\). This is the abosolute value, and the sign will be determined later.
- Draw a random number \(r\) from
\([0, 1]\):
- If \(r <
\frac{|E_{\mathrm{left}}^{\mathrm{d}}|}{|E_{\mathrm{left}}^{\mathrm{d}}|
+ |E_{\mathrm{internal}}^{\mathrm{d}}| +
|E_{\mathrm{right}}^{\mathrm{d}}|}\), the phonon is drawn from
the left boundary:
- Draw the phonon mode according to its contribution to \(E_{\mathrm{left}}^{\mathrm{d}}\), thus the probability is proportional to \(c_{\boldsymbol{q}s}v_{\boldsymbol{q}s}\), but make sure that \(v_{gx}>0\).
- The sign of the energy is determined by the sign of \(E_{\mathrm{left}}^{\mathrm{d}}\).
- If \(\frac{|E_{\mathrm{left}}^{\mathrm{d}}|}{|E_{\mathrm{left}}^{\mathrm{d}}|
+ |E_{\mathrm{internal}}^{\mathrm{d}}| +
|E_{\mathrm{right}}^{\mathrm{d}}|} \leq r <
\frac{|E_{\mathrm{left}}^{\mathrm{d}}| +
|E_{\mathrm{internal}}^{\mathrm{d}}|}{|E_{\mathrm{left}}^{\mathrm{d}}| +
|E_{\mathrm{internal}}^{\mathrm{d}}| +
|E_{\mathrm{right}}^{\mathrm{d}}|}\), the phonon is drawn from
the internal area:
- Draw the phonon mode according to its contribution to \(E_{\mathrm{internal}}^{\mathrm{d}}\), thus the probability is proportional to \(c_{\boldsymbol{q}s}\).
- Draw another number to determine the initial position of the phonon.
- The sign of the energy is determined by the sign of \(E_{\mathrm{internal}}^{\mathrm{d}}\).
- If \(r \geq
\frac{|E_{\mathrm{left}}^{\mathrm{d}}| +
|E_{\mathrm{internal}}^{\mathrm{d}}|}{|E_{\mathrm{left}}^{\mathrm{d}}| +
|E_{\mathrm{internal}}^{\mathrm{d}}| +
|E_{\mathrm{right}}^{\mathrm{d}}|}\), the phonon is drawn from
the right boundary:
- Draw the phonon mode according to \(E_{\mathrm{right}}^{\mathrm{d}}\), thus the probability \(\propto c_{\boldsymbol{q}s}v_{\boldsymbol{q}s}\), but make sure that \(v_{gx}<0\).
- If \(r <
\frac{|E_{\mathrm{left}}^{\mathrm{d}}|}{|E_{\mathrm{left}}^{\mathrm{d}}|
+ |E_{\mathrm{internal}}^{\mathrm{d}}| +
|E_{\mathrm{right}}^{\mathrm{d}}|}\), the phonon is drawn from
the left boundary:
- Determine the initial time of the phonon, which is drawn uniformly from \([0, t_{\mathrm{total}}]\).
Phonon advection and scattering
After initialized a phonon, we need to advect it and scatter it. We first have to determine when the phonon will be scattered. Due to the scattering time \(\tau_{\boldsymbol{q}s}\), we draw a random scattering time \(t_{\mathrm{scat}}\) following the Poisson process. Draw a random number \(r\) from \([0, 1]\), and the scattering time is:
\[ t_{\mathrm{scat}} = -\tau_{\boldsymbol{q}s} \log r \]
Before the scattering, the phonon is advected according to its group velocity. Once the phonon goes beyond the boundary, it is removed from the simulation. If the phonon is scattered, we need to draw a new phonon mode, but keeps the same amount and sign of energy. The new phonon mode is drawn according to the scattering probability:
\[ P_{\boldsymbol{q}s} \propto \frac{c_{\boldsymbol{q}s}}{\tau_{\boldsymbol{q}s}} \]
The key is to understand why we draw phonon from \(c_{\boldsymbol{q}s}\) inside the area and from \(c_{\boldsymbol{q}s}v_{\boldsymbol{q}s}\) at the boundaries for initialization, but draw phonon from \(c_{\boldsymbol{q}s}/\tau_{\boldsymbol{q}s}\) for scattering. For the initialization, always keep in mind that each sample represents a fixed amount of energy, thus we need to draw phonon according to their contribution to the total deviated energy. For the scattering, the key is to understand the source term in the BTE, which is \(\left(-\frac{e^{\mathrm{d}} - (e^{\mathrm{eq}}(T)-e^{\mathrm{eq}}(T_{\mathrm{eq}}))}{\tau_{\boldsymbol{q}s}}\right)\), and the first term controls the annihilation, the second term thus controls the creation, which is proportional to \(c_{\boldsymbol{q}s}/\tau_{\boldsymbol{q}s}\).
It is also important to understand why each phonon is deleted once it goes beyond the boundary. This is because, (take the left boundary as an example) we have already considered all energy goes from left to right, and if the phonon bounces back or diffuses back, it will be counted twice.
Results
\(L=200\,\mathrm{nm}\): strong size effect
The "Fourier" lines are calculated by solving the thermal diffusion
equation using fipy. The diffusivity is calculated by the
ab initio data.
In this case, \(k_{\text{eff}} = \frac{|q|}{|T_{\text{left}} - T_{\text{right}}|/L} = 48.268\,\mathrm{W/mK}\) (34.56% of bulk). And it can be seen that Fourier's law propagates heat way too fast.
\(L=10\,\mathrm{\mu m}\): almost Fourier
Here \(k_{\text{eff}} = 119.370\,\mathrm{W/mK}\) (85.58% of bulk). Thermal diffusion predicted by Fourier's law is still faster, but not that much.
\(L=100\,\mathrm{\mu m}\): Fourier
Here \(k_{\text{eff}} = k_{\text{bulk}}\). Fourier's law is fully recovered.
Effective thermal conductivity vs. film thickness
Here I compare the effective thermal conductivity calculated by the present VRMC and the results reported by Jeong et al. using MD.
The results are remarkably close. But notice that the predicted bulk thermal conductivities are different by 10%. Anyway, it can be at least concluded that the trend is the same.
References
- Péraud and Hadjiconstantinou, PRB 2011, 10.1103/PhysRevB.84.205331.
- Péraud and Hadjiconstantinou, APL 2012, 10.1063/1.4757607.
- Péraud et al., Annual Review of Heat Transfer 2014, 10.1615/AnnualRevHeatTransfer.2014007381.
- 1D_KMC_matlab: https://github.com/jeanphilippeperaud/1D_KMC_matlab.
- Jeong et al., JAP 2012, 10.1063/1.4710993.