== Role of Normal Processes in Thermal Conductivity of Silicon ==
James Loy, Prabhakar Marepalli
ME 503 Final Project Report
''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  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. , . 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 . Common methods to obtain relaxation rates include fitting the rate expressions to experimental values or perturbation theory , . 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  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  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  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 ==
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 . 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.
'''Figure 1''': Three-phonon scattering events 
'''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 . 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  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 . 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:
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
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
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
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 , .
One of the papers that made improvements to Callaway’s model is Armstrong’s paper using two-fluid model . 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
where is the equilibrium distribution, is the phonon occupation number, is the displaced distribution, is the group velocity of a given polarization, and are the characteristic relaxation times, is the temperature gradient across the material.
Recently, several rigorous computational simulations have been performed by relaxing many assumptions that were made to simplify the equation. For example, Yunfei Chen et al.,  used Monte Carlo (MC) simulation to solve the BTE to compute the thermal conductivity of silicon nanowire. In this work, N-processes are included using genetic algorithm to generate the phonons that satisfy momentum conservation and tend towards displaced equilibrium distribution defined by Callaway.
'''4. Thermal conductivity of Silicon using Callaway’s analysis'''
As a first approximation to our simulation, we start with applying Callaway’s model to compute the thermal conductivity of silicon. Both the terms of thermal conductivity expression are retained. We use the same relaxation rate expressions used by Callaway given as
where represents isotopic scattering;
includes the Umklapp processes with containing the exponential temperature factor where is the Debye temperature; and
represents the boundary scattering with
being the correction factor due to both smoothness of the surface and the finite length to thickness ratio of the sample . For N-processes we use
. The combined relaxation rate is then defined as
The above relaxation rates are used in the thermal conductivity expression to fit the constants , , and . Though ''B1'' implicitly has temperature dependence, we neglect that dependence while fitting. Using the experimental results shown in Table 1 ,  we obtain the following values of fitting constants:
''A'' = .22e-44 sec^3^
''B1+B2'' = 2.9e-24 sec/deg^3^
''F'' = .8
||T(K)||k (W/cmK)|||T (K)||K (W/cmK)||
||2|| .44|| 200|| 2.66||
||4|| 3.11|| 300|| 1.56||
||6|| 8.99|| 400|| 1.05||
||8|| 16.4|| 500|| 0.80||
||10|| 24.0|| 600|| 0.64||
||20|| 47.7|| 700|| 0.52||
||30|| 44.2|| 800|| 0.43||
||40|| 36.6|| 900|| 0.36||
||50|| 28.0|| 1000|| 0.31||
||100|| 9.13|| 1100|| 0.28||
'''Table1''': Thermal conductivity measurements of Silicon , .
We compare the thermal conductivity obtained by substituting fitted constants into equation (1.4) and observe a good comparison with experiments at low temperatures, see Figure 3. The fitted equation is really helpful to quickly investigate the effects of individual scattering events. For example, we can fix all the constants and just vary the constant to examine the change of thermal conductivity with various isotope scattering events. The value of A for a given isotopic concentration can be found from Klemens expression  for isotopic scattering. This is how several researchers exploit this equation to better understand the thermal conductivity behavior of different materials.
'''Figure 3''': Callaway’s model for thermal conductivity. The constants of relaxation rates are fitted to match the model with experiment at low temperatures.
While the above means of using Callaway’s expression would be extremely helpful, our focus in this project is to better understand the model itself. A small discussion on correction term ( ) and the justification for neglecting it might shed some light in this direction. As we discussed earlier, Callaway used a combined relaxation rate
in his derivation. But as N-processes do not actually contribute to thermal resistance, the term under predicts the thermal conductivity. Most of the times the under prediction is very low and hence can be safely neglected. But in some cases, it might the prediction might differ by around 20%. An example is when isotopic scattering is absent. In this case, the predicted thermal conductivity with and without the correction term is more than 20% of its overall value. This can be qualitatively seen in Figure 4. (Note that the y-axis here is logarithmic). So it can be explained that the correction term compensates for the fact that using provides more resistance to the conductivity.
'''Figure 4''': Meaning of correction term in Callaway’s model for thermal conductivity. The correction term compensates the fact that the overall relaxation rate used in Callaway’s model underpredicts the thermal conductivity
'''5. Numerical Method'''
All of the governing equations are solved using the Finite Volume method . The BTE is a partial differential equation involving two vector spaces (physical space and wave vector space). We therefore must discretize both spaces accordingly. With this, we will arrive at a linear system with at least one equation per discretized physical space and per discretized wave vector space.
The physical domain is discretized into several arbitrary convex polyhedral. Here we have chosen a square domain which can easily be discretized using a non-uniform structured grid, shown below:
'''Figure 5''': 80x80 grid used as the physical domain. The side length is unity and is scaled to the appropriate domain size.
To discretize the BTE, we integrate over a finite control volume in physical space and wave vector space then apply the divergence theorem to the convective operator as shown below:
We now apply our discretization. In the wave vector space, a central difference approximation is made on both sides of the equation, thus the finite volume of k-space appears on both sides and is dropped. We arrive at the following discrete equation set:
The energy conservation equation can be rearranged as follows:
When we apply a second order finite volume discretization we arrive at the following:
Because of the tight inter-equation coupling, caused by the in-scattering term of the BTE, it is advantageous to visit each physical cell and solve all points in wave vector space and the energy conservation equation in a coupled fashion. This procedure is very similar to that shown in , thus we have adopted the same name, the Coupled Ordinates METhod (COMET). In COMET, the success hinges on the ability to solve all the BTE’s and the energy conservation equation. To do this, we must determine the common variable which couples all said equations. The lattice temperature shows up in the function for the equilibrium distribution function, albeit nonlinear. To extract the lattice temperature equation, we simply linearize the equilibrium function using a Taylor series expansion:
The derivative used is a very familiar value, the specific heat of the specific phonon frequency at the previous iterations temperature. For this procedure, we will construct a matrix which solves for the correction to the previous values of and ''T''. In this way, the residual of the previous iteration act as a source for the correction equation. With this, the BTE and the energy conservation equations become:
where is with ''old'' and ''new'' representing the current and previous iteration values of simulation. ''R,,BTE,,'' is the residual for the BTE while ''R,,energy,,'' is the residual for the energy equation. These equations form a linear system of order , where is the number if points in wave vector space. The shape of the matrix is especially convenient, forming an arrowhead pointing to the lower right (all diagonals populated, last row populated, last column populated) which can be solved directly in ''O(N)'' operations.
To keep the advantageous shape of the arrowhead matrix, we do not include the momentum conservation equation in the coupled solve. The shifted equilibrium value used is the prevailing value. After the BTE and the energy conservation equation are solved, we use the new value of the lattice temperature and calculate the lambda vector which we then use to update the shifted equilibrium energy.
The solution procedure for COMET begins with the initialization of all values. The first cell is visited whereby we solve the coupled BTE and energy conservation equation. We use the new lattice temperature to solve the momentum conservation equation, giving us a new lambda vector which we use to update the shifted equilibrium energy. Each cell is visited in turn. The residuals are collected for each cell, and the heat balance is assessed. If the residual has reached the prescribed tolerance and heat balance is sufficiently satisfied, the procedure exits, otherwise the procedure begins again. A visual representation of the solution procedure is shown below.
'''Figure 6''': Flow chart for the COMET solution procedure.
We will use an isotropic Brillouin zone, using the dispersion of silicon in the  direction. The dispersion is taken using the environment dependent interatomic potential . A plot of the dispersion is shown in Figure 7. To discretize the Brillouin zone, we will use a spherical coordinate system and divide the Brillouin zone sphere into control volumes. is the number of discretizations in the polar angle, is the number of discretizations in the azimuthal angle, and is the number of discretizations in the wave vector magnitude. See Figure 8 for a visual representation.
'''Figure 7''': Dispersion relation for silicon in the  direction .
The Umklapp scattering rates take the following form: . Here the constants ''B'' and ''C'' are 1.73x10-19 ''s/K'' and 137.39 K [Mingo et al.]. The Normal scattering rates take the following form: . Here the constant ''B,,l,,'' is 2x10-24 ''s/K^3^''.
The domain is shown in Figure 9. We will have temperature boundary conditions on all sides of the domain. Three of the boundaries will have the same temperature, while the fourth boundary will be 1 Kelvin higher. A schematic of the domain is shown below.
'''Figure 8''': Schematic of k-space discretization
'''Figure 9''': Sketch of the simulation domain. .
To extract the thermal conductivity, we will use the analytical solution for the heat rate through the bottom wall, shown as follows: