[[File(JamesAndPrabhakar_ProjectReport.pdf)]]

== Role of Normal Processes in Thermal Conductivity of Silicon ==

James Loy, Prabhakar Marepalli

ME 503 Final Project Report

'''Abstract'''
''In the past decade, with the miniaturization of electronic circuits and growing interest in microscale heat transfer, Boltzmann transport equation (BTE) emerged as an important model for predicting the thermal behavior at small scales. Recent advancements in numerical methods and computational power enabled to solve this equation rigorously by relaxing several assumptions. But the scattering term of BTE is incredibly complex, still requiring some simplifications to make the solution possible while preserving relevant physics. One approximation is the single mode relaxation time (SMRT) approximation. This approximation assumes that each phonon scatters only to a lattice equilibrium energy. For many materials near equilibrium and at high enough temperatures, this works well. However, for low temperatures, and for low dimensional materials like graphene, this model fails to include phonon scattering due to normal processes, which play an indirect role in thermal resistance. In this project, we include normal phonon scattering in the BTE through the use of a shifted equilibrium distribution proposed by Callaway. We solve the non-gray BTE using finite volume method (FVM) with coupled ordinates method (COMET) to compute the thermal conductivity and temperature distribution of silicon at different temperatures and length scales. We implement full phonon dispersion for all polarizations under isotropic assumption. The effect of including normal processes on the thermal conductivity predictions is rigorously analyzed. Our results show that ignoring the normal process overpredicts the thermal conductivity – which is a physically intuitive result. We also observe that the thermal conductivity increases and trends towards an asymptotic value as the length scale increases. By analyzing the temperature distribution, we also show that inclusion of normal processes diffuses the energy across the material – which is an expected result.''

== 1. Introduction==

With increasing miniaturization of integrated electronic circuits (ICs) following the Moore’s law [1] several challenges pop up in trying to keep up with the trend. One of the major bottlenecks is the high density localized heat generation in ICs that impair the device performance. In order to better understand the thermal behavior at those scales it is essential to create robust models that can accurately predict the device failure. As the device size in ICs these days are less than 32nm where the ballistic behavior dominates, Fourier law cannot be used to make accurate predictions. One of the widely used alternatives to use the Boltzmann Transport Equation (BTE) which models the phonon distribution – the major heat carrier in semiconductors – for a given macroscopic device conditions. Phonon is a quantum of thermal vibration which is considered as particle in BTE.

Several researchers have developed simplified versions of BTE making intuitive assumptions to make it analytically or numerically solvable. Using these assumptions they were able to predict the thermal behavior of various materials. One of the most celebrated applications of BTE developed in early 1950s is to compute the thermal conductivity of a given material as a function of temperature, composition, and geometry, etc. [2], [3]. These analyses consider various phonon scattering mechanisms responsible for the thermal conductivity of a given material and models the effect of these scattering rates on the phonon distribution. Widely considered scattering mechanisms are isotopic scattering, boundary scattering, and three-phonon scattering (Umklapp (U) and Normal (N) scattering). Isotopic scattering is caused due to the fact that any given material, by nature, has various isotopic compositions in it. Boundary scattering is dominant at low temperatures and length scales where the phonons hit the physical boundary of the material thereby causing resistance to heat flow. Three-phonon scattering takes place when three phonons interact and results in frequency modification. This process is also called intrinsic or inelastic scattering which occurs due to the anharmonic nature of interatomic potential. These mechanisms explain why even perfectly pure crystals do not have infinite thermal conductivity (In most of the modeling procedures the interatomic potential is assumed harmonic which fail to show the three-phonon scattering).

Over the last decade, with the improvement in computational power and numerical methods, many assumptions are relaxed and the BTE is solved more rigorously using entire phonon dispersion and relaxation time approximation [4]. Common methods to obtain relaxation rates include fitting the rate expressions to experimental values or perturbation theory [2], [5]. While the influence of all scattering terms on overall thermal conductivity at different temperatures and geometries has been well analyzed, the three-phonon N-processes have been neglected in most of the computations. The reason for neglecting normal processes is the premise that they conserve phonon momentum and hence do not offer thermal resistance. While the above fact is partially true, N-processes populate phonons in that region that can participate in U-processes thereby indirectly contributing to the overall thermal conductivity. Neglecting N-processes still provided reasonable comparison with experiment for materials like Si and Ge because these are 3D materials and the number of phonons that participate in U-processes are comparably higher than N-process phonons. But in case of lower dimensional (2D) materials like graphene, the population of phonons in large-wavevector region is very small and neglecting N-process provides false prediction of diverging thermal conductivity. This is because of the fact that the scarce phonon population makes them travel ballistically without many collisions. On the other hand, it was shown in [6] that on including N-processes the conductivity asymptotes to a constant value.

In this paper, we simulate the thermal conductivity of Silicon by including N-processes with full BTE model using relaxation time approximation and full phonon dispersion. N-process scattering formulation developed by Callaway [2] is used by strictly enforcing momentum conservation for N-process phonons and energy conservation for N, and U process phonons. We only consider isotopic and three-phonon scattering mechanisms in this project. The simulation is performed on Silicon owing to easy availability of its dispersion curves and relaxation rates. Isotropic assumption is made in k-space which is reasonable for Si. In our simulation, we solve the energy form of non-gray BTE [7] simultaneously with overall energy conservation (N+U processes) to extract energy distribution and temperatures respectively in a coupled fashion. This provides quick convergence compared to sequential solution of these equations (non-gray BTE and energy conservation). Then we perform a detailed study of the effect of temperature, geometry, and mainly N-process scattering on the overall bulk thermal conductivity of Si.

We organize the rest of the paper as follows. In Section 2, we explain the physics of Normal scattering processes and discuss the situations when their effect would be significant. In Section 3, we provide a literature review of different models and assumptions used to simulate and analyze the effect of N-processes. In Section 4, we use Callaway’s thermal conductivity model to make a first order prediction of thermal conductivity of silicon. Here we briefly discuss the assumptions made in the model and their implications. In Section 5, we discuss the numerical method and solution procedure we used for simulating the thermal conductivity. Here we provide a detailed description of equations involved and the formulation used. In Section 6, we present and discuss the results of our simulations. We conclude in Section 7 by providing a direction for possible future work.

==2. Normal process scattering ==

'''Three-phonon scattering'''

A three phonon scattering process results in frequency modification of the resultant phonons. They are also called as inelastic scattering events. These processes can be described by the energy and momentum relations shown in Figure 1 [8]. As shown in the figure, these processes can be classified into Normal and Umklapp processes. A Normal process conserves energy and momentum whereas Umklapp process only conserves energy. Another illustration in Figure 2 shows why U-processes do not conserve momentum. The Brillouin zone of the given material is shown in gray. Incoming phonons of wave vector k1 and k2 combine to form a single phonon of wave vector k3. The left part of the figure shows N-process in which the resultant phonon lies inside the Brillouin zone; whereas the resultant phonon for U-process (figure on the right side) has such a high wave vector that it is knocked out of the Brillouin zone. By mapping it back into Brillouin zone using reciprocal lattice vector G, we can see that the resultant phonon of wave vector k3 is in the direction opposite to that of k1 and k2. This explains why U-processes impede phonon momentum and thereby the heat flow. On the other hand, as N-processes do not impede phonon momentum they do not impede the heat flow directly. But they contribute indirectly by redistributing overall phonon population which can further participate in U-processes.

[[Image(three_phonon_events.bmp, 600px)]]

'''Figure 1''': Three-phonon scattering events [8]

[[Image(Phonon_nu_process.png, 600px)]]

'''Figure 2''': Illustration of momentum conservation by N and U processes [Wikipedia]

'''Importance of N-process'''

As discussed earlier, given the nature of N-processes it would be worthwhile to think of situations when N-process contribution is indeed significant. First, we consider how low-frequency modes interact with high frequency-modes near the Brillouin zone boundary. Considering the selection rules (see Figure 1) of U-processes, only a mode of some minimum frequency ωi can participate in them [9]. This prohibits the interaction of low-frequency modes with that of Brillouin zone boundary. But intuitively we know that these modes should somehow contribute to thermal resistance. This can be explained by the premise that N-processes that involve these low-frequency modes generate the modes this minimum frequency ωi, which can then participate in U-processes, thereby providing thermal resistance. One of the other observations on importance of N-processes is discussed in [6] where it is shown that, with N-processes, graphene’s thermal conductivity diverges with increasing flake diameter thereby providing length dependence. But by including N-processes they showed that the conductivity asymptotes to a constant value.

'''3. Literature review of Normal processes '''

In this section, we provide a brief literature review of the N-process analysis and their findings. We begin with Callaway’s phenomenological model for lattice thermal conductivity at low temperatures [2]. In this work, Callaway uses a relaxation time approximation for the scattering term of BTE, and assumes that all momentum destroying processes (isotopic, boundary, and Umklapp scattering) tend toward an equilibrium Planck distribution, whereas N processes lead to a displaced Planck distribution. Using this approximation, the scattering terms can be written as:

$\left\{\left\{\left\left( \frac\left\{\partial N\right\}\left\{\partial t\right\} \right\right)\right\}_\left\{c\right\}\right\}=\frac\left\{N\left(\lambda \right)-N\right\}\left\{\left\{\tau \right\}_\left\{N\right\}\right\}+\frac\left\{\left\{N\right\}_\left\{0\right\}-N\right\}\left\{\left\{\tau \right\}_\left\{u\right\}\right\}$
-
(1)

where is the distribution function, is the relaxation time for all normal processes and is the relaxation time for all other momentum destroying processes, is the equilibrium Planck’s distribution and is displaced Planck’s distribution defined as

$N\left(\lambda \right)=\left\{\left\{\left\left[ \exp \left\left( \frac\left\{\hbar \omega -\lambda .\mathsf\left\{k\right\}\right\}\left\{\left\{k\right\}_\left\{B\right\}T\right\} \right\right) \right\right]\right\}^\left\{-1\right\}\right\}$
-
(2)

The term is along the direction of temperature gradient and defines the amount of energy redistributed by N-processes. He computed by enforcing momentum conservation for all N-process phonons using

$\int\left\{\left\{\left\left( \frac\left\{\partial N\right\}\left\{\partial t\right\} \right\right)_\left\{N\right\}\right\}\text\left\{k\right\}\right\}\left\{\left\{d\right\}^\left\{3\right\}\right\}k=\int\left\{\frac\left\{N\left(\lambda \right)-N\right\}\left\{\left\{\tau \right\}_\left\{N\right\}\right\}\text\left\{k\right\}\left\{\left\{d\right\}^\left\{3\right\}\right\}k=0\right\}$
-
(3)

Using the assumptions that only acoustic phonons contribute to thermal conductivity all acoustic modes can be averaged using a single group velocity, and the relaxation rates for all processes can be expressed as a function of frequency and temperature he computed a simplified expression for thermal conductivity as

$k=\frac\left\{\left\{k\right\}_\left\{B\right\}\right\}\left\{2\left\{\pi \right\}^\left\{2\right\}c\right\}\left(\left\{I\right\}_\left\{1\right\}+\beta \left\{I\right\}_\left\{2\right\}\right)$

$\left\{\left\{I\right\}_\left\{1\right\}\right\}=\int_\left\{0\right\}^\left\{\left\{\left\{k\right\}_\left\{B\right\}\right\}\Theta /\hbar \right\}\left\{\left\{\left\{\tau \right\}_\left\{c\right\}\right\}\frac\left\{\left\{\left\{\hbar \right\}^\left\{2\right\}\right\}\left\{\left\{\omega \right\}^\left\{2\right\}\right\}\right\}\left\{\left\{\left\{k\right\}_\left\{B\right\}\right\}^\left\{2\right\}\left\{\left\{T\right\}^\left\{2\right\}\right\}\right\}\frac\left\{\left\{\left\{e\right\}^\left\{\hbar \omega /\left\{\left\{k\right\}_\left\{B\right\}\right\}T\right\}\right\}\right\}\left\{\left\{\left\{\left\left( \left\{\left\{e\right\}^\left\{\hbar \omega /\left\{\left\{k\right\}_\left\{B\right\}\right\}T\right\}\right\}-1 \right\right)\right\}^\left\{2\right\}\right\}\right\}\left\{\left\{\omega \right\}^\left\{2\right\}\right\}d\omega \right\}$

$\left\{\left\{I\right\}_\left\{2\right\}\right\}=\int_\left\{0\right\}^\left\{\left\{\left\{k\right\}_\left\{B\right\}\right\}\Theta /\hbar \right\}\left\{\frac\left\{\left\{\left\{\tau \right\}_\left\{c\right\}\right\}\right\}\left\{\left\{\left\{\tau \right\}_\left\{N\right\}\right\}\right\}\frac\left\{\left\{\left\{\hbar \right\}^\left\{2\right\}\right\}\left\{\left\{\omega \right\}^\left\{2\right\}\right\}\right\}\left\{\left\{\left\{k\right\}_\left\{B\right\}\right\}^\left\{2\right\}\left\{\left\{T\right\}^\left\{2\right\}\right\}\right\}\frac\left\{\left\{\left\{e\right\}^\left\{\hbar \omega /\left\{\left\{k\right\}_\left\{B\right\}\right\}T\right\}\right\}\right\}\left\{\left\{\left\{\left\left( \left\{\left\{e\right\}^\left\{\hbar \omega /\left\{\left\{k\right\}_\left\{B\right\}\right\}T\right\}\right\}-1 \right\right)\right\}^\left\{2\right\}\right\}\right\}\left\{\left\{\omega \right\}^\left\{2\right\}\right\}d\omega \right\}$

where is the combined relaxation time given by , is the frequency, is the temperature, and is the normalized Plancks constant.

Using the above expression, he computed the thermal conductivity of germanium and compared it to various experimental observations to obtain fitting constants for relaxation rates. Once the constants are obtained, he found a striking match between theory and experiments. These constants are then used to examine the effects of individual processes. The second term in the thermal conductivity expression is commonly referred as correction term that accounts for correction due to N-processes. He computed that the value of this term is 10% of the overall conductivity in case of a pure germanium sample. This model for thermal conductivity has been rigorously used to compute thermal conductivities for various materials and found to provide good match at low temperatures [3], [9].

One of the papers that made improvements to Callaway’s model is Armstrong’s paper using two-fluid model [9]. In this work Armstrong divides phonons into two groups – propagating and reservoir modes. He assumes that both low and high-frequency phonons participate in N-processes. High-frequency phonons account for the fact that they can also participate in N-processes by splitting into phonons of lower frequencies. Instead of Callaway’s displaced Planck’s distribution, he used two displaced distribution functions for N-process phonons. This is because the high frequency phonons participating in N-process cannot equilibriate to very low value of distribution function that low-frequency phonons equilibriate to. He also considers the effect of different polarizations. Using this assumptions he solves the modified BTE called Boltzmann-Peierls equation for phonons written as

$-\frac\left\{\left\{\left\{N\right\}_\left\{q\right\}\right\}-N\left(\beta \right)\right\}\left\{\left\{\left\{\tau \right\}_\left\{NN\right\}\right\}\right\}-\frac\left\{\left\{\left\{N\right\}_\left\{q\right\}\right\}-\left\{\left\{N\right\}_\left\{0\right\}\right\}\right\}\left\{\left\{\left\{\tau \right\}_\left\{R\right\}\right\}\right\}=\left\{\left\{c\right\}_\left\{q\right\}\right\}.\nabla T\frac\left\{d\left\{\left\{N\right\}_\left\{q\right\}\right\}\right\}\left\{dT\right\}$

where $\left\{\left\{N\right\}_\left\{0\right\}\right\}$ is the equilibrium distribution, $\left\{\left\{N\right\}_\left\{q\right\}\right\}$ is the phonon occupation number, $N\left(\beta \right)$ is the displaced distribution, $\left\{\left\{c\right\}_\left\{q\right\}\right\}$ is the group velocity of a given polarization, $\left\{\left\{\tau \right\}_\left\{NN\right\}\right\}$ and $\left\{\left\{\tau \right\}_\left\{R\right\}\right\}$ are the characteristic relaxation times,$\nabla T$ is the temperature gradient across the material.