Prediction of thermal conductivity in dielectrics using fast, spectrally-resolved phonon transport simulations

TL;DR

This paper introduces a novel spectral phonon transport simulation method based on the Boltzmann equation to accurately predict the thermal conductivity of dielectric materials like silicon at various length scales.

Contribution

The paper develops a new computational approach using a Self-Adjoint Angular Flux formulation and Richardson iteration for spectral phonon transport simulations, enhancing prediction accuracy.

Findings

Accurate thermal conductivity predictions for silicon at room temperature.

Effective modeling of spectral phonon interactions and disorder effects.

Validation of the method across multiple length scales.

Abstract

We present a new method for predicting effective thermal conductivity ($\kappa_{\textrm{eff}}$) in materials, informed by ${ab\,initio}$ material property simulations. Using the Boltzmann transport equation in a Self-Adjoint Angular Flux formulation, we performed simulations in silicon at room temperatures over length scales varying from 10 nm to 10 $\mu$m and report temperature distributions, spectral heat flux and thermal conductivity. Our implementation utilizes a Richardson iteration on a modified version of the phonon scattering source. In this method, a closure term is introduced to the transport equation which acts as a redistribution kernel for the total energy bath of the system. This term is an effective indicator of the degree of disorder between the spectral phonon radiance and the angular phonon intensity of the transport system. We employ polarization, density of states and full dispersion spectra to resolve thermal conductivity with numerous angular and spatial discretizations.

Figures

Tables

Table 1. Table 1: Silicon material properties.

G	$\Lambda$ [nm]	$\overline{\omega }\cdot {10}^{13}\left[s^{-1}\right]$	$\overline{v}\left[\text{m}\cdot {\text{s}}^{-1}\right]$	$\overline{𝔻}\cdot {10}^{27}\left[{\text{m}}^{-3}\right]$	$\overline{\tau }\left[\text{ps}\right]$	$\overline{\omega v𝔻}\left[{\text{m}}^{-2}\cdot {\text{s}}^{-2}\right]$
1 (LA)	3120	1.49	8079	0.0736	386	$12.9\cdot {10}^{42}$
2 (LA)	155	4.25	6926	0.515	22.4	$156\cdot {10}^{42}$
3 (LA)	28.7	6.46	5152	1.4	5.6	$458\cdot {10}^{42}$
4 (TA)	898	0.99	5076	0.147	177	$9.4\cdot {10}^{42}$
5 (TA)	79.5	2.3	2063	1.03	38.5	$42.2\cdot {10}^{42}$
6 (TA)	11.1	2.75	721	2.8	15.4	$51.8\cdot {10}^{42}$
7 (O)	3.15	9.43	786	0.22	4	$24.9\cdot {10}^{42}$
8 (O)	3.8	8.92	1651	1.55	2.3	$219\cdot {10}^{42}$
9 (O)	4.58	8.43	1126	4.2	4.1	$397\cdot {10}^{42}$

Table 2. Table 2: Difference in κ eff subscript 𝜅 eff \kappa_{\textrm{eff}} with and without optic phonons for S 8 , F subscript 𝑆 8 𝐹 S_{8,F} simulation. Thermal conductivity is in units of [ W ⋅ m − 1 ⋅ K − 1 ] delimited-[] ⋅ W superscript m 1 superscript K 1 \left[\textrm{W}\cdot\textrm{m}^{-1}\cdot\textrm{K}^{-1}\right] .

$L$ [nm]	${\kappa }_{\text{LA,TA,O}}$	${\kappa }_{\text{LA,TA}}$	$\Delta \kappa (\%)$
10	6.1	$5.2$	14.8
100	30.3	$28.2$	6.9
1000	90.8	$88.4$	2.6
10000	168.6	$166.7$	1.1

Table 3. Table 3: S 8 , F subscript 𝑆 8 𝐹 S_{8,F} with LA, TA, O phonons.

$L$ [nm]	GMRES Iterations	$\min \left({\zeta }_{p,\eta }\right)$	$\max \left({\zeta }_{p,\eta }\right)$
10	29	$0.0032$	3.175
100	59	$0.032$	31.75
1000	1276	$0.32$	317.5
10000	46600	$3.2$	3175

Table 4. Table 4: S 8 , F subscript 𝑆 8 𝐹 S_{8,F} with LA and TA phonons.

$L$ [nm]	GMRES Iterations	$\min \left({\zeta }_{p,\eta }\right)$	$\max \left({\zeta }_{p,\eta }\right)$
10	28	$0.0032$	0.901
100	41	$0.032$	9.01
1000	114	$0.32$	90.1
10000	3306	$3.2$	901

Equations

$$\kappa =\frac{1}{8{\pi }^3}\sum _p{\int }_{\mathrm{BZ}}{\tau }_{p\mathbf{k}}{\mathbf{v}}_{p\mathbf{k}}\otimes {\mathbf{v}}_{p\mathbf{k}}{\omega }_{p\mathbf{k}}\hslash n\left({\omega }_{p\mathbf{k}},T\right)d{\mathbf{k}}^3,$$

$$\frac{df\left(\mathbf{r},\mathbf{k},p\right)}{dt}=-\mathbf{v}\left(\mathbf{k},p\right)\cdot \nabla f\left(\mathbf{r},\mathbf{k},p\right)+{\left[\frac{df}{dt}\right]}_{\text{coll.}}.$$

$$\left|\mathbf{v}\left(\hat{\mathbf{\Omega }},\eta ,p\right)\right|\hat{\mathbf{\Omega }}\cdot \nabla f\left(r,\hat{\mathbf{\Omega }},\eta ,p\right)={\left[\frac{df\left(r,\hat{\mathbf{\Omega }},\eta ,p\right)}{dt}\right]}_{\text{coll.}},$$

$$\hat{\mathbf{\Omega }}\cdot \nabla f=\frac{f^0\left(T\left(\mathbf{r}\right),\hat{\mathbf{\Omega }},\eta ,p\right)-f}{\Lambda },$$

$$f^0=f^0\left(T\left(\mathbf{r}\right),\hat{\mathbf{\Omega }},\eta ,p\right)=⟨n_{\text{BE}}⟩=\frac{1}{\exp \left[\frac{\hslash \omega \left(\hat{\mathbf{\Omega }},\eta ,p\right)}{k_{\text{B}}T\left(\mathbf{r}\right)}\right]-1}.$$

$$\psi \left(\mathbf{r},\hat{\mathbf{\Omega }},\eta ,p\right)=\hslash \omega \mathbf{v}f\mathbb{D}d\eta .$$

$$3N_v=\sum _p\int d\eta {\int }_{4\pi }\mathbb{D}\left(\hat{\mathbf{\Omega }},\eta ,p\right)d\Omega ,$$

$$\hat{\mathbf{\Omega }}\cdot \nabla \psi =\frac{{\varphi }^0-\psi }{\Lambda },$$

$${\varphi }^0={\varphi }^0\left(T\left(\mathbf{r}\right),\hat{\mathbf{\Omega }},\eta ,p\right)=\frac{\hslash \omega \mathbf{v}\mathbb{D}d\eta }{\exp \left[\frac{\hslash \omega }{k_{\text{B}}T\left(\mathbf{r}\right)}\right]-1},$$

$$U_T\left(\mathbf{r}\right)$$

$${\dot{U}}_I\left(\mathbf{r}\right)=\sum _p\int d\eta {\int }_{4\pi }d\Omega \frac{{\varphi }^0-\psi }{\Lambda }.$$

$$\hat{\mathbf{\Omega }}\cdot \nabla \psi =\frac{{\varphi }^0-\psi }{\Lambda }-\beta \frac{{\varphi }^0}{\left|\mathbf{v}\right|},$$

$$\beta =\beta \left(\mathbf{r}\right)=\frac{{\dot{U}}_I}{U_T}.$$

$$\hat{\mathbf{\Omega }}\cdot \nabla \psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)+\frac{1}{\Lambda }\psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)=\frac{1}{4\pi \Lambda }{\varphi }^{\text{T}}\left(\mathbf{r}\right),$$

$$-\hat{\mathbf{\Omega }}\cdot \nabla \left[\Lambda \hat{\mathbf{\Omega }}\cdot \nabla \psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)\right]+\frac{1}{\Lambda }\psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)=\frac{1}{4\pi }\left[\frac{1}{\Lambda }{\varphi }^{\text{T}}\left(\mathbf{r}\right)-\hat{\mathbf{\Omega }}\cdot \nabla {\varphi }^{\text{T}}\left(\mathbf{r}\right)\right].$$

$${\varphi }^{\text{T}}={\varphi }^{\text{T}}\left(\mathbf{r}\right)={\int }_{4\pi }\psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)d\Omega ,$$

$${\Phi }^{\text{T}}={\Phi }^{\text{T}}\left(\mathbf{r}\right)=\sum _p{\int }_{\eta }{\int }_{4\pi }\psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)d\Omega d\eta ,$$

$$\left|\mathbf{v}\right|\hat{\mathbf{\Omega }}\cdot \nabla \psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)=\frac{\frac{1}{4\pi }{\varphi }^0\left(T\left(\mathbf{r}\right)\right)-\psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)}{\tau }-\frac{\beta \left(T\left(\mathbf{r}\right)\right){\varphi }^0\left(T\left(\mathbf{r}\right)\right)}{4\pi },$$

$$\hat{\mathbf{\Omega }}\cdot \nabla \psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)+\frac{1}{\Lambda }\psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)=\frac{1}{4\pi \Lambda }{\varphi }^0\left(T\left(\mathbf{r}\right)\right)-\frac{\beta \left(T\left(\mathbf{r}\right)\right)}{4\pi \left|\mathbf{v}\right|}{\varphi }^0\left(T\left(\mathbf{r}\right)\right),$$

$$\psi \left({\mathbf{r}}_b\right)=\{\begin{array}{l}{\psi }^{\text{vac}}\left({\mathbf{r}}_b,\hat{\mathbf{\Omega }}\right),\hat{\mathbf{\Omega }}\cdot {\bar{\mathbf{n}}}_b<0\\ {\psi }^{\text{ref}}\left({\mathbf{r}}_b,{\hat{\mathbf{\Omega }}}_r\right),\hat{\mathbf{\Omega }}\cdot {\bar{\mathbf{n}}}_b<0\end{array},$$

$${\hat{\mathbf{\Omega }}}_r=\hat{\mathbf{\Omega }}-2\left(\hat{\mathbf{\Omega }}\cdot {\bar{\mathbf{n}}}_b\right){\bar{\mathbf{n}}}_b.$$

$$\psi =\frac{{\varphi }^0}{4\pi }-\frac{\Lambda \beta {\varphi }^0}{4\pi \left|\mathbf{v}\right|}-\Lambda \hat{\mathbf{\Omega }}\cdot \nabla \psi ,$$

$$\hat{\mathbf{\Omega }}\cdot \nabla \left(\frac{{\varphi }^0}{4\pi }-\frac{\Lambda \beta {\varphi }^0}{4\pi \left|\mathbf{v}\right|}-\Lambda \hat{\mathbf{\Omega }}\cdot \nabla \psi \right)+\frac{1}{\Lambda }\psi =\frac{{\varphi }^0}{4\pi \Lambda }-\frac{\beta {\varphi }^0}{4\pi \left|\mathbf{v}\right|},$$

$$-\hat{\mathbf{\Omega }}\cdot \nabla \left[\Lambda \hat{\mathbf{\Omega }}\cdot \nabla \psi \right]+\frac{1}{\Lambda }\psi =\frac{1}{4\pi }\left[\frac{{\varphi }^0}{\Lambda }-\hat{\mathbf{\Omega }}\cdot \nabla {\varphi }^0-\frac{\beta {\varphi }^0}{\left|\mathbf{v}\right|}+\hat{\mathbf{\Omega }}\cdot \nabla \frac{\Lambda \beta {\varphi }^0}{\left|\mathbf{v}\right|}\right].$$

$$f\left(\mathbf{r},\hat{\mathbf{\Omega }}\right)\to f\left(\mathbf{r},{\hat{\mathbf{\Omega }}}_m\right)\equiv f_m\left(\mathbf{r}\right),m=1,\dots ,M.$$

$$\psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)\approx \frac{1}{4\pi }\sum _{k=0}^{\infty }\sum _{n=-k}^{n=k}{\psi }_{k,n}\left(\mathbf{r}\right)Y_{k,n}\left(\hat{\mathbf{\Omega }}\right),$$

$${\psi }_{k,n}\left(\mathbf{r}\right)={\int }_{4\pi }{Y}_{k,n}^{\ast }\left(\hat{\mathbf{\Omega }}\right)\psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)d\Omega .$$

$${\varphi }_{l,n}\left(\mathbf{r}\right)=\frac{1}{4\pi }\sum _{l=0}^L\sum _{n=-l}^l\sum _{m=1}^M{w}_m{\psi }_m\left(\mathbf{r}\right)Y_{l,n}^{\ast }\left({\hat{\mathbf{\Omega }}}_m\right)Y_{l,n}\left({\hat{\mathbf{\Omega }}}_m\right),$$

$$\sum _{m=1}^M{w}_m=2\pi .$$

$$\mathbf{q}\left(\mathbf{r}\right)=\sum _p{\int }_{\eta }{\int }_{4\pi }\psi \left(\mathbf{r},\hat{\mathbf{\Omega }}\right)\hat{\mathbf{\Omega }}d\Omega d\eta ,$$

Full text

I Introduction

In many applications a material’s propensity for conducting heat is an important factor in its overall performance. For thermal management applications, having high or low thermal conductivity is of principal importance, but thermal conductivity makes an auxiliary contribution to a material’s figure of merit in many other applications such as batteries, high frequency electronics, and thermoelectric materials.

In order to accelerate the discovery, design, and selection of application-driven materials, it is desirable to efficiently compute a thermal conductivity from first principles. The past decade has seen considerable progress towards this goal in semiconductors and ceramics, where heat transport is dominated by phonons [1, 2]. Methods have been developed which make routine the computing of second order interatomic stiffness matrices of a crystal from density functional theory (DFT) calculations, and from them, the phonon dispersion [3, 4, 5, 6, 7, 8]. More recently, methods have been developed for computing third and fourth order stiffness matrices from which the rates of intrinsic three- and four-phonon scattering processes may be derived. The use of these methods is becoming commonplace, as together they enable one to compute the intrinsic phonon thermal conductivity, $\kappa$, using the equation derived from the Boltzmann transport equation for a phonon gas:

[TABLE]

which integrates the contributions to thermal transport from all phonons at all polarizations $p$ and wave vectors $\mathbf{k}$ in the crystal’s first Brillouin zone. Here, $\omega$ and $\mathbf{v}_{p,\mathbf{k}}$ are the phonon angular frequency and group velocity, respectively, which are obtained from the phonon dispersion relation.

The equation above gives a material’s inherent capacity for transporting heat under the assumption that the perturbation in the phonon population that gives rise to a net heat flux is spatially homogeneous. However, in most practical applications materials contain defects and impurities in addition to micro and nano-structured morphology; these all present additional resistance to phonon transport. Predicting the macroscale thermal conductivity in these heterogeneous materials is more complicated, but is of profound technological importance as the hierarchy of defect structure and morphology can be central to the material’s functionality. For example, in thermoelectric devices, dopants, impurities, inclusions, and nano-structuring are engineered into materials in order to optimize the electrical and thermal properties. In other applications where materials are under irradiation, it is essential to predict the evolution of thermal conductivity during the accumulation of radiation-induced damage.

In materials where atomistic to micro-scale heterogeneity impacts thermal conductivity, the phonon population is also not spatially uniform, and so predicting the macroscopic thermal conductivity requires an understanding of the full phonon population distribution within some representative microstructure. This can be done by simulating the spectrally-dependent phonon distribution as governed by the Boltzmann transport equation (BTE) using phonon transport properties derived from DFT calculations of atomic interactions. The BTE describes the evolution of a distribution of particles as a function of their position, frequency and direction of travel, making the BTE a six-dimensional equation (or seven-dimensional, in the time-dependent case).

Advances have been made in the simulation of spectrally-dependent phonon transport, providing important insights into the spectral contribution to heat transport, but to date the computational burden associated with BTE solutions have limited the size, dimensionality, and/or physics of the systems considered. Allu et al. [9] investigated frequency dependence of the BTE in silicon slabs and reported spectral thermal conductivity based on Knudsen number (the mean free path divided by the domain length). However, the trade-off of this approach was their reported high computational expense, and the discrete phonon groups had no coupling. Mazumder and Majumdar [10] demonstrated polarization dependence of phonons in Monte Carlo transport simulations and showed the importance of the branch dependence; polarization transitions occur for background scattering as well as anharmonic interactions [11]. Hua derived the time and frequency dependent phonon transport equation using the method of degenerate kernels in an attempt to increase computational efficiency compared to traditional integral discretization methods [12]. This was only done in 1-D, did not provide any spatial discretization, and was primarily analytical. As simulating frequency dependent phonon transport allows for the characterization of the frequency spectra, it allows for mode dependent contributions to the specific heat capacity and phonon velocity. Minnich et al. [13] demonstrated this mode dependence in their simulations of aluminum and silicon thin films, and calculated the spectral heat capacity from each frequency interval, explicitly illustrating the contribution of each mode to thermal conductivity. Rather than discriminating phonons by frequency and polarization, Romano has developed an approach that separates phonon modes into groups by discretizing the spectrum of phonon mean free path (MFP) [14, 15]. This provides considerable computational savings as the MFP is the only spectrally relevant phonon property that appears in the BTE. Most recently, Zhang et al. [16] developed a coupled approach using a reference temperature formulation which efficiently simulates temperature distribution and thermal conductivity over the nano- to micro-scale length range. However, they did not consider optical phonons with short mean free paths in their simulations. We discuss the results of Zhang et al. in greater detail below, but taken together, these prior works highlight that for simulating phonons it is imperative that the phonons are spectrally-resolved, and that however one defines the phonon spectrum, one cannot cavalierly disregard portions of it. This makes fully self-coupled calculations of large systems computationally challenging. In this manuscript we present an alternative approach for solving the spectrally-dependent phonon BTE that draws on numerical methods developed for neutron transport in the software package Rattlesnake [17]. The method solves the BTE in the self-adjoint angular flux (SAAF) formulation, which is computationally efficient, making it possible to simulate phonon transport in relatively large systems.

The numerical procedures for simulating the collective transport of phonons occupying a spectrum of vibrational modes is similar to the procedures used in nuclear engineering for simulating multi-group neutron or radiation transport [18, 19]. With neutron or radiation transport, the energy spectrum is partitioned into discrete groups and the transport simulation accounts for within-group and group-to-group scattering of neutrons based on the interaction cross sections. Rattlesnake [17] employs a multi-group energy discretization scheme to solves the BTE for neutrons. It is developed in the Multi-physics Object Oriented Simulation Environment (MOOSE) framework [20] at Idaho National Laboratory (INL) and solves the SAAF-BTE using a finite element spatial discretization and a discrete ordinates angular discretization. Having previously adapted Rattlesnake to model one-group (grey) phonon transport [21, 22, 23] in both homogeneous and heterogeneous media, in the work presented here we have further developed Rattlesnake to simulate spectrally-dependent phonon transport.

The remainder of the manuscript is organized as follows: In the next section we develop the methods used to derive the SAAF-BTE for spectral transport, temperature coupling between discrete phonon groups, and the energy redistribution function used in closing the system of equations. Following this, we describe the methods used to compute a material’s phonon properties and the methods we have developed to convolve these data into the collective properties of a set of discrete transport groups simulated in the SAAF-BTE. Finally, to demonstrate the computational efficiency of our multi-group, temperature coupled approach, simulations are presented of phonon transport in silicon of varying thickness. For these we report spectrally-resolved temperature, heat flux, and thermal conductivity. We discuss each of these quantities, how they are affected by spatial and angular resolution, the effect of the new closure term, and how convergence is affected by the simulation of all phonon modes.

II Methods

The generalized Boltzmann transport equation is used widely by the transport community to model the evolution of a distribution of particles or carriers. The modes of atomic vibration of crystal are traveling waves with well defined wave vectors. These waves extend across the entirety of the crystal; however, if uncertainty is applied to the wavevector the vibrations become propagating localized wavepackets. As the energy of vibrations comes in discrete quanta proportional to their vibrational frequency, individual wavepackets can be treated as quasiparticles – phonons – that collectively follow Bose-Einstein statistics. The total vibrational energy of the crystal can thus be modeled as a phonon quasiparticle gas and so can be described by the Boltzmann transport equation for uncharged particles:

[TABLE]

Here $f\left(\mathbf{r},\mathbf{k},p\right)$ is the number of phonons at location $\mathbf{r}$ occupying the phonon mode with wavevector $\mathbf{k}$ and polarization $p$. The term $\mathbf{v}$ is the group velocity of a phonon wavepacket and is given by $\mathbf{v}=\mathbf{v}\left(\mathbf{k},p\right)=\nabla_{\mathbf{k}}\omega\left(\mathbf{k}\right)$, the derivative of the vibration’s angular frequency, $\omega$, with respect to $\mathbf{k}$. The group velocity is a vector normal to the isofrequency surface for the dispersion relation which is not necessarily parallel to $\mathbf{k}$. While the natural parameter for enumerating the phonon modes of a crystal is the wavevector and polarization (4 quantum numbers are required), for transport simulations it is more natural to parameterize the particle/carrier population by the particles’ direction of travel, and then by their speed, and to treat these two parameters separately. This is particularly true for neutron transport as their dispersion relation is isotropic. To put the phonon transport problem in to a form where it can leverage the numerical power of existing transport solvers such as Rattlesnake, and restricting our interest to the time independent case, in what follows we write the Boltzmann transport

[TABLE]

where $f=f(\mathbf{r},\hat{\mathbf{\Omega}},\eta,p)$ has been parameterized by the phonons’ direction of propagation, $\mathbf{\hat{\Omega}}$, and spectral variable $\eta$. The spectral metric $\eta$ could be energy, frequency, wave number, mean free path, or any other descriptor or quantum number for enumerating the phonon modes.

The term on the right hand side (RHS) of Eq. (3) is the rate of change of the population of due to the creation/annihilation of phonons and events that change a phonon’s direction of travel. We model this using the single mode relaxation time (SMRT) approximation [24] which assumes that the net rate of loss/gain of phonons in a single mode is described by a single relaxation time, $\tau=\tau(\hat{\mathbf{\Omega}},\eta,p)$, and is proportional to the deviation of the population from thermodynamic equilibrium. Defining the mean free path (MFP) of phonons $\Lambda(\hat{\mathbf{\Omega}},\eta,p)=\left|\mathbf{v}\right|\tau$, the simplified BTE is written as

[TABLE]

where the equilibrium phonon distribution $f^{0}$ is described by Bose-Einstein statistics

[TABLE]

Here $T\left(\mathbf{r}\right)$ is the local temperature, $\hbar$ is the reduced Planck constant and $k_{B}$ is Boltzmann’s constant.

The population $f$ gives the number of phonons in the mode characterized by $\mathbf{\hat{\Omega}}$, $\eta$, and $p$ and is thus dimensionless. As we are interested transport of heat rather than phonon number is more useful to weight the population by the modes’ contribution to energy flux, and solve for the angular flux intensity, $\psi$, which has dimensions of power per area per steradian, and is defined as

[TABLE]

Here the term, $\mathbb{D}=\mathbb{D}(\hat{\mathbf{\Omega}},\eta,p)$ is the phonon density of states (DOS) per increment of solid angle $\mathbf{\hat{\Omega}}$ and spectral descriptor. It has dimensions of number per unit volume per steradian per $\eta$, and it satisfies the relationship

[TABLE]

where $N_{v}$ is the number of of atoms per unit volume. The BTE in terms of angular radiance is

[TABLE]

where $\phi^{0}$ is the radiance of the mode at equilibrium given by

[TABLE]

and is discretized into partitions with dependency on $\hat{\mathbf{\Omega}}$ within the Brillouin zone; this method is discussed in a proceeding section. The local temperature is defined as the temperature that a system at equilibrium would have in order to have the same thermal energy, $U_{T}$, as the phonon bath

[TABLE]

The numerical approach to extract $T(\mathbf{r})$ from this relationship is described in detail later in the manuscript, but we first address the issue of energy conservation.

The SMRT approximation for the collision term that appears in the BTE in Eqs. (4) and (7) assumes that modes which are out of equilibrium exchange energy with an implied phonon reservoir at the local temperature to relax them back into equilibrium. However, as the rate of this energy exchange is different for all modes for any arbitrary phonon population, it is likely that there will be a net flow of energy to or from the phonon bath, resulting in a lack of local energy conservation. The problem is that the flow of energy into the implied phonon bath in the SMRT is not reflected in the explicit phonon bath used to compute temperature in Eq. (9), and so we must include an additional source term to correct this.

The net energy flow per unit volume, $\dot{U}_{I}$, to the implied phonon bath is

[TABLE]

To conserve energy, we must return this energy to the explicit phonon bath, and to do that we must first decide how to distribute the energy over the phonon modes. The most obvious choice for this is to add the energy back proportional to the equilibrium distribution of energy across the modes. Applying this correction to the BTE to close the coupling of SMRT with the phonon bath gives:

[TABLE]

where $\beta$ is a rate and is given by

[TABLE]

This residual energy projection is a measure of the balance of the equilibrium phonon population to the transport system population; $\beta$ is tightly coupled to temperature, the transport and equilibrium systems, and is the closure required to ensure energy conservation. In our framework, both the temperature and $\beta$ are functionals of the entire phonon population. The energy flow $\dot{U}_{I}$ can be either positive or negative depending on the phonon population. The most ballistic groups tend to have the strongest non-equilibrium condition occurring near isothermal boundaries and these are the regions in which we find $\beta$ to be the largest.

II.1 SAAF-BTE for phonons

In our previous work [21, 22, 23], we modified the SAAF-BTE for transport of grey neutrons to make it treat grey phonon transport, and then solved it to model heat conduction in UO2 containing Xe bubbles. In grey transport the particles of interest are resolved in direction but the spectral properties are averaged over the entire distribution making a “grey” spectrum. In spectrally-resolved transport the grey approximation is equivalent to coupling the particle distribution evolution over angle, but having no coupling between the different portions of the property spectrum as in the first and second order equations

[TABLE]

[TABLE]

Here the coupling term, $\phi^{\textrm{T}}$, is the zeroth angular moment

[TABLE]

which we refer to as the transport scalar flux. The total scalar flux is defined as the integral of the zeroth angular moment over all spectral groups and polarizations

[TABLE]

and $\phi^{\textrm{T}}$ and $\Phi^{\textrm{T}}$ both have units of $\textrm{W}\cdot\textrm{rad}^{-1}\cdot\textrm{m}^{-2}$. In the BTE, the phonon angular intensity is defined along a direction of travel $\hat{\mathbf{\Omega}}$, a function of polar and azimuthal angles $\theta$ and $\phi$; $\hat{\mathbf{\Omega}}$ is an independent variable of $\psi$, which itself is the dependent variable, $\psi(\mathbf{r},\hat{\mathbf{\Omega}})$.

The weak form of Eq. (12) was then solved with Rattlesnake. In this current work, we introduce a new form of the SAAF-BTE specifically for phonons, which relies on implicit temperature coupling to better represent the subtle physics of phonon transport. We return to Eq. (4), with phonon state variables multiplied through and the addition of the non-equilibrium bath source correction to the right hand side

[TABLE]

and rearrange

[TABLE]

where the terms on the left are the streaming and collision operators, respectively, with source terms on the right hand side. We have included temperature dependent notation $\left(T\left(\mathbf{r}\right)\right)$ in Eq. (14) as a reference point to show the terms where temperature is coupled in our approach. For much of the remainder of this section, we drop the independent variables $(\mathbf{r},\hat{\mathbf{\Omega}})$.

The vacuum and reflecting boundary conditions for this transport equation are defined as:

[TABLE]

here, $\bar{\mathbf{n}}_{b}$ is the outward unit normal at a point $\mathbf{r}_{b}$ on the boundary. In neutron transport space, $\psi^{\textrm{vac}}$ implies a vacuum boundary ($\psi(\mathbf{r}_{b},\hat{\mathbf{\Omega}})=0$). However, in this implementation, we leverage the vacuum boundary as an adiabatic boundary to specify an incident source of phonons, e.g., $\psi(\mathbf{r}_{b},\hat{\mathbf{\Omega}})=\phi^{0}\left(\mathbf{r}_{b}\right)$. The reflective angle $\hat{\mathbf{\Omega}}_{r}$ in $\psi^{\textrm{ref}}$ is

[TABLE]

The derivation of the SAAF form of the phonon BTE is based on a straightforward algebraic technique. From a computational perspective the SAAF formulation is advantageous as the full angular flux intensity is the unknown. Upwinding is a common numerical discretization of advection operators in partial differential equations, and is frequently used with the traditional first-order form of the Boltzmann transport equation, particularly in one and two spatial dimensions. The (primarily) lower triangular linear algebra problem in each quadrature direction is solved via a “transport sweep” - forward propagation of information through spatial cells in the direction of particle travel. The SAAF equations have an elliptic streaming term, and the resulting linear algebraic equations for the intensity in each angular ordinate are then solved via preconditioned GMRES methods on unstructured spatial grids – this is the primary numerical advantage to using the SAAF equations; boomerAMG-preconditioned GMRES is massively parallelizable. In contrast, sweeping algorithms (especially on unstructured grids in 3D) are challenging to implement in parallel; the development of an efficient parallel sweep algorithm is an active area of interest in the neutron transport community [25]. The use of reflecting boundary conditions is easier with the availability of the full angular flux, as the incoming and outgoing directions are coupled in the same manner as the first order form of the transport equation [26]. Through the application of a continuous finite element (CFEM) spatial discretization, the matrices are symmetric positive definite (SPD), which allows for the use of solution techniques such as the preconditioned Krylov family of solvers. To obtain the SAAF form of the BTE, a simple algebraic approach is followed, solving Eq. (14) for the angular intensity:

[TABLE]

and substitute this expression back into the streaming term in Eq. (14)

[TABLE]

rearranging yields the SAAF form of the phonon BTE

[TABLE]

We use the discrete ordinates method [27] to discretize the angular variable. In the discrete ordinates approach, we represent the independent variable $\hat{\mathbf{\Omega}}$ by a discrete set of directions, $\{\hat{\mathbf{\Omega}}_{m},m=1,\ldots,M\}$. The functions of $\hat{\mathbf{\Omega}}$ are represented only by their values at each of their directions on the spatial mesh, e.g.,

[TABLE]

Consider the angular flux intensity as a function of the direction variable $\hat{\mathbf{\Omega}}$. There exists a set of complete orthogonal functions in $\hat{\mathbf{\Omega}}$, much like Legendre polynomials; these are the spherical harmonics functions. To compute the total scalar flux in a group at the end of an iteration, the angular flux intensity is integrated over direction using Eq. (13). However, we make the approximation that the angular flux intensity $\psi(\mathbf{r},\hat{\mathbf{\Omega}})$ can be represented by these spherical harmonics functions (termed $Y$) and we expand the angular flux in terms of these functions

[TABLE]

where

[TABLE]

This approach involves expanding the angular dependence of the flux intensity in a finite series of spherical harmonics $Y_{l,m}(\hat{\mathbf{\Omega}})=Y_{l,m}\left(\theta,\varphi\right)$, a familiar sight in quantum mechanics. Coupling for the angular flux intensity and the transport flux exists in the angular moments alone

[TABLE]

here, $L$ is the truncated spherical harmonics, $l$ is the degree of spherical harmonics and $n$ is the order, and weights sum to $2\pi$ in 2D

[TABLE]

The total heat flux $\mathbf{q}\left(\mathbf{r}\right)$ is the sum of the first angular moments ($\phi_{1,0}$ and $\phi_{1,1}$ for two spatial dimensions) over all spectral groups and polarizations

[TABLE]

and has units of $\textrm{W}\cdot\textrm{m}^{-2}$. In this way we solve the transport equation for each angular direction $\hat{\mathbf{\Omega}}_{m}$. At the end of each iteration, the numerical integration is performed to represent the scalar flux with these angular moments. Then, $\phi^{\textrm{T}}$ is used to determine the local temperatures in the transport system, and are then passed into the equilibrium radiance relations. Now that we have represented the angular flux intensity by this quadrature, the discrete form of Eq. (16) for any direction $\hat{\mathbf{\Omega}}_{m}$ with an associated weight $w_{m}$ in the angular quadrature set is:

[TABLE]

where we have assigned iteration indices to Eq. (20); we solve for the value of angular flux intensity, $\psi^{\left(\ell+1\right)}$ using sources from the previous iteration $\left(\ell\right)$. This resembles source iteration, a well known solution technique in neutron transport [27, 28, 29]. This procedure is outlined in Algorithm 1.

The transport calculation assumes that all of the spectral groups are decoupled; each system of equations is solved separately and no information is passed between them. We have demonstrated this capability previously with Rattlesnake [22, 23]. It is possible to solve each system of equations separately and then compute an average material temperature for the domain. However, this only considers isolated group contributions to temperature and heat flux. Explicit simulation of group-to-group phonon scattering foregoes the SMRT approximation, is computationally demanding, and is an active research area [13, 30]; this work does not consider explicit scattering processes.

The physics dictate that phonon modes are coupled as they become excited through phonon collisions and temperature changes. The amount of heat flowing through each phonon mode in the system must ‘see’ the others – even forsaking an explicit collision model, individual modes must feel the effects of the material temperature; Bose-Einstein statistics assert this, and the phonon groups must be globally coupled in some manner. A material temperature is determined by prescribing a temperature gradient between spatial boundaries, giving rise to blackbody phonon emission into the material. Phonon scattering in each of the groups are driven by this temperature gradient, but again no coupling exists. This decoupled approach fails to consider the temperature in the Bose-Einstein statistics which must be extended to all phonon modes – the true equilibrium distribution is an integration coupling all spectral zones. Temperature is defined by the phonon phase space distributions, spanning all modes and polarizations. Phonons are coupled via the internal energy of the crystal system, and we can leverage this detail to derive a set of self-consistent equations where an average material temperature is shared between each of the discrete phonon groups.

The total phonon radiance of the system is defined as

[TABLE]

which is just the contribution of $\phi^{0}$ over all modes. We make the assumption that the total equilibrium radiance $\Phi^{00}$ can be related to the total transport scalar flux in the domain, such that the two must be in radiative equilibrium. This treatment uses the Bose-Einstein distribution as the phonon source for each of the spectral groups as opposed to the transport scalar flux; we still perform an integration of the angular intensity at the end of each linear iteration to obtain transport scalar flux. Equating $\Phi^{00}$ to the sum of the transport flux, we have

[TABLE]

where the $\left(p,\eta\right)$ index denotes a double summation over the polarizations and spectral property indices. Now we may carry out the solution for $T\left(\mathbf{r}\right)$, which involves summations over the angularly dependent material properties (index $m$) and group index $g$, both described in proceeding Sec. II.2.

[TABLE]

Equation (23) is used to compute an average temperature for the entire spatial domain. $T\left(\mathbf{r}\right)$ is what couples the phonon groups together – they now depend on a global average temperature. The resulting value of $T\left(\mathbf{r}\right)$ is used to compute a new value of $\phi^{0}$ using Eq. (8).

Equation (20) in the discrete ordinate method generates a linear system of equations arising from the spatial discretization of an elliptic operator, and is solved for angular intensity $\psi_{p,\eta,m}$ in each quadrature direction. This means that software for solving the diffusion approximation to transport can be exploited to implement acceleration techniques, which will be part of a future study. The Galerkin weak form of Eq. (20) using weighted residual formulation with CFEM is then solved with Rattlesnake. The derivation of the weak form of the SAAF-BTE for phonons is similar to that of the weak form of the SAAF-BTE for neutrons and can be followed in the Rattlesnake theory [17]. Algorithm 1 is outlines the technique used to solve the system of equations.

II.2 Phonon Transport Properties of Silicon

For the purposes of demonstrating and testing the use of Rattlesnake to simulate spectrally-resolved phonon transport we use silicon as a model material system. The thermal conductivity of Si is isotropic and well characterized, and there are a wealth of experimental and computational data against which to compare our model. Silicon has cubic crystal symmetry; to simplify our model and to enable us to use a standard set of transport ordinates akin to neutron transport (which does not have cubic symmetry), we simplify the transport spectrum to make it spherically symmetric. The discussion that follows describes first the process we used to compute ab initio the phonon properties of Si, and the metric we have developed for grouping the phonon modes of Si into a discrete set of transport groups with aggregated transport properties.

II.2.1 Calculation of phonon properties

The frequency and group velocity of phonons with a given wave vector and polarization are computed from the dispersion relation obtained from the dynamical matrix of second order force constants. The probability of three phonon scattering processes, the dominant intrinsic scattering process at high temperature, are computed from the matrix of third order stiffness interactions. The \nth2 and \nth3 order force constants for Si at 0 K were obtained from the AlmaBTE database [5]. In this data set the \nth2 order force constants were computed via the frozen phonon method as implemented using phonopy [7] using a large $5\times 5\times 5$ supercell in order to accurately capture the contributions from long ranged interatomic interactions. The \nth3 order force constants were calculated using a $5\times 5\times 5$ supercell using the package thirdorder.py [6] and considered all three-atom interactions of atoms less than five neighbours removed from one another. The complexity of the unit cell does not directly affect computational cost, but it can have indirect influence on the computational cost. The main source of computational expense comes from how finely one resolves the phonon spectrum, which depends on the span of frequencies and mean free paths. However, the CPU time is impacted by the acoustic thickness of a particular angle-group; if the unit cell has many optical modes with low velocity and short lifetime, these diffusive modes will take longer to converge. We discuss this further in a proceeding section. The full details of the DFT calculations are given in reference [6].

Starting with these stiffness matrices we used AlmaBTE to compute the full phonon dispersion and group velocities on a $24\times 24\times 24$ Brillouin zone mesh. The 0 K velocities and frequencies were assumed to be unchanged at higher temperatures. The phonon-phonon scattering times were computed on the same mesh at a temperature of 300 K. All of the phonon modes in the Brillouin zone are represented in the transport groups, albeit in an approximate way. In this work, this was done to make the transport spectrum spherically symmetric; however it is no extra work to include all the modes in the Brillouin zone in a non-approximate way, and to model anisotropic materials. For this, one would need to use a set of ordinates which match the crystal symmetry, and we would not be able to decouple the transport group’s properties and directions.

II.2.2 Decomposition of phonon spectrum into transport groups

Discretization of the BTE in Eq. (20) means solving for the heat flux in a set of transport groups where each group represents the combined contribution to transport from a collection of phonon modes with similar direction of propagation and spectral character. The task is to determine the effective transport properties of each group such that all phonon modes are accounted for, and that the temperature in Eq. (23) may be inverted.

In discretized form, the angular and spectral distribution of a quantity $A(\hat{\mathbf{\Omega}},\eta,p)$ of the population is approximated as a sum of piece-wise basis functions:

[TABLE]

Here $m$ is the index over discrete ordinates in $\hat{\mathbf{\Omega}}$, and $g$ is the index over control points (groups) in $p$ and $\eta$. The functions $\chi_{m}(\hat{\mathbf{\Omega}})$ and $\chi_{g}\left(\eta\right)$ are element shape functions — non overlapping boxcar functions that satisfy the relationships

[TABLE]

where $\Delta\eta_{g}$ is the width of the $g$-th spectral group, and

[TABLE]

The transport group quantities $A_{m,g}$ represent the combined contribution to transport from all of the phonon modes that fall into the $m,g$-th group. To perform this calculation we define the group integration operator

[TABLE]

Using this operator we define the groups’ mean free path as the average of the mean free paths of the phonon modes in the group weighted by their contribution to the phonon flux

[TABLE]

Where $\overline{T}$ is the expected average temperature of the system to be simulated. Similarly, for computing radiance from temperature we define:

[TABLE]

The only material properties in the discretized BTE (Eq. (20)) are the groups’ effective mean free path and radiance, however the transport simulation needs to be able to compute the local temperature and so it also requires the weighting terms to obtain this from the sum over fluxes. To compute the energy in a group from the angular flux, we define the group velocity as the averaged velocity weight by mode energy

[TABLE]

For the local temperature calculation we need to know the total of the product of frequency, density of states and velocity

[TABLE]

and the average frequency

[TABLE]

The approach developed up to this point is quite general; we have not yet specified a choice of spectral parameterization, and the method can be applied to anisotropic materials if one uses a set of transport ordinates with the appropriate symmetry. In the remainder of the manuscript we limit ourselves to one choice for parameterization of the dispersion in Si, in which we approximate the dispersion relation to be spherically symmetric enabling us to use a standard set of non-symmetric transport ordinates. In this approach we used just the dispersion for Si along the (100) crystal direction, and we collapsed the three optical branches and the two transverse acoustic branches into a single branch with averaged frequency as shown in Fig. 1. The dispersion was parameterized in wave number, $k$, and this condensed (100) dispersion was used as the dispersion along all directions. The density of states was defined as

[TABLE]

where $d_{p}$ is the degeneracy that accounts for the merged branches and $a_{\mathrm{Si}}$ is the lattice parameter of the Si unit cell. The integral of this accounts for all the phonon modes of Si

[TABLE]

The approximated dispersion was split into between three and 20 groups in $k$ of equal size.

Using these data, we can compute radiant emission sources on the boundaries, driven by a prescribed temperature boundary condition,

[TABLE]

Equation (24) is used by the vacuum condition at a boundary $\mathbf{r}_{b}$ to supply an incident source of phonons for a given group and polarization. We compute an effective thermal conductivity along a direction by taking the ratio of the total heat flux to the end-to-end temperature gradient (which includes boundary effects) in the system

[TABLE]

III Results

We performed spectrally resolved phonon transport simulations in two-dimensional planes of silicon of varying geometric size. We report heat flux, thermal conductivity, and the equilibrium temperature distribution. We used full phonon dispersion and density of states computed at room temperature from ab initio DFT simulations. Spatial domain sizes varied from 10 nm to 10 $\mu$m and were spatially discretized using coarse (C) and fine (F) triangular finite element meshes (Fig. 3). We employed $S_{4}$, $S_{8}$, and $S_{16}$ Gauss-Legendre angular quadratures. We simulate a 1 K temperature gradient along the $x$-axis, with boundary temperatures of $T_{\textrm{L}}=301\,\text{K}$, $T_{\textrm{R}}=300\,\text{K}$. Reflecting conditions are placed on remaining boundaries. We use AMG-preconditioned GMRES [31, 32, 33] to solve the linear system of equations, with convergence criteria set to $\epsilon=10^{-8}$. The selected four cases we discuss in this section were all simulated using $S_{8}$ quadrature with the fine spatial mesh.

In centerline temperature profiles (Fig. 4), it is clear that non-equilibrium behavior arises when the domain is small, and this occurs because radiative equilibrium cannot be established. The incident phonon radiance from opposing sides encounter interference due to their proximity – the distance between hot and cold sources is less than the mean free path of the majority of phonons. This feature is prominent in smaller geometric domain sizes, exacerbated by the presence of ballistic phonons, which may undergo very few collisions before reaching the opposite side of the domain. An equilibrium solution to the temperature distribution exists upon numerical convergence of the simulation as the temperature distribution in each phonon mode is identical. However, if an equivalent transient simulation was conducted, the modal temperature distributions may not be identical. The temperature profiles produced from this work may provide a basis for the benchmarking of temperature distributions in molecular dynamics (MD) simulations, which rely on the careful selection of fitting parameters; the temperature profiles computed by our method may aid in the accurate fitting of these parameters [34].

The closure term $\beta$ (Fig. 5) provides a glimpse into the non-equilibrium behavior of the simulation at various domain sizes; $\beta$ tightly couples the total phonon radiance ($\Phi^{00}$) and total flux ($\Phi^{\textrm{T}}$). We have shown $\beta$ multiplied by the average relaxation time $\overline{\tau}$ in order to understand the fraction of total energy exchanged into the bath in one scattering event. The total amount of energy exchanged in $\beta$ between $\Phi^{\textrm{T}}$ and $\Phi^{00}$ is infinitesimal relative to the total energy of the system, but without the presence of $\beta$ conservation is broken, and a parabolic, rather than constant heat flux is observed.

$\Phi^{\textrm{T}}$ and $\Phi^{00}$ are not equivalent at small domain sizes; this occurs due to the disparity between the boundary emission source and the localized transport flux, and it is no surprise that the majority of the curvature exhibited by $\beta$ exists near the boundaries of Fig. 5. This artifact is a direct consequence of the system size, where emitting boundaries are in competition with each other. $\beta$ exhibits a sign change, beginning negative at $T_{\textrm{L}}$ and turning positive at $T_{\textrm{R}}$. In cases where $\beta<0$, $\phi^{\textrm{T}}>\phi^{0}$; $\phi^{\textrm{T}}$ is strongly influenced by the boundary conditions, where the incident angular intensity $\psi$ is specified by the phonon radiance relation in Eq. (24). As phonons flow from hot to cold sides, the sign of $\beta$ tends positive as $\phi^{0}>\phi^{\textrm{T}}$, after the spatial midpoint. This is because the local radiance $\phi^{0}$ is now stronger than local transport flux $\phi^{\textrm{T}}$ – the cold boundary can not overcome the flow of phonons from the hot side, which results in a positive $\beta$, but also elevates the temperature away from the cold boundary shown on the profiles in Fig. 4. Blackbody phonon emission is very strong in proportion to the transport flux, and although overall convergence is achieved the localized non-equilibrium effect is very apparent, especially in smaller domain sizes.

Centerline total heat flux profiles for the $S_{8,F}$ case are shown in Fig. 6. As system size increases, heat flux decreases. Because of the flat heat flux profiles observed in all cases, it is clear that including $\beta$ to close the transport equation is a necessity. Without $\beta$, energy leaks out of the system and yields a non-conserved heat flux. However, inspecting group heat flux paints a clear image of which groups carry the most heat in the system. It is, of course, the groups with the largest $\Lambda$, and those which remain the most ballistic over the entire geometric domain range. Figure 7 shows group heat flux for the $S_{8,F}$ case for all domain sizes. It is clear that larger values of $\Lambda_{p,\eta}$ are responsible for a higher heat flux; ballistic phonons carry energy further between collisions. The LA phonons are the dominant carriers due to their higher velocities and mean free path whereas TA and O phonons contribute lesser quantities. Acoustic phonons are generated by atomic displacements moving in phase, which is responsible for their higher velocity. Conversely, optic phonon motion is out of phase and are shorter range carriers. This phenomenon can be traced to Fig. 2, where larger magnitudes of the derivatives of the dispersion curves equate to higher wave propagation speed. Heat flux flattens out as geometric size increases, shown in Fig. 7; in this way we can analyze the degree to which the heat-carrying groups are affected by the change in domain size. If phonon groups were decoupled, the relation between heat flux and domain length would remain fixed in each group, as individual groups are not influenced by the average energy of the domain. However, in a coupled simulation the equilibrium distribution of each phonon channel is impacted by the average material temperature. Results for heat flux and thermal conductivity vary by approximately 3% between all angular and spatial resolutions. Thus, for a homogeneous material, a moderate angular and spatial resolution is sufficient for accurate results. This will likely change in heterogeneous environments such as porous materials, bulk material with dopants or inclusions (e.g. UO2 with Xe bubbles), where ray effects have been observed in grey simulations [23].

Normalized spectral heat flux is shown in Fig. 8, giving insight into how rapidly each of the groups approach the ballistic limit. It is clear that $q_{1}$ (with $\Lambda_{1}=3120$ nm) carries the dominant portion of energy throughout the entire range of domain sizes, as it barely begins to approach the asymptotic limit at a length of 10 microns. However, it is overtaken early on by groups 2-5, until the domain sizes reaches about 1000 nm. As expected, the diffuse groups do not contribute appreciably to the overall heat flux.

Thermal conductivity for all cases is shown in Fig 9. The increasing angular refinements reduce heat flux and thermal conductivity by minute amounts, caused by additional information present due to a higher number of streaming directions. The sizes of the domains span 10 nm to 10 $\mu$m and $\kappa_{\textrm{eff}}$ is proportional to this size. As domain size increases, $\kappa_{\textrm{eff}}$ asymptotically approaches a bulk value of about $165\,\textrm{W}\cdot\textrm{m}^{-1}\cdot\textrm{K}^{-1}$, and many researchers have reported values for $\kappa_{\textrm{eff}}$ with variation on the order of 100% [35]. As spatial and angular discretizations change, so do the values of spectral $\kappa_{\textrm{eff}}$. On smaller spatial meshes, $\Lambda_{p,\eta}$ tends to be ballistic per cell, even if they would be diffuse on the global domain. The coarse spatial discretization is acoustically thicker (more diffuse) where the fine spatial discretization is acoustically thinner (more ballistic). It is only as global domain increases in size that these groups become acutely diffuse, and their contribution to total heat flux becomes quite negligible; this trend is shown in Figs. 7 and 8. Table 2 contains the contributions from optic phonons to $\kappa_{\textrm{eff}}$, and they have a stronger effect in smaller global domains. For simulations with large (greater than 1 $\mu$m) domain sizes, omitting optic phonons would cause a dramatic decrease in the amount of iterations required to converge the solution, with negligible loss of accuracy (Table 2). Overall, it is more important to finely resolve spatial cell (finite element) size in comparison to angular resolution.

Comparisons of $\Phi^{00}$ and $\Phi^{\textrm{T}}$ are shown in Fig. 10. The distance between $\Phi^{00}$ and $\Phi^{\textrm{T}}$ are related to the magnitude of $\beta$, and intuitively, as domain size increases and the magnitude of $\beta$ decreases, $\Phi^{00}$ and $\Phi^{\textrm{T}}$ should converge. However, for the 100 nm domain this is not the case; the two fluxes are further apart. Looking back to Fig. 7, there is a slight bump for groups 1 and 4, where heat flux experiences a slight rises rather than a monotonic decrease. Modal dominance of heat flux has an affect on the magnitude between $\Phi^{00}$ and $\Phi^{\textrm{T}}$ at 100 nm, where the system is still sufficiently far from equilibrium. The same effect is present in domains up to 600 nm, after which the heat flux contribution from group 4 at 700 nm finally becomes less than the heat flux at 50 nm. This lag is caused by phonons in group 1 finally overtaking those from group 4; this transition happens in between 600 and 700 nm, and is observed in Fig. 8. As the geometric domain increases to the micrometer scale, $\Phi^{00}$ and $\Phi^{\textrm{T}}$ become increasingly convergent and their respective difference reduces to $<1\%$. These criteria lead to the system being thrust into a state of equilibrium, compared to differences in the distributions on the nanometer scale where the system is far from equilibrium.

Phonon transport suffers from boundary scattering effects, and in systems which have strong non-equilibrium behavior, the boundaries can have a profound influence on solution behavior a finite distance into the geometric domain. In our method, $\beta$ has a strong gradient near any emitting boundary; both the equilibrium radiance and transport scalar flux are rapidly changing a finite distance away from a boundary. The magnitude of $\Phi^{\textrm{T}}$ near a boundary is a direct result of the boundary emission, where each outgoing ordinate of the angular intensity at an emissive boundary is assigned an incident phonon radiance (the adiabatic boundary condition). In turn, the magnitude of $\Phi^{00}$ depends on the solution of $\Phi^{\textrm{T}}$, but will not be equivalent near a boundary at short geometric domains due to the influence of an opposing boundary emission. To combat the strong boundary effects, we applied adaptive mesh refinement (AMR) using a gradient jump indicator. The gradient jump criteria takes the gradient norm of any specified solution quantity (in this case, we chose $\beta$) at the face between a mesh element and its’ neighbor element. The elements are sorted by increasing error, and the elements which fall into a specified tolerance are refined. As we used triangular finite elements, a single level of refinement splits one element into four smaller elements of equivalent area, and so on; we specified a maximum of two refinement levels. All refinement took place a finite distance from the emissive boundaries, shown in Fig. 11, where the gradient of $\beta$ changes rapidly. The strongest boundary effects are observed in the 10 nm and 100 nm simulations, which is expected as phonons will nearly all be in the ballistic scattering regime at those domain sizes.

The planar temperature distribution for four length scale (10 nm, 100 nm, 1 $\mu$m, 10 $\mu$m) cases are shown in Fig. 12; domain size increases from (a) - (d). In small domain sizes, equilibrium is difficult to establish. Phonon flux on the hot side is suppressed by phonons incident from the cold side, which reduces temperature. As the geometric domain increases in size, the Fourier limit is recovered and heat transport exists in an almost purely diffusive regime.

In Figure 13 the distribution of $\beta$ is shown for various domain sizes. The relationship between a system unable to establish equilibrium (Fig. 13a) to one in an equilibrate state (Fig. 13d) is stark – $\beta$ is nearly constant in large systems as the difference between $\Phi^{00}$ and $\Phi^{\textrm{T}}$ reduces. $\beta$ also shows the effect of the fixed boundary temperature and the slight curvature effect present in the temperature distribution. The magnitude and sign of $\beta$ also affect the convergence of the transport problem. Phonon transport simulations are classified as “purely scatterering” meaning the ratio of the mean free paths on the left and right hand sides of Eq. (20) are unity, and these physics cause convergence to lag in acoustically thick problems. The presence of $\beta$ affects the balance of Eq. (20), and in simulations where $\beta$ is not approximately constant in the bulk, the disparity between $\Phi^{00}$ and $\Phi^{\textrm{T}}$ is also large, in correlation with Fig. 10. This energy redistribution was required to close the transport system – without its presence, we observed a non-conserved heat flux. The inclusion of $\beta$ projects the residual energy; the phonon groups are analogous to networks of siphons drawing from a central plenum, the energy bath.

The ratio $\zeta_{p,\eta}=\frac{\mathcal{D}}{\Lambda_{p,\eta}}$ is the acoustic thickness of the domain ($\mathcal{D}$) with respect to mean free path per polarization and spectral group. The value of $\zeta_{p,\eta}$ changes in each group and polarization; certain group and mode polarizations are acoustically thick and affect the convergence of the numerical simulation. Table 3 contains required iterations for each simulation and bounding values of $\zeta_{p,\eta}$ for the polarization and mode simulated. An increasing domain size requires more iterations to converge the solution, shown in Table 3. This is a consequence of the acoustic thickness; in the set $\left\{p,\Lambda_{\eta}\right\}$, the minimum value of $\Lambda_{\eta}$ is the limiting factor. Due to the coupled nature of the transport simulation, the convergence performance of ballistic phonons is hindered by the diffuse phonons; in a decoupled simulation ballistic phonons converge much faster than diffuse phonons. The acoustic thickness in each polarization and mode are of primary interest in the numerical transport simulation. Depending on mode and polarization, $\Lambda_{p,\eta}$ can vary greatly due to its associated $v_{p,\eta}$ and $\tau_{p,\eta}$, affecting the acoustic thickness and causing significant degradation of numerical convergence, shown in Table 3. As optic phonons do not contribute significantly to thermal conductivity in these simulations, elimination of optic modes would reduce the overall amount of iterations; the simulation converges only as fast as its’ slowest converging group. For the TA phonons, the implementation of an acceleration scheme such as diffusion synthetic acceleration (DSA) could reduce the required number of iterations.

One of the drawbacks of the coupled temperature method is the implicit dependency of all phonon channels; that is, if one or more of the channels is very diffuse, the simulations require many iterations to converge, if convergence is achieved at all. We have demonstrated convergence on domains ranging from 10 nm to 10 $\mu$m, but encountered significant decrease in efficiency when domain size increases to thousands of times the smallest value of $\Lambda_{\eta}$ (large $\zeta_{p,\eta}$). This is expected, as phonons in diffuse channels are highly scattering. We have observed similar, slowly convergent behavior in previous decoupled transport simulations where $\Lambda$ was small compared to the domain size [22, 23]. The distinction in a coupled implementation is that ballistic phonons are throttled by the presence of diffuse phonons, and convergence behavior is significantly affected. The involvement of a numerical acceleration scheme such as nonlinear diffusion acceleration (NDA) may improve convergence properties; with NDA, the solution of a lower order (with no angular variable) diffusion problem informs a source term in a higher order transport problem. This is a proven technique in neutron transport [17, 27, 28, 36], but may require significant revision for use in a phonon transport simulation and is planned in a future study.

Zhang et al. simulate a similar temperature coupling through informing each discrete phonon group by the overall temperature of the domain, and were able to simulate large spatial domains, up to 100 $\mu$m [16]. However, they neglected optic phonons in their coupled simulations, and so convergence properties of their approach may be artificially positive as convergence of phonon transport problems is limited by the smallest value of $\Lambda$, e.g., the most diffuse phonon groups. Multiple researchers [2, 37, 38, 39] have cited the necessity in simulating optic phonons due to the acoustic-optic coupling effect found in various materials – all phonons are coupled through the total energy of the domain. The inclusion of optic phonons effectively places a muzzle on the TA and LA phonons and prevents them from carrying their maximal potential of heat. While we study Si, as did Zhang et al., the development of this transport framework is desired to be general enough in its approach to be used in simulating materials which may have strong acousto-optic coupling effects.

We performed simulations isolating LA and TA phonons to gain insight into how much optic phonons contribute to $\kappa_{\textrm{eff}}$ depending on the size of the geometric domain. Table 2 shows the difference in $\kappa_{\textrm{eff}}$ with and without optic phonon contribution. There is a nontrivial difference in $\kappa_{\textrm{eff}}$ when domain size is small, but this effect diminishes as the simulation becomes more diffuse with increasing acoustic thickness. It is important to recognize the contribution of optic phonons in smaller spatial domains; at these sizes, optic phonons can exist “agnostically” between the ballistic and diffuse regimes. The presence of optic phonons causes significant degradation in the performance of the solver when comparing iterations between Table 3 and 4. While optic phonons may not be required in large simulation cells due to overshadowing influence of LA and TA phonons, it is nonetheless important to consider their contribution $\kappa_{\textrm{eff}}$, and subsequent effect on solver performance.

IV Conclusions

We have developed a temperature coupled, spectral method for simulating effective thermal conductivity from nano- to micro-scale using material properties for homogeneous silicon obtained from ab initio density functional theory. In this study we simulated 2D geometric domains to characterize algorithmic and transport performance. The extension of this method to simulate three-dimensional structures is trivial, but will demand additional spatial and angular degrees of freedom. Angular degrees of freedom scale as $N(N+2)$ for 3D and $N(N+2)/2$ for 2D, where $N$ is the chosen order of the discrete ordinates angular discretization. Spatial degrees of freedom scale with the nodal count in the selected finite element type, e.g., 3 nodes in 2D for a triangular element, 4 nodes in 3D for a tetragonal element. Computational cost will scale proportional to degrees of freedom in the equation system. We demonstrated correlation between the geometric dependency of heat flux and thermal conductivity, and our results for $\kappa_{\textrm{eff}}$ in Si are in good agreement with those available in the open literature. Convergence properties are negatively impacted by the presence of optic phonons; the diffuse characteristics of optic phonons increase the acoustic thickness of simulations. A closure term, $\beta$, was developed to conserve the local system energy, and this fraction is projected as a source term in each phonon group, providing additional systemic coupling through the energy bath. The range and sign of $\beta$ indicates the degree of disorder between the transport system and the total energy of the phonon bath. In small geometric domains, the importance of using BTE simulations to model heat transport cannot be overstated; the difference between $\Phi^{\textrm{T}}$ and $\Phi^{00}$ is large where equilibrium cannot be established, and spatial effects dramatically change temperature distributions and heat flux. Previous works have shown optic phonons to contribute little to overall heat flux in silicon [40, 41, 16, 42]. However, geometric domains less than $100$ nm in size (Table 2), and in devices where precision is critical, optic phonons should be included in simulation to accurately compute $\kappa_{\textrm{eff}}$. This work focused on purely homogeneous materials to characterize a new method, but our approach can be extended to bulk systems with voids (porous materials for thermoelectric devices) and heterogeneous systems (systems with inclusions, nuclear fuel with fission product defects) quite readily with the addition of an interface physics model. Subsequent improvement on this work would include simulation of transference between polarization branches, and the additional modeling of 3-phonon scattering through the relaxation time parameter included in the total relaxation time via Matthiessen’s rule. We are currently developing diffusion acceleration techniques adapted from neutron transport methods in order to speed up convergence properties as domains become extremely acoustically thick.

Detailed thermal transport requires development of high precision tools with which to perform modeling and simulation efforts. In turn, these efforts can guide the design of nano- and micro-scale devices, and the prediction of thermal behavior at the engineering scale. Perturbations at the atomic-molecular scale drive changes in the microstructure; ab initio driven BTE simulations are a promising development in a scale-bridging framework, representing one element in the larger problem of informing engineering scale simulations.

Acknowledgements.

The authors thank Sebastian Schunert and Yaqi Wang at INL for their invaluable knowledge and assistance. This research made use of the resources of the High-Performance Computing Center at Idaho National Laboratory, which is supported by the Office of Nuclear Energy of the U.S. Department of Energy and the Nuclear Science User Facilities under Contract No. DE-AC07-05ID14517. We also acknowledge the support of Pacific Northwest National Laboratory.