[[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 nongray 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 threephonon 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. Threephonon 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 threephonon 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 threephonon Nprocesses 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, Nprocesses populate phonons in that region that can participate in Uprocesses thereby indirectly contributing to the overall thermal conductivity. Neglecting Nprocesses 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 Uprocesses are comparably higher than Nprocess phonons. But in case of lower dimensional (2D) materials like graphene, the population of phonons in largewavevector region is very small and neglecting Nprocess 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 Nprocesses the conductivity asymptotes to a constant value.
In this paper, we simulate the thermal conductivity of Silicon by including Nprocesses with full BTE model using relaxation time approximation and full phonon dispersion. Nprocess scattering formulation developed by Callaway [2] is used by strictly enforcing momentum conservation for Nprocess phonons and energy conservation for N, and U process phonons. We only consider isotopic and threephonon 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 kspace which is reasonable for Si. In our simulation, we solve the energy form of nongray 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 (nongray BTE and energy conservation). Then we perform a detailed study of the effect of temperature, geometry, and mainly Nprocess 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 Nprocesses. 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 ==
'''Threephonon 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 Uprocesses 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 Nprocess in which the resultant phonon lies inside the Brillouin zone; whereas the resultant phonon for Uprocess (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 Uprocesses impede phonon momentum and thereby the heat flow. On the other hand, as Nprocesses 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 Uprocesses.
[[Image(three_phonon_events.bmp, 600px)]]
'''Figure 1''': Threephonon scattering events [8]
[[Image(Phonon_nu_process.png, 600px)]]
'''Figure 2''': Illustration of momentum conservation by N and U processes [Wikipedia]
'''Importance of Nprocess'''
As discussed earlier, given the nature of Nprocesses it would be worthwhile to think of situations when Nprocess contribution is indeed significant. First, we consider how lowfrequency modes interact with high frequencymodes near the Brillouin zone boundary. Considering the selection rules (see Figure 1) of Uprocesses, only a mode of some minimum frequency ωi can participate in them [9]. This prohibits the interaction of lowfrequency 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 Nprocesses that involve these lowfrequency modes generate the modes this minimum frequency ωi, which can then participate in Uprocesses, thereby providing thermal resistance. One of the other observations on importance of Nprocesses is discussed in [6] where it is shown that, with Nprocesses, graphene’s thermal conductivity diverges with increasing flake diameter thereby providing length dependence. But by including Nprocesses 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 Nprocess 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:
$\{\{\backslash left(\; \backslash frac\{\backslash partial\; N\}\{\backslash partial\; t\}\; \backslash right)\}\_\{c\}\}=\backslash frac\{N(\backslash lambda\; )N\}\{\{\backslash tau\; \}\_\{N\}\}+\backslash frac\{\{N\}\_\{0\}N\}\{\{\backslash tau\; \}\_\{u\}\}$

(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(\backslash lambda\; )=\{\{\backslash left[\; \backslash exp\; \backslash left(\; \backslash frac\{\backslash hbar\; \backslash omega\; \backslash lambda\; .\backslash mathsf\{k\}\}\{\{k\}\_\{B\}T\}\; \backslash right)\; \backslash right]\}^\{1\}\}$

(2)
The term is along the direction of temperature gradient and defines the amount of energy redistributed by Nprocesses. He computed by enforcing momentum conservation for all Nprocess phonons using
$\backslash int\{\{\backslash left(\; \backslash frac\{\backslash partial\; N\}\{\backslash partial\; t\}\; \backslash right)\_\{N\}\}\backslash text\{k\}\}\{\{d\}^\{3\}\}k=\backslash int\{\backslash frac\{N(\backslash lambda\; )N\}\{\{\backslash tau\; \}\_\{N\}\}\backslash text\{k\}\{\{d\}^\{3\}\}k=0\}$

(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=\backslash frac\{\{k\}\_\{B\}\}\{2\{\backslash pi\; \}^\{2\}c\}(\{I\}\_\{1\}+\backslash beta\; \{I\}\_\{2\})$
$\{\{I\}\_\{1\}\}=\backslash int\_\{0\}^\{\{\{k\}\_\{B\}\}\backslash Theta\; /\backslash hbar\; \}\{\{\{\backslash tau\; \}\_\{c\}\}\backslash frac\{\{\{\backslash hbar\; \}^\{2\}\}\{\{\backslash omega\; \}^\{2\}\}\}\{\{\{k\}\_\{B\}\}^\{2\}\{\{T\}^\{2\}\}\}\backslash frac\{\{\{e\}^\{\backslash hbar\; \backslash omega\; /\{\{k\}\_\{B\}\}T\}\}\}\{\{\{\backslash left(\; \{\{e\}^\{\backslash hbar\; \backslash omega\; /\{\{k\}\_\{B\}\}T\}\}1\; \backslash right)\}^\{2\}\}\}\{\{\backslash omega\; \}^\{2\}\}d\backslash omega\; \}$
$\{\{I\}\_\{2\}\}=\backslash int\_\{0\}^\{\{\{k\}\_\{B\}\}\backslash Theta\; /\backslash hbar\; \}\{\backslash frac\{\{\{\backslash tau\; \}\_\{c\}\}\}\{\{\{\backslash tau\; \}\_\{N\}\}\}\backslash frac\{\{\{\backslash hbar\; \}^\{2\}\}\{\{\backslash omega\; \}^\{2\}\}\}\{\{\{k\}\_\{B\}\}^\{2\}\{\{T\}^\{2\}\}\}\backslash frac\{\{\{e\}^\{\backslash hbar\; \backslash omega\; /\{\{k\}\_\{B\}\}T\}\}\}\{\{\{\backslash left(\; \{\{e\}^\{\backslash hbar\; \backslash omega\; /\{\{k\}\_\{B\}\}T\}\}1\; \backslash right)\}^\{2\}\}\}\{\{\backslash omega\; \}^\{2\}\}d\backslash omega\; \}$
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 Nprocesses. 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 twofluid model [9]. In this work Armstrong divides phonons into two groups – propagating and reservoir modes. He assumes that both low and highfrequency phonons participate in Nprocesses. Highfrequency phonons account for the fact that they can also participate in Nprocesses by splitting into phonons of lower frequencies. Instead of Callaway’s displaced Planck’s distribution, he used two displaced distribution functions for Nprocess phonons. This is because the high frequency phonons participating in Nprocess cannot equilibriate to very low value of distribution function that lowfrequency phonons equilibriate to. He also considers the effect of different polarizations. Using this assumptions he solves the modified BTE called BoltzmannPeierls equation for phonons written as
$\backslash frac\{\{\{N\}\_\{q\}\}N(\backslash beta\; )\}\{\{\{\backslash tau\; \}\_\{NN\}\}\}\backslash frac\{\{\{N\}\_\{q\}\}\{\{N\}\_\{0\}\}\}\{\{\{\backslash tau\; \}\_\{R\}\}\}=\{\{c\}\_\{q\}\}.\backslash nabla\; T\backslash frac\{d\{\{N\}\_\{q\}\}\}\{dT\}$

(5)
where $\{\{N\}\_\{0\}\}$ is the equilibrium distribution, $\{\{N\}\_\{q\}\}$ is the phonon occupation number, $N(\backslash beta\; )$ is the displaced distribution, $\{\{c\}\_\{q\}\}$ is the group velocity of a given polarization, $\{\{\backslash tau\; \}\_\{NN\}\}$ and $\{\{\backslash tau\; \}\_\{R\}\}$ are the characteristic relaxation times,$\backslash nabla\; T$ 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., [10] used Monte Carlo (MC) simulation to solve the BTE to compute the thermal conductivity of silicon nanowire. In this work, Nprocesses 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
$\{\{\backslash tau\; \}\_\{u\}\}^\{1\}=A\{\{\backslash omega\; \}^\{4\}\}+\{\{B\}\_\{1\}\}\{\{T\}^\{3\}\}\{\{\backslash omega\; \}^\{2\}\}+c/LF$

(6)
where $A\{\{\backslash omega\; \}^\{4\}\}$ represents isotopic scattering;
$\{\{B\}\_\{1\}\}\{\{T\}^\{3\}\}\{\{\backslash omega\; \}^\{2\}\}$ includes the Umklapp processes with $\{\{B\}\_\{1\}\}$ containing the exponential temperature factor $\{\{e\}^\{\backslash Theta\; /aT\}\}$ where $\backslash Theta$ is the Debye temperature; and
$c/LF$ represents the boundary scattering with
$F$ being the correction factor due to both smoothness of the surface and the finite length to thickness ratio of the sample [3]. For Nprocesses we use
$\{\{\backslash tau\; \}\_\{N\}\}^\{1\}=\{\{B\}\_\{2\}\}\{\{T\}^\{3\}\}\{\{\backslash omega\; \}^\{2\}\}$. The combined relaxation rate is then defined as
$\{\{\backslash tau\; \}\_\{c\}\}^\{1\}=A\{\{\backslash omega\; \}^\{4\}\}+(\{\{B\}\_\{1\}\}+\{\{B\}\_\{2\}\})\{\{T\}^\{3\}\}\{\{\backslash omega\; \}^\{2\}\}+c/LF$

(7)
The above relaxation rates are used in the thermal conductivity expression to fit the constants $A$, $(\{\{B\}\_\{1\}\}+\{\{B\}\_\{2\}\})$, and $F$. Though ''B1'' implicitly has temperature dependence, we neglect that dependence while fitting. Using the experimental results shown in Table 1 [11], [12] we obtain the following values of fitting constants:
''A'' = .22e44 sec^3^
''B1+B2'' = 2.9e24 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
1504.10
'''Table1''': Thermal conductivity measurements of Silicon [11], [12].
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 [13] for isotopic scattering. This is how several researchers exploit this equation to better understand the thermal conductivity behavior of different materials.
[[Image(Callaways_model_fitting.png)]]
'''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
$\{\{\backslash tau\; \}\_\{c\}\}^\{1\}=\{\{\backslash tau\; \}\_\{N\}\}^\{1\}+\{\{\backslash tau\; \}\_\{u\}\}^\{1\}$ in his derivation. But as Nprocesses do not actually contribute to thermal resistance, the $\{\{\backslash tau\; \}\_\{c\}\}$ 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 yaxis here is logarithmic). So it can be explained that the correction term compensates for the fact that using $\{\{\backslash tau\; \}\_\{c\}\}$ provides more resistance to the conductivity.
[[Image(Callaways_model_with_and_wo_correction.png)]]
'''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 [14]. 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 nonuniform structured grid, shown below:
[[Image(80x80grid.tif)]]
'''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:
$\backslash int\backslash limits\_\{\backslash Delta\; \backslash mathsf\{K\}\}\{\backslash int\backslash limits\_\{\backslash Delta\; V\}\{\backslash nabla\; \backslash cdot\; \backslash left(\; \backslash mathsf\{v\}e\text{'}\text{'}\; \backslash right)dV\{\{d\}^\{3\}\}\backslash mathsf\{K\}\}\}=\backslash int\backslash limits\_\{\backslash Delta\; \backslash mathsf\{K\}\}\{\backslash int\backslash limits\_\{\backslash Delta\; \backslash mathsf\{A\}\}\{e\text{'}\text{'}\backslash mathsf\{v\}\backslash cdot\; d\backslash mathsf\{A\}\{\{d\}^\{3\}\}\backslash mathsf\{K\}=\}\}\backslash int\backslash limits\_\{\backslash Delta\; \backslash mathsf\{K\}\}\{\backslash int\backslash limits\_\{\backslash Delta\; V\}\{\backslash frac\{\{\{e\}^\{0\}\}e\text{'}\text{'}\}\{\{\{\backslash tau\; \}\_\{U\}\}\}+\backslash frac\{e\_\{\backslash mathsf\{\backslash lambda\; \}\}^\{0\}e\text{'}\text{'}\}\{\{\{\backslash tau\; \}\_\{N\}\}\}dV\{\{d\}^\{3\}\}\backslash mathsf\{K\}\}\}$

(8)
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 kspace appears on both sides and is dropped. We arrive at the following discrete equation set:
$\backslash sum\backslash limits\_\{f\}\{\{\{e\}\_\{f\}\}\text{'}\text{'}\}v\backslash cdot\; \backslash Delta\; \{\{\backslash mathsf\{A\}\}\_\{f\}\}+e\text{'}\text{'}\backslash Delta\; V\backslash left(\; \backslash tau\; \_\{N\}^\{1\}+\backslash tau\; \_\{U\}^\{1\}\; \backslash right)\{\{e\}^\{0\}\}\backslash frac\{\backslash Delta\; V\}\{\{\{\backslash tau\; \}\_\{U\}\}\}e\_\{\backslash mathsf\{\backslash lambda\; \}\}^\{0\}\backslash frac\{\backslash Delta\; V\}\{\{\{\backslash tau\; \}\_\{N\}\}\}=0$

(9)
The energy conservation equation can be rearranged as follows:
$\backslash int\{\backslash frac\{\{\{e\}^\{0\}\}\}\{\{\{\backslash tau\; \}\_\{U\}\}\}\{\{d\}^\{3\}\}\backslash mathsf\{K\}=\}\backslash int\{\backslash frac\{e\text{'}\text{'}\}\{\{\{\backslash tau\; \}\_\{U\}\}\}\{\{d\}^\{3\}\}\backslash mathsf\{K\}\}\backslash int\{\backslash frac\{e\_\{\backslash mathsf\{\backslash lambda\; \}\}^\{0\}e\text{'}\text{'}\}\{\{\{\backslash tau\; \}\_\{N\}\}\}\{\{d\}^\{3\}\}\backslash mathsf\{K\}\}$

(10)
When we apply a second order finite volume discretization we arrive at the following:
$\backslash sum\{\{\{e\}^\{0\}\}\backslash frac\{\backslash Delta\; \backslash mathsf\{K\}\}\{\{\{\backslash tau\; \}\_\{U\}\}\}\}=\backslash sum\{e\text{'}\text{'}\backslash left(\; \backslash tau\; \_\{N\}^\{1\}+\backslash tau\; \_\{U\}^\{1\}\; \backslash right)\backslash Delta\; \backslash mathsf\{K\}\}\backslash sum\{e\_\{\backslash mathsf\{\backslash lambda\; \}\}^\{0\}\backslash frac\{\backslash Delta\; \backslash mathsf\{K\}\}\{\{\{\backslash tau\; \}\_\{U\}\}\}\}$

(11)
Because of the tight interequation coupling, caused by the inscattering 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 [15], 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:
$\{\{e\}^\{0,new\}\}=\{\{e\}^\{0,old\}\}+\backslash left(\; \{\{T\}^\{new\}\}\{\{T\}^\{old\}\}\; \backslash right)\{\{\backslash left(\; \backslash frac\{\backslash partial\; \{\{e\}^\{0\}\}\}\{\backslash partial\; T\}\; \backslash right)\}^\{old\}\}$

(12)
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 $e\text{'}\text{'}$ 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:
$\backslash begin\{align\}$
& \sum\limits_{f}{\Delta {{e}_{f}}''}\mathsf{v}\cdot \Delta {{\mathsf{A}}_{f}}+\Delta e''\Delta V\left( \tau _{N}^{1}+\tau _{U}^{1} \right)\Delta T\frac{\Delta V}{{{\tau }_{U}}}{{\left( \frac{\partial {{e}^{0}}}{\partial T} \right)}^{old}}={{R}_{BTE}} \\
& \Upsilon \sum{\Delta e''}\left( \tau _{N}^{1}+\tau _{U}^{1} \right)\Delta \mathsf{K}+\Delta T={{R}_{energy}} \\
& {{\Upsilon }^{1}}={{\sum{\left( \frac{\partial {{e}^{0}}}{\partial T} \right)}}^{old}}\frac{\Delta \mathsf{K}}{{{\tau }_{U}}} \\
\end{align}

(13)
where $\backslash Delta\; e$ is $\{\{e\}^\{new\}\}\{\{e\}^\{old\}\}$ 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 $\{\{N\}\_\{K\}\}+1$, where $\{\{N\}\_\{K\}\}$ 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.