Unified theory of thermal transport in crystals and disordered solids

TL;DR

This paper develops a unified theoretical framework for thermal transport in both crystals and disordered solids, bridging the gap between existing models and accurately predicting thermal properties across different regimes.

Contribution

It introduces a comprehensive transport equation that encompasses both Peierls and Allen-Feldman limits, accounting for anharmonicity and disorder effects simultaneously.

Findings

The derived equation reduces to known limits in respective regimes.

Successfully predicts thermal conductivity of ultralow-conductivity thermoelectric materials.

Provides a unified understanding of heat conduction mechanisms in diverse solids.

Abstract

Crystals and glasses exhibit fundamentally different heat conduction mechanisms: the periodicity of crystals allows for the excitation of propagating vibrational waves that carry heat, as first discussed by Peierls; in glasses, the lack of periodicity breaks Peierls' picture and heat is mainly carried by the coupling of vibrational modes, often described by a harmonic theory introduced by Allen and Feldman. Anharmonicity or disorder are thus the limiting factors for thermal conductivity in crystals or glasses; hitherto, no transport equation has been able to account for both. Here, we derive such equation, resulting in a thermal conductivity that reduces to the Peierls and Allen-Feldman limits, respectively, in anharmonic-and-ordered or harmonic-and-disordered solids, while also covering the intermediate regimes where both effects are relevant. This approach also solves the long-standing problem of accurately predicting the thermal properties of crystals with ultralow or glass-like thermal conductivity, as we show with an application to a thermoelectric material representative of this class.

Tables

Table 1. Table 1: Lattice parameters of CsPbBr 3 in the orthorhombic phase, obtained by first-principles simulations or experiments lee2017ultralow ; Miyatae1701217 . Inset: unit cell of orthorhombic CsPbBr 3 visualized with VESTA momma2011vesta ; the directions x ^ , y ^ ^ 𝑥 ^ 𝑦 \hat{x},\hat{y} and z ^ ^ 𝑧 \hat{z} discussed e.g. in figure ( 1 ) refer respectively to lattice vectors 𝒂 1 subscript 𝒂 1 \bm{a}_{1} , 𝒂 2 subscript 𝒂 2 \bm{a}_{2} and 𝒂 3 subscript 𝒂 3 \bm{a}_{3} .

	$a_1$ [Å]	$a_2$ [Å]	$a_3$ [Å]
Present (th)	7.963	8.389	11.632
Ref.lee2017ultralow (th)	7.990	8.416	11.674
Ref.exp_cell_param (exp)	8.198	8.244	11.735
Ref.Miyatae1701217 (exp)	8.223	8.243	11.761

Equations

$$\mathrm{a}(\bm{q})_{b\alpha}{=}\tfrac{1}{\sqrt{2\hbar}}\Big{(}\sqrt[4]{\mathsfit{D}^{{-}1}(\bm{q})}_{b\alpha,b^{\prime}\alpha^{\prime}}{\mathrm{P}}(\bm{q})_{b^{\prime}\alpha^{\prime}}{-}i\sqrt[4]{\mathsfit{D}(\bm{q})}_{b\alpha,b^{\prime}\alpha^{\prime}}{\mathrm{U}}^{\dagger}(\bm{q})_{b^{\prime}\alpha^{\prime}}\Big{)}.$$

$$\begin{split}&\mathrm{H}^{\rm har}{=}\hbar\sum\limits_{\bm{q}}\sqrt{\mathsfit{D}(\bm{q})}_{b\alpha,b^{\prime}\alpha^{\prime}}\big{(}\mathrm{a}^{\dagger}(\bm{q})_{b\alpha}\mathrm{a}(\bm{q})_{b^{\prime}\alpha^{\prime}}{+}\frac{1}{2}\delta_{b,b^{\prime}}\delta_{\alpha,\alpha^{\prime}}\big{)}.\hskip 19.91692pt\end{split}$$

$$\frac{{\partial}{\rho}}{\partial t}(t)+\frac{i}{\hbar}\Big{[}\mathrm{H}^{\rm har},{\rho}(t)\Big{]}=\frac{\partial\rho(t)}{\partial t}\bigg{|}_{\mathrm{H}^{\rm col}};$$

$$\begin{split}&\frac{\partial}{\partial t}{\bm{\rho_{1}}}(\bm{q},\bm{q^{\prime}},t)+i\Big{(}\sqrt{\bm{\mathsfit{D}}(\bm{q})}\bm{\rho_{1}}(\bm{q},\bm{q^{\prime}},t)-\bm{\rho_{1}}(\bm{q},\bm{q^{\prime}},t)\sqrt{\bm{\mathsfit{D}}(\bm{q^{\prime}})}\Big{)}\\ &\hskip 142.26378pt=\frac{\partial}{\partial t}{\bm{\rho_{1}}({\bm{q}},{\bm{q}^{\prime}},t)}\Big{|}_{\mathrm{H}^{\rm col}}.\end{split}$$

$$\begin{split}{\mathsfit{W}}{(}{\bm{R}},{\bm{q}},{t}{)}_{{b}{\alpha}{,}{b^{\prime}}{\alpha^{\prime}}}{=}{\sum\limits_{\bm{q^{\prime\prime}}}}{\rho_{1}}{\big{(}}{\bm{q}}{+}{\bm{q^{\prime\prime}}},{\bm{q}}{-}{\bm{q^{\prime\prime}}},{t}{\big{)}}_{{b}{\alpha}{,}{b^{\prime}}{\alpha^{\prime}}}e^{2i\bm{q}^{\prime\prime}{\cdot}\bm{R}}.\hskip 19.91692pt\end{split}$$

$$\begin{split}&\frac{\partial}{\partial t}\bm{\mathsfit{W}}({\bm{R}},{\bm{q}},t)+i\Big{[}\sqrt{\bm{\mathsfit{D}}(\bm{q})},\bm{\mathsfit{W}}(\bm{R},\bm{q},t)\Big{]}+\hskip 56.9055pt\\ &+\frac{1}{2}\Big{\{}\vec{\nabla}_{\bm{q}}\sqrt{{\bm{D}}(\bm{q})},\cdot\vec{\nabla}_{\bm{R}}\bm{\mathsfit{W}}(\bm{R},\bm{q},t)\Big{\}}=\frac{\partial}{\partial t}\bm{\mathsfit{W}}(\bm{R},{\bm{q}},t)\Big{|}_{\mathrm{H}^{\rm col}};\end{split}$$

$$N(\mathbf{R},\mathbf{q},t{)}_{s,s^{\prime }}={\mathcal{E}}^{\ast }(\mathbf{q}{)}_{s,b\alpha }{W(\mathbf{R},\mathbf{q},t)}_{b\alpha ,b^{\prime }{\alpha }^{\prime }}\mathcal{E}(\mathbf{q}{)}_{s^{\prime },b^{\prime }{\alpha }^{\prime }},$$

$$V^{\beta }(\mathbf{q}{)}_{s,s^{\prime }}={\mathcal{E}}^{\ast }(\mathbf{q}{)}_{s,b\alpha }{{\nabla }_{\mathbf{q}}^{\beta }\sqrt{D(\mathbf{q})}}_{b\alpha ,b^{\prime }{\alpha }^{\prime }}\mathcal{E}(\mathbf{q}{)}_{s^{\prime },b^{\prime }{\alpha }^{\prime }}$$

$$\frac{\partial}{\partial t}\bm{\mathsfit{N}}({\bm{R}},{\bm{q}},t)+i\Big{[}\bm{\Omega}(\bm{q}),\bm{\mathsfit{N}}(\bm{R},\bm{q},t)\Big{]}+\frac{1}{2}\Big{\{}\vec{\bm{\mathsfit{V}}}(\bm{q}),\cdot\vec{\nabla}_{\bm{R}}\bm{\mathsfit{N}}(\bm{R},\bm{q},t)\Big{\}}=\frac{\partial}{\partial t}\bm{\mathsfit{N}}(\bm{R},{\bm{q}},t)\bigg{|}_{\mathrm{H}^{\rm col}}.$$

$$\bm{J}(\bm{R},t){=}\frac{1}{2}\frac{1}{\mathcal{V}N_{\rm c}}\sum_{\bm{q},s}{\Big{(}\big{\{}\vec{\bm{\mathsfit{V}}}(\bm{q}),\bm{\mathsfit{N}}(\bm{R},\bm{q},t)\big{\}}\hbar\bm{\Omega}(\bm{q})\Big{)}}_{s,s}.$$

$$\begin{split}\frac{\partial}{\partial t}&\mathsfit{N}({\bm{R},\bm{q},t})_{s{,}s^{\prime}}\Big{|}_{\mathrm{H}^{\rm col}}{=}\hskip 2.84526pt-(1-\delta_{s,s^{\prime}})\frac{\Gamma(\bm{q})_{s}+\Gamma(\bm{q})_{s^{\prime}}}{2}\mathsfit{N}(\bm{R},\bm{q},t)_{s,s^{\prime}}\\ &{-}\frac{\delta_{s,s^{\prime}}}{\mathcal{V}N_{\rm c}}{\sum\limits_{{s^{\prime\prime}}{\bm{q}^{\prime\prime}}}}\mathsfit{{A}}^{T}(\bm{q},{\bm{q}^{\prime\prime}})_{s,{s^{\prime\prime}}}\big{(}\mathsfit{N}(\bm{R},{\bm{q}^{\prime\prime}},t)_{{s^{\prime\prime}},{s^{\prime\prime}}}{-}\bar{\mathsfit{N}}^{T}({\bm{q}^{\prime\prime}})_{{s^{\prime\prime}}}\big{)};\end{split}$$

$$\begin{array}{rl}{\kappa }^{\alpha \beta }={\kappa }_{\mathrm{P}}^{\alpha \beta }+\frac{{\hslash }^2}{k_B{T}^2}\frac{1}{\mathcal{V}{N}_{\mathrm{c}}}& \sum _{\mathbf{q}}\sum _{s\ne s^{\prime }}\frac{\omega (\mathbf{q}{)}_s+\omega (\mathbf{q}{)}_{s^{\prime }}}{2}{V}^{\alpha }(\mathbf{q}{)}_{s,s^{\prime }}{V}^{\beta }(\mathbf{q}{)}_{s^{\prime },s}\times \\ & \times \frac{\omega (\mathbf{q}{)}_s{\bar{N}}^T(\mathbf{q}{)}_s[{\bar{N}}^T(\mathbf{q}{)}_s+1]+\omega (\mathbf{q}{)}_{s^{\prime }}{\bar{N}}^T(\mathbf{q}{)}_{s^{\prime }}[{\bar{N}}^T(\mathbf{q}{)}_{s^{\prime }}+1]}{4[\omega (\mathbf{q}{)}_{s^{\prime }}-\omega (\mathbf{q}{)}_s{]}^2+[\Gamma (\mathbf{q}{)}_s+\Gamma (\mathbf{q}{)}_{s^{\prime }}{]}^2}[\Gamma (\mathbf{q}{)}_s+\Gamma (\mathbf{q}{)}_{s^{\prime }}];\end{array}$$

$$\begin{split}\mathrm{a}(\bm{R})_{b\alpha}{=}&\frac{1}{\sqrt{N_{\rm c}}}\sum_{\bm{q}}\mathrm{a}(\bm{q})_{b\alpha}e^{+i\bm{q}\cdot\bm{R}}{=}\\ {=}&\frac{1}{\sqrt{2\hbar}}\sum_{\bm{R^{\prime\prime}}}\bigg{(}\sqrt[4]{\mathsfit{G}^{{-}1}}_{\bm{R}b\alpha,(\bm{R}{+}\bm{R^{\prime\prime}})b^{\prime}\alpha^{\prime}}\frac{{\mathrm{p}}(\bm{R}{+}\bm{R^{\prime\prime}})_{b^{\prime}\alpha^{\prime}}}{\sqrt{M_{b^{\prime}}}}\\ &{-}i\sqrt{M_{b^{\prime}}}{\mathrm{u}}(\bm{R}{+}\bm{R^{\prime\prime}})_{b^{\prime}\alpha^{\prime}}\sqrt[4]{\mathsfit{G}}_{(\bm{R}{+}\bm{R^{\prime\prime}})b^{\prime}\alpha^{\prime},\bm{R}b\alpha}\bigg{)};\end{split}$$

$$G_{\mathbf{R}b\alpha ,(\mathbf{R}+{\mathbf{R}}^{\prime \prime })b^{\prime }{\alpha }^{\prime }}=\frac{{\Phi }_{\mathbf{R}b\alpha ,(\mathbf{R}+{\mathbf{R}}^{\prime \prime })b^{\prime }{\alpha }^{\prime }}}{\sqrt{M_b{M}_{b^{\prime }}}};$$

$$\begin{split}&\mathcal{H}(\bm{R})_{b\alpha}{=}{\frac{\mathrm{p}^{2}(\bm{R})_{b\alpha}}{2{M_{b}}}}+\\ &+\frac{1}{2}\bigg{(}\sum_{\scriptscriptstyle{\bm{R^{\prime}}}}{\sqrt{\mathsfit{G}}_{\bm{R}b\alpha,(\bm{R}{+}\bm{R^{\prime}})b^{\prime}\alpha^{\prime}}}\sqrt{M_{b^{\prime}}}\mathrm{u}\left(\bm{R}{+}\bm{R^{\prime}}\right)_{b^{\prime}\alpha^{\prime}}\bigg{)}{\times}\\ &{\times}\bigg{(}\sum_{\scriptscriptstyle{\bm{R^{\prime\prime}}}}{\sqrt{\mathsfit{G}}_{\bm{R}b\alpha,(\bm{R}{+}\bm{R^{\prime\prime}})b^{\prime\prime}\alpha^{\prime\prime}}}\sqrt{M_{b^{\prime\prime}}}\mathrm{u}\left(\bm{R}{+}\bm{R^{\prime\prime}}\right)_{b^{\prime\prime}\alpha^{\prime\prime}}\bigg{)}.\end{split}$$

$$\begin{split}&\mathcal{H}(\bm{R})_{b\alpha}{=}\hbar\sum_{\bm{R^{\prime}},\bm{R^{\prime\prime}}}\sqrt[4]{\mathsfit{G}}_{\bm{R}b\alpha,(\bm{R}{+}\bm{R^{\prime}})b^{\prime}\alpha^{\prime}}\sqrt[4]{\mathsfit{G}}_{\bm{R}b\alpha,(\bm{R}{+}\bm{R^{\prime\prime}})b^{\prime\prime}\alpha^{\prime\prime}}{\times}\\ &{\times}\Big{(}\mathrm{a}^{\dagger}(\bm{R}{+}\bm{R^{\prime}})_{b^{\prime}\alpha^{\prime}}\mathrm{a}(\bm{R}{+}\bm{R^{\prime\prime}})_{b^{\prime\prime}\alpha^{\prime\prime}}+\frac{1}{2}\delta_{\bm{R^{\prime}},\bm{R^{\prime\prime}}}\delta_{b^{\prime},b^{\prime\prime}}\delta_{\alpha^{\prime},\alpha^{\prime\prime}}\Big{)}.\end{split}$$

$$\begin{split}E(\bm{R}{+}\bm{R^{\prime\prime}})_{{\bm{R}b\alpha}}{=}&\sum_{{b^{\prime\prime}\alpha^{\prime\prime}}}\big{<}0\big{|}\mathrm{a}(\bm{R})_{b\alpha}\mathcal{H}(\bm{R}{+}\bm{R^{\prime\prime}})_{b^{\prime\prime}\alpha^{\prime\prime}}\mathrm{a}^{\dagger}(\bm{R})_{b\alpha}\big{|}0\big{>}\\ &-\sum_{{b^{\prime\prime}\alpha^{\prime\prime}}}\big{<}0\big{|}\mathcal{H}(\bm{R}{+}\bm{R^{\prime\prime}})_{b^{\prime\prime}\alpha^{\prime\prime}}\big{|}0\big{>}\\ =&\hbar\sum_{{b^{\prime\prime}\alpha^{\prime\prime}}}\Big{(}\sqrt[4]{\mathsfit{G}}_{\bm{R^{\prime\prime}}b^{\prime\prime}\alpha^{\prime\prime},{\bm{0}}{b}{\alpha}}\Big{)}^{2},\end{split}$$

$$\mathrm{a}(\mathbf{q}{)}_{s^{\prime }}$$

$$\tilde{\mathrm{a}}(\mathbf{q}{)}_s=\mathrm{a}(\mathbf{q}{)}_s{e}^{i\mathbf{q}\cdot {\mathbf{R}}_{\Delta }};$$

$$\begin{array}{rl}& \mathrm{a}(\mathbf{R}{)}_s=\frac{1}{\sqrt{N_{\mathrm{c}}}}\sum _{\mathbf{q}}\mathrm{a}(\mathbf{q}{)}_s{e}^{+i\mathbf{q}\cdot \mathbf{R}};\\ & \tilde{\mathrm{a}}(\mathbf{R}{)}_s=\frac{1}{\sqrt{N_{\mathrm{c}}}}\sum _{\mathbf{q}}\tilde{\mathrm{a}}(\mathbf{q}{)}_s{e}^{+i\mathbf{q}\cdot \mathbf{R}}=\mathrm{a}(\mathbf{R}+{\mathbf{R}}_{\Delta }{)}_s.\end{array}$$

$$N(\mathbf{x},\mathbf{q}{)}_{s,s^{\prime }}={\bar{N}}^T(\mathbf{q}{)}_s{\delta }_{s,s^{\prime }}+{\bar{n}}^T(\mathbf{x},\mathbf{q}{)}_s{\delta }_{s,s^{\prime }}+{\vec{n}}^{(1)}(\mathbf{q}{)}_{s,s^{\prime }}\cdot \nabla T(\mathbf{x});$$

$$\begin{split}\vec{\mathsfit{V}}&(\bm{q})_{s,s}\cdot\left(\frac{\partial\mathsfit{N}^{T}(\bm{q})_{s}}{\partial T}\nabla T(\bm{x})\right)=\\ &=-\frac{1}{\mathcal{V}N_{\rm c}}\sum_{{s^{\prime\prime}}{\bm{q}^{\prime\prime}}}\mathsfit{A}^{T}(\bm{q},{\bm{q}^{\prime\prime}})_{s,{s^{\prime\prime}}}\big{(}\vec{\mathsfit{n}}^{(1)}({{\bm{q}^{\prime\prime}}})_{{s^{\prime\prime}},{s^{\prime\prime}}}{\cdot}\nabla T(\bm{x})\big{)};\end{split}$$

$$\begin{split}&\frac{1}{2}\left(\frac{\partial\mathsfit{N}^{T}(\bm{q})_{s}}{\partial T}+\frac{\partial\mathsfit{N}^{T}(\bm{q})_{s^{\prime}}}{\partial T}\right)\vec{\mathsfit{V}}(\bm{q})_{s,s^{\prime}}\cdot\nabla T(\bm{x})=\\ &-\left(i\big{(}\omega(\bm{q})_{s}-\omega(\bm{q})_{s^{\prime}}\big{)}+\frac{\Gamma(\bm{q})_{s}+\Gamma(\bm{q})_{s^{\prime}}}{2}\right)\vec{\mathsfit{n}}^{(1)}({\bm{q}})_{s,s^{\prime}}{\cdot}\nabla T(\bm{x}).\end{split}$$

$$\begin{array}{rl}& n^{(1),\alpha }(\mathbf{q}{)}_{s,s^{\prime }}=-\frac{\hslash }{k_B{T}^2}{V}^{\alpha }(\mathbf{q}{)}_{s,s^{\prime }}\times \\ & \frac{\omega (\mathbf{q}{)}_s{\bar{N}}^T(\mathbf{q}{)}_s[{\bar{N}}^T(\mathbf{q}{)}_s+1]+\omega (\mathbf{q}{)}_{s^{\prime }}{\bar{N}}^T(\mathbf{q}{)}_{s^{\prime }}[{\bar{N}}^T(\mathbf{q}{)}_{s^{\prime }}+1]}{2i[\omega (\mathbf{q}{)}_s-\omega (\mathbf{q}{)}_{s^{\prime }}]+[\Gamma (\mathbf{q}{)}_s+\Gamma (\mathbf{q}{)}_{s^{\prime }}]}.\end{array}$$

$${\kappa }_{{}_{\mathrm{P},\mathrm{SMA}}}^{\alpha \beta }=\frac{{\hslash }^2}{k_B{T}^2}\frac{1}{\mathcal{V}{N}_{\mathrm{c}}}\sum _{\mathbf{q},s}\frac{{\omega }^2(\mathbf{q}{)}_s{\bar{N}}^T(\mathbf{q}{)}_s[{\bar{N}}^T(\mathbf{q}{)}_s+1]}{\Gamma (\mathbf{q}{)}_s}{V^{\alpha }(\mathbf{q})}_{s,s}{V}^{\beta }(\mathbf{q}{)}_{s,s}.$$

$$\begin{array}{rl}{\kappa }_{{}_{\mathrm{SMA}}}^{\alpha \beta }=\frac{{\hslash }^2}{k_B{T}^2}\frac{1}{\mathcal{V}{N}_{\mathrm{c}}}\sum _{\mathbf{q}}\sum _{s,s^{\prime }}& \frac{\omega (\mathbf{q}{)}_s+\omega (\mathbf{q}{)}_{s^{\prime }}}{2}{V}^{\alpha }(\mathbf{q}{)}_{s,s^{\prime }}{V}^{\beta }(\mathbf{q}{)}_{s^{\prime },s}\times \\ & \times \frac{\omega (\mathbf{q}{)}_s{\bar{N}}^T(\mathbf{q}{)}_s[{\bar{N}}^T(\mathbf{q}{)}_s+1]+\omega (\mathbf{q}{)}_{s^{\prime }}{\bar{N}}^T(\mathbf{q}{)}_{s^{\prime }}[{\bar{N}}^T(\mathbf{q}{)}_{s^{\prime }}+1]}{4[\omega (\mathbf{q}{)}_{s^{\prime }}-\omega (\mathbf{q}{)}_s{]}^2+[\Gamma (\mathbf{q}{)}_s+\Gamma (\mathbf{q}{)}_{s^{\prime }}{]}^2}[\Gamma (\mathbf{q}{)}_s+\Gamma (\mathbf{q}{)}_{s^{\prime }}];\end{array}$$

$${\kappa }_A=\frac{1}{V_{\mathrm{t}}}{C}_i(T)D_i(T);$$

$$C_i(T)=\left(\frac{\hslash {\omega }_i^2}{V_{\mathrm{t}}T}\right)\left(-\frac{\partial {\bar{N}}_i}{\partial {\omega }_i}\right);$$

$$D_i(T)=\frac{\pi V_{\mathrm{t}}^2}{3{\hslash }^2{\omega }_i^2}\frac{{\hslash }^2}{V^2}\sum _{j\ne i}{\left|{\vec{V}}_{\text{AA}}(\mathbf{q}{)}_{i,j}\frac{{\omega }_i+{\omega }_j}{2}\right|}^2\delta ({\omega }_i-{\omega }_j).$$

$$V_{\text{AA}}^{\beta }(\mathbf{q}{)}_{s,s^{\prime }}=\frac{\omega (\mathbf{q}{)}_s+\omega (\mathbf{q}{)}_{s^{\prime }}}{2\sqrt{\omega (\mathbf{q}{)}_s\omega (\mathbf{q}{)}_{s^{\prime }}}}{V}^{\beta }(\mathbf{q}{)}_{s,s^{\prime }}.$$

$$-\frac{\partial \bar{N}(\mathbf{q}{)}_s}{\partial \omega (\mathbf{q}{)}_s}\frac{\hslash \omega (\mathbf{q}{)}_s}{V_{\mathrm{t}}T}=\frac{1}{\mathcal{V}{N}_{\mathrm{c}}}\frac{{\hslash }^2}{k_B{T}^2}\omega (\mathbf{q}{)}_s\bar{N}(\mathbf{q}{)}_s(\bar{N}(\mathbf{q}{)}_s+1);$$

Full text

Unified theory of thermal transport in        

crystals and disordered solids

Michele Simoncelli1

  

Nicola Marzari1

  

Francesco Mauri2⋆

**Crystals and glasses exhibit fundamentally different heat conduction mechanisms: the periodicity of crystals allows for the excitation of propagating vibrational waves that carry heat, as first discussed by Peierlspeierls1929kinetischen; in glasses, the lack of periodicity breaks Peierls’ picture and heat is mainly carried by the coupling of vibrational modes, often described by a harmonic theory introduced by Allen and Feldmanallen1989thermal. Anharmonicity or disorder are thus the limiting factors for thermal conductivity in crystals or glasses; hitherto, no transport equation has been able to account for both. Here, we derive such equation, resulting in a thermal conductivity that reduces to the Peierls and Allen-Feldman limits, respectively, in anharmonic-and-ordered or harmonic-and-disordered solids, while also covering the intermediate regimes where both effects are relevant. This approach also solves the long-standing problem of accurately predicting the thermal properties of crystals with ultralow or glass-like thermal conductivitysun2010lattice; li2015ultralow; wang2018cation; lee2017ultralow; chen2015twisting; voneshen2013suppression; lory2017direct; PhysRevB.96.214202, as we show with an application to a thermoelectric material representative of this class. **

In 1929 Peierlspeierls1929kinetischen developed a semi-classical theory for heat conduction in terms of a Boltzmann transport equation (BTE) for propagating phonon wave packets. Nowadays, modern algorithms and computing systems allow to solve its linearized form (LBTE) either approximately (in the so-called single mode approximation (SMA)PhysRevLett.106.045901) or exactly, using iterativeomini1995iterative; carrete2017almabte, variationalfugallo2013ab, or exact diagonalizationPhysRevLett.110.265506; cepellotti2016thermal methods; its accuracy has been highlighted in many studiesfugallo2013ab; cepellotti2016thermal; cepellotti2015phonon; PhysRevLett.110.265506. Nevertheless, these cases are characterized by having few, well-separated phonon branches and anharmonicity-limited thermal conductivity; we will refer to these in the following as “simple” crystals. In 1963 Hardy was able to express the thermal conductivity in terms of the phonon velocity operatorhardy1963energy and showed that its diagonal elements are the phonon group velocities entering the Peierls’ BTE, while the off-diagonal terms, missing from it, are actually negligible in simple crystalshardy1963energy. In 1989 Allen and Feldmanallen1989thermal envisioned that these off-diagonal elements, neglected so far, could become dominant in disordered regimes, where Peierls’ picture breaks down due to the impossibility of defining phonons and group velocities. As a consequence, a harmonic theory of thermal transport in glasses was introduced, where disorder limits thermal conductivity and heat is carried by couplings of vibrational modes arising from the off-diagonal elements of the velocity operator (diffusons and loconsallen1999diffusons; lv2016non). Recently, it has been argued that the diffuson conduction mechanism, typical of glasses, can also emerge in a third class of materials, termed “complex” crystalschen2015twisting, characterized by large unit cells and many quasi-degenerate phonon branches, where it coexists with phonon transport. Conversely, crystal-like propagation mechanisms have been suggested also for glasses (propagonsallen1999diffusons) — albeit without a formal justification — in order to explain the experimental resultsallen1999diffusons.

Here, we show that a more general theory exists, encompassing the emergence and coexistence of all these excitations, and reducing to the BTE in the semi-classical limit of a simple crystal, or to the Allen-Feldman formulation in the case of a harmonic glass.

We start by considering a crystal in the Born-Oppenheimer approximation, where we add perturbatively to a harmonic Hamiltonianziman1960electrons $\mathrm{H}^{\rm har}$ a collision term $\mathrm{H}^{\rm col}$ accounting for anharmonicityziman1960electrons and isotopic disordertamura1983isotope. Heat transfer is driven by a temperature gradient in real space (i.e. a space-dependent lattice vibrational energy) and in order to derive a more general transport equation we want to track excitations both in real and Fourier space. With this goal (see Methods for details) we introduce the bosonic operator $\mathrm{a}(\bm{q})_{b\alpha}$ that annihilates vibrations with wave-vector $\bm{q}$ for all atoms $b$ in direction $\alpha$ (the Einstein summation convention is implied over repeated indexes)

[TABLE]

Here ${\mathrm{U}}(\bm{q})_{b\alpha}{=}\sqrt{{{M_{b}}}/{{N_{\rm c}}}}\sum_{\bm{R}}\mathrm{u}(\bm{R})_{b\alpha}e^{{+}i\bm{q}\cdot\bm{R}}$ and ${\mathrm{P}}(\bm{q})_{b\alpha}{=}({\sqrt{M_{b}N_{\rm c}}})^{-1}\sum_{\bm{R}}\mathrm{p}(\bm{R})_{b\alpha}e^{-i\bm{q}\cdot\bm{R}}$ are the Fourier representations of the displacement and momentum operators $\mathrm{u}(\bm{R})_{b\alpha}$ and $\mathrm{p}(\bm{R})_{b\alpha}$ of atom $b$ in unit cell $\bm{R}$, with mass $M_{b}$ and in a lattice containing $N_{\rm c}$ unit cells, and $\bm{\mathsfit{D}}(\bm{q})$ is the dynamical matrix of the crystal with eigenvalues $\omega^{2}(\bm{q})_{s}$ ($s$ labels the $3{\cdot}N_{\rm at}$ phonon branches, $N_{\rm at}$ being the number of atoms in the unit cell). The usual commutation relations $[\mathrm{a}(\bm{q})_{b\alpha},\mathrm{a}^{\dagger}(\bm{q}^{\prime})_{b^{\prime}\alpha^{\prime}}]{=}\delta_{\bm{q,q^{\prime}}}\delta_{b,b^{\prime}}\delta_{\alpha,\alpha^{\prime}}$ are satisfied, and in terms of these operators $\mathrm{H}^{\rm har}$ assumes the simple form:

[TABLE]

When a temperature gradient is imposed with couplings to external baths, the system’s state — described by the density matrix $\rho$ — undergoes an irreversible Markovian evolution described by the master equationgebauer2004kinetic

[TABLE]

where the collision term \frac{\partial\rho(t)}{\partial t}\big{|}_{\mathrm{H}^{\rm col}} accounts for the transitions between the eigenstates of $\mathrm{H}^{\rm har}$ driven by the perturbation $\mathrm{H}^{\rm col}$. Here we are interested in the vibrational energy, which is a one-body operator at the harmonic leading order. We can thus move to a formulation based on the one-body density matrix $\bm{\rho_{1}}$, which is defined in terms of a trace over the Fock spacegebauer2004kinetic: \rho_{1}(\bm{q},\bm{q^{\prime}},t)_{b\alpha,b^{\prime}\alpha^{\prime}}{=}{\rm Tr}\big{(}{\rho}(t)\;{a}^{\dagger}(\bm{q^{\prime}})_{b^{\prime}\alpha^{\prime}}\;{a}(\bm{q})_{b\alpha}\big{)}. It follows from equation (3) that $\bm{\rho_{1}}$ evolves according to:

[TABLE]

For systems that are homogeneous in real space $\bm{\rho_{1}}$ is diagonal: $\bm{\rho_{1}}(\bm{q},\bm{q^{\prime}},t){=}\bm{\rho_{1}}(\bm{q},\bm{q^{\prime}},t)\delta_{\bm{q},\bm{q^{\prime}}}$; to study non-homogeneous systems it becomes useful to introduce the Wigner distribution $\bm{\mathsfit{W}}$, obtained by applying the Wigner-Weyl transformation to $\bm{\rho_{1}}$frensley1990boundary:

[TABLE]

It becomes apparent from equation (5) that a non-homogeneity in real space is related to a one-body density matrix with off-diagonal elements in reciprocal space. Now, we consider the out-of-equilibrium regime of a system perturbed by a small temperature gradient, appreciable over a length scale $L$ much larger than the norms of the direct lattice vectors $|\bm{a}_{i}|$ ($i{=}1,2,3$). Applying the transformation (5) to equation (4), we can obtain the evolution for the Wigner distribution. Such evolution equation can be greatly simplified by using the hypothesis of having non-homogeneities only at mesoscopic scales or above, i.e. of having a one-body density matrix \bm{\rho_{1}}{\big{(}}{\bm{q}}{+}{\bm{q^{\prime\prime}}},{\bm{q}}{-}{\bm{q^{\prime\prime}}},{t}{\big{)}} that is sharply peaked and significantly different from zero only for $|\bm{q^{\prime\prime}}|{\ll}{2\pi}|\bm{a}_{i}|^{-1}$. This allows us to perform a Taylor expansion for \sqrt{\bm{\mathsfit{D}}\big{(}{\bm{q}}{\pm}\bm{q^{\prime\prime}}\big{)}} around $\bm{q^{\prime\prime}}{\to}0$, obtaining the simplified equation:

[TABLE]

where the gradients appearing in the anticommutator are contracted using the scalar product. The spatial gradient $\vec{\nabla}_{\bm{R}}$ acts on the variable $\bm{R}$ appearing in the exponential $e^{2i\bm{q^{\prime\prime}}\cdot\bm{R}}$, whose characteristic length scale is $L{\sim}|\bm{q}^{\prime\prime}|^{-1}{\gg}|\bm{a}_{i}|$. Over such mesoscopic scale $L$ the discrete variable $\bm{R}$ can be considered as continuous and the gradient $\vec{\nabla}_{\bm{R}}$ well-defined. At this point, thanks to the locality in $\bm{q}$ of equation (6), we can recast it in the phonon representation by applying the unitary transformation that diagonalizes ${\bm{\mathsfit{D}}(\bm{q})}$: $\mathcal{E}^{*}(\bm{q})_{{s},b\alpha}{\mathsfit{D}(\bm{q})}_{b\alpha,b^{\prime}\alpha^{\prime}}\mathcal{E}(\bm{q})_{s^{\prime},b^{\prime}\alpha^{\prime}}{=}\omega^{2}(\bm{q})_{s}\delta_{s,s^{\prime}}$. We thus introduce the following quantities

[TABLE]

where the atomic and Cartesian indexes are replaced by the phonon branch index $s$. We define the diagonal matrix $\Omega(\bm{q})_{s,s^{\prime}}{=}\omega(\bm{q})_{s}\delta_{s,s^{\prime}}$ and denote with \frac{\partial\mathsfit{N}(\bm{R},{\bm{q}},t)_{s,s^{\prime}}}{\partial t}\big{|}_{\mathrm{H}^{\rm col}} the collision operator in the basis of the phonon modes. Equations (7) and (8) yield respectively a generalization of the phonon distributions and group velocities appearing in the BTEpeierls1929kinetischen; fugallo2013ab; in the absence of degeneracies (i.e. for $\omega(\bm{q})_{s}{\neq}\omega(\bm{q})_{s^{\prime}}\forall s{\neq}s^{\prime}$) the diagonal elements $s{=}s^{\prime}$ of equation (7) and of equation (8) coincide, respectively, with the phonon distributions appearing in Peierls’ equation and with their group velocities. Equation (6) in the phonon representation reads

[TABLE]

This equation is the first main result of this paper, and offers a more general phonon transport equation, where phonon wave packets are not only allowed to propagate particle-like in space but also to tunnel, wave-like, from one branch to another. We note that the left-hand side of equation (9) does not contain $\hbar$, and it does not distinguish between fermions or bosons. Quantum statistics is instead enforced by the collision term at the right-hand side (see e.g. equation (11) later). The vibrational energy density can now be determined from the solution of equation (9): E(\bm{R},t){=}\frac{\hbar}{\mathcal{V}N_{\rm c}}\sum_{\bm{q},s}\big{(}\bm{\mathsfit{N}}(\bm{R},\bm{q},t)\bm{\Omega}(\bm{q})\big{)}_{s,s}. Computing the time derivative of $E(\bm{R},t)$ and using equation (9), one obtains the harmonic heat flux $\bm{J}(\bm{R},t)$:

[TABLE]

To compute the thermal conductivity we focus on the steady-state regime and search for a solution of equation (9) linear in the temperature gradient. The right-hand side collision term can be written in a closed form approximation expressing the many-body collision operator (3) in terms of one-body quantities and linearizing it around the equilibrium temperature $T$frensley1990boundary; gebauer2004kinetic; cohen2004atom

[TABLE]

where $\mathcal{V}$ is the unit cell volume, \bar{\mathsfit{N}}^{T}({\bm{q}})_{{s}}{=}\big{(}{\exp\big{[}\frac{\hbar\omega(\bm{q})_{{s}}}{k_{B}T}\big{]}{-}1}\big{)}^{-1} is the equilibrium Bose-Einstein distribution, $\Gamma(\bm{q})_{s}{=}(\mathcal{V}N_{\rm c})^{-1}\mathsfit{A}^{T}(\bm{q},{\bm{q}})_{s,s^{\prime\prime}}\delta_{s^{\prime\prime},s}$ is the phonon linewidth and $\mathsfit{A}^{T}(\bm{q},{\bm{q}^{\prime\prime}})_{s,{s^{\prime\prime}}}$ is the standard collision operator that accounts for anharmonicity and isotopic disorder (for the third-order anharmonicities considered here, see equation (13) of Ref.fugallo2013ab). Under these conditions, equation (9) can be decoupled in its diagonal and off-diagonal parts, as shown in Methods; these two equations describe the evolution of the populations and coherences, respectively, using the standard nomenclature in this field. The equation for the populations turns out to be the standard Peierls’ LBTE that can be solved exactlyfugallo2013ab; omini1995iterative; carrete2017almabte; PhysRevLett.110.265506; cepellotti2016thermal; the equation for the coherences can be solved straightforwardly (see Methods). Once the solution $\mathsfit{N}(\bm{R},\bm{q})_{s,s^{\prime}}$ for both populations ($s{=}s^{\prime}$) and coherences ($s{\neq}s^{\prime}$) is known, one can compute the heat flux ${J}^{\alpha}(\bm{R})$ and the thermal conductivity $\kappa^{\alpha\beta}$ from ${J}^{\alpha}(\bm{R}){=}{-}\kappa^{\alpha\beta}\nabla T^{\beta}(\bm{R})$. The second main result of this work follows, where the thermal conductivity is generalized to an expression including both the populations’ and coherences’ contributions:

[TABLE]

where $\kappa^{\alpha\beta}_{\rm P}$ is the standard Peierls contribution to the conductivity and the additional tensor, which we name $\kappa^{\alpha\beta}_{\rm C}$, derives from the coherences’ equation. The relevance of equation (12) is its capability to account for heat transfer associated to both the diagonal (populations) and off-diagonal (coherences) Wigner distribution elements. The former is associated to the particle-like propagation of phonon wave packets discussed by Peierls’ semiclassical picture; the latter is related to the wave-like tunneling and loss of coherence between different vibrational eigenstates. Equation (12) highlights the conditions under which $\kappa^{\alpha\beta}_{\rm C}$ is relevant, namely when the difference between the phonon frequencies is of the same order of the average phonon linewidth and, at the same time, the velocity operator has non-negligible off-diagonal elements — we will show later that this is what precisely happens in complex crystals. In simple crystals, instead, the present approach and the standard LBTE provide very similar results since, as we will show later, $\kappa^{\alpha\beta}_{\rm P}{\gg}\kappa^{\alpha\beta}_{\rm C}$. Since the thermal conductivity (12) derives from an exact linear-order solution of equation (9), it remains valid also in the hydrodynamic regimecepellotti2015phonon. In the kinetic regimecepellotti2015phonon instead, a good estimation of $\kappa^{\alpha\beta}_{\rm P}$ can be obtained using the SMA approximation, i.e. neglecting the scattering events that repopulate the diagonal $\mathsfit{N}(\bm{R},{\bm{q}},t)_{{s},{s}}$fugallo2013ab. In practice, the SMA consists in replacing the collision matrix for the populations of equation (11) with its diagonal elements only: $\mathsfit{A}^{T}(\bm{q},{\bm{q^{\prime\prime}}})_{s,s^{\prime\prime}}{\simeq}(\mathcal{V}N_{\rm c})\Gamma(\bm{q})_{s}\delta_{\bm{q},\bm{q^{\prime\prime}}}\delta_{s,s^{\prime\prime}}$. The equation for the coherences is unaffected by this approximationcohen2004atom; it follows that the total SMA conductivity is given by an expression equal to the last term of equation (12), but with the summation running in full over $s$ and $s^{\prime}$ (equation (26) in Methods). Last, we show in Methods how in the limit of a harmonic crystal with an infinitely large unit cell (i.e. with the Brillouin zone reduced to the $\Gamma$-point) the normalized trace of the coherences’ conductivity tensor $\frac{1}{3}\kappa^{\alpha\alpha}_{\rm C}$ reduces to the Allen-Feldman expression for glassesallen1989thermal.

The perovskite CsPbBr3, which belongs to a family of materials promising for thermoelectric energy conversionwang2018cation; lee2017ultralow, is a good case study for this more general framework, as it features many quasi-degenerate phonon branches and its thermal propertieswang2018cation; lee2017ultralow are not correctly described by the LBTElee2017ultralow. We compute from first-principles the harmonic Hamiltonian and the collision operator (see details in Methods), calculating the complete thermal conductivity of equation (12) as a function of temperature.

The results are shown in figure (1) (we show here the SMA for $\kappa^{xx}$; this is almost identical to the full solution (see Methods)): at 300 K the populations term contributes just $30\%$ of the total conductivity, while the coherences term provides an additional $70\%$, leading to much closer agreement with experimentswang2018cation that becomes even more relevant in the high-temperature asymptotics. Conversely, at low temperature Peierls’ term becomes dominant (at 50 K it provides 78$\%$ of the total conductivity). It can be seen that in the high temperature limit $\kappa^{\alpha\beta}_{\rm P}{\propto}\;T^{-1}$, as predicted by Peierls’ theorysun2010lattice; lory2017direct; li2015ultralow; ziman1960electrons; this in broad disagreement with experiments. Instead, the present theory predicts a decay of $\kappa^{xx}$ much milder than $T^{-1}$ as shown here for CsPbBr3, and as also present in many other complex crystalsli2015ultralow; lory2017direct; PhysRevB.46.6131; PhysRevB.96.214202; chen2015twisting; voneshen2013suppression. We note in passing that the $\kappa^{\alpha\beta}_{\rm P}$ calculated here and the one presented in Ref.lee2017ultralow are very similar, confirming that all computational parameters are correctly accounted for. Figure (2a-c) shows that phonon branches with large group velocities and low anharmonicity provide a large contribution to the populations term; conversely, highly anharmonic flat bands contribute to the coherences term. The populations contribution originates mainly from low-frequency modes, while all the spectrum contributes to the coherences term. At 50 K, the 2-dimensional density of states for the coherences’ thermal conductivity ($\kappa^{xx}_{\rm{C,DOS}(2)}$, figure (2d)) shows couplings between quasi-degenerate states similarly to the case of harmonic glasses ($\Gamma(\bm{q})_{s}{\to}0$). At 300 K (figure (2e)) $\kappa^{xx}_{\rm{C,DOS}(2)}$ instead includes couplings between phonon modes having very different frequencies, driven by the large anharmonicity — therefore the corresponding heat conduction mechanism is intrinsically different from the one of harmonic glasses. In order to verify that the present theory reduces to the Peierls regime in the limit of simple crystals, we computed the room-temperature (300 K) thermal conductivity (12) of silicon and diamond. We found that for these crystals the coherences term is negligible: for silicon $\kappa_{\rm P}^{xx}{=}142.83\;\rm{W}/(\rm{m}{\cdot}\rm{K})$ vs. $\kappa_{\rm C}^{xx}{=}0.23\;\rm{W}/(\rm{m}{\cdot}\rm{K})$; for diamond $\kappa_{\rm P}^{xx}{=}2780.07\;\rm{W}/(\rm{m}{\cdot}\rm{K})$ vs. $\kappa_{\rm C}^{xx}{=}0.03\;\rm{W}/(\rm{m}{\cdot}\rm{K})$.

In summary, we have derived from the Wigner phase space formulation of quantum mechanics a general transport equation (9) able to account rigorously and simultaneously for anharmonicity and disorder, leading to a thermal conductivity (12) that adds an interbranch tunneling term for the coherences; this equation (12) reduces to the Peierls’ BTE in the limit of simple crystals, and it provides a rigorous derivation of the BTE for phonon populations first heuristically derived by Peierlspeierls1929kinetischen. Also, in the limit of harmonic glasses, the thermal conductivity (12) reduces to the Allen-Feldman limit; equation (12) is more general and encompasses all intermediate cases, applying also to anharmonic glasses. The LBTE approach, which accounts for populations only, greatly underestimates the thermal conductivity of complex crystalsli2015ultralow; lee2017ultralow; chen2015twisting; PhysRevB.96.214202; lory2017direct that ubiquitously display a glass-like high-temperature behavior (i.e. decay milder than the universal $T^{-1}$ trend predicted by Peierls’ theorysun2010lattice; lory2017direct; li2015ultralow; ziman1960electrons) and that up to now has been only heuristically explained by several modelsPhysRevB.46.6131; agne2018minimum; Clarke2003; chen2015twisting. In the present approach, the correct conductivity is recovered, and its high-temperature decay (figure (1)).

We note in passing that a complementary approach to these transport equations is to compute the thermal conductivity via the Green-Kubo methodmarcolongo2016microscopic; seyf2017rethinking; lv2016non; carbogno2017ab; puligheddu2017first, which is based on molecular dynamics and is exact, although computationally more expensive, provided one is in the classical Maxwell-Boltzmann high-temperature limit (i.e. at $T{\gg}T_{D}$, where $T_{D}$ is the Debye temperature). Coherences’ effects have been thoroughly studied in the context of electronic transportgebauer2004kinetic; PhysRevB.86.155433; frensley1990boundary and it is worth mentioning that the term $\kappa^{\alpha\beta}_{\rm C}$ has some analogies with the electronic current related to the Zener interband transitionsPhysRevB.86.155433. However, thermal transport is fundamentally different from electronic transport since the full phonon spectrum contributes to the heat flow, while for electrons the current originates only from the bands close to the Fermi level. This novel formulation is also going to be relevant to applied research, as it allows to predict the low thermal conductivity of target materials for thermoelectric energy conversionlee2017ultralow; voneshen2013suppression.

**Acknowledgements

**

**Author contributions

**The project was conceived by all authors. M.S. performed the numerical calculations and prepared the figures with inputs from N.M and F.M. All authors contributed to the redaction of the manuscript.

**Competing interests

**The authors declare no competing interests.

**Additional information

Correspondence and request of materials** should be addressed to F.M.

**Methods

Bosonic operators in the Cartesian basis.** In order to describe heat transfer driven by a temperature gradient (i.e. a spatially-dependent lattice vibrational energy), we need to keep track of the vibrations’ location in real space. This is possible using the Fourier transform of the bosonic operator defined in equation (1):

[TABLE]

since, as we will show now, it is related to vibrations localized around the lattice site indexed by $(\bm{R},b,\alpha)$. The tensor ${\mathsfit{G}}_{\bm{R}b\alpha,{(\bm{R}{+}\bm{R^{\prime\prime}})}b^{\prime}\alpha^{\prime}}$ appearing in equation (13) is related to the interatomic harmonic force constants in real spaceziman1960electrons, $\Phi_{\bm{R}b\alpha,(\bm{R}{+}\bm{R^{\prime\prime}})b^{\prime}\alpha^{\prime}}$:

[TABLE]

it is translation-invariant and decays to zero for $\bm{R^{\prime\prime}}{\to}\infty$ziman1960electrons; Simoncelli_full_2018. The harmonic Hamiltonian (2) can be written in terms of the vibrational energy field $\mathrm{H}^{\rm har}{=}\sum\limits_{\bm{R},b,\alpha}\mathcal{H}(\bm{R})_{b,\alpha}$, where

[TABLE]

In terms of the bosonic operators $\mathrm{a}^{\dagger}(\bm{R})_{b\alpha},\;\mathrm{a}(\bm{R})_{b\alpha}$, the energy field (15) becomes:

[TABLE]

From the decay to zero of ${\mathsfit{G}}_{\bm{R}b\alpha,{(\bm{R}{+}\bm{R^{\prime\prime}})}b^{\prime}\alpha^{\prime}}$ for $\bm{R^{\prime\prime}}{\to}\infty$, it follows that the excitation created by $\mathrm{a}^{\dagger}(\bm{R})_{b\alpha}$ has an energy distribution $E(\bm{R}{+}\bm{R^{\prime\prime}})_{{\bm{R}b\alpha}}$ localized around the position $\bm{R}$Simoncelli_full_2018:

[TABLE]

where we have used the translational invariance of the tensor $\bm{\mathsfit{G}}$ to simplify the last equation: $\sqrt[4]{\mathsfit{G}}_{(\bm{R}{+}\bm{R^{\prime\prime}})b^{\prime\prime}\alpha^{\prime\prime},{\bm{R}}{b}{\alpha}}{=}\sqrt[4]{\mathsfit{G}}_{\bm{R^{\prime\prime}}b^{\prime\prime}\alpha^{\prime\prime},{\bm{0}}{b}{\alpha}}$. We have therefore demonstrated that the operator (13) creates vibrations localized around the lattice position $\bm{R}$; hence, performing the Fourier transform of the operator (1) it is possible to keep track of the vibrations’ locations in real space. In contrast, the more common phonon-eigenmodes basis, obtained by the transformation (discussed in equations (7,8) of the main text):

[TABLE]

does not allow such a real space-based description. In fact, the operator $\mathrm{a}(\bm{q})_{s^{\prime}}$ does not have a direct correspondence with the location of vibrations in real space, due to the phase indeterminacy of the eigenstates $\mathcal{E}(\bm{q})_{{s^{\prime}},b^{\prime}\alpha^{\prime}}$ of the dynamical matrix $\bm{\mathsfit{D}}(\bm{q})$. This becomes apparent considering the following gauge transformation, allowed by the phase freedom on $\mathcal{E}(\bm{q})_{{s^{\prime}},b^{\prime}\alpha^{\prime}}$ at each $\bm{q}$-point:

[TABLE]

where $\bm{R}_{\Delta}$ is a direct lattice vector. Despite $\tilde{\mathrm{a}}(\bm{q})_{s}$ and $\mathrm{a}(\bm{q})_{s}$ being equivalent for a description in Fourier space, their real space representations are shifted by $\bm{R}_{\Delta}$:

[TABLE]

It is worth mentioning that this mirrors the electronic case, where the phase indeterminacy of the Bloch orbitals is reflected in the non-uniqueness of the transformation into Wannier functionsRevModPhys.84.1419. In this work, the use of the Cartesian basis — which has no phase freedom — renders the real space description and hence the Wigner distribution (5) well-defined.

Decoupling of equation (9). We consider the steady-state regime, in which a small, slowly varying temperature gradient $\nabla T$ is established throughout the system. We look for a solution linear in the temperature gradient $\nabla T$fugallo2013ab on the distribution corresponding to the local equilibrium temperaturehardy1970phonon $T(\bm{x})$:

[TABLE]

where $\bar{\mathsfit{n}}^{T}({\bm{x}},{\bm{q}})_{s}\delta_{s,s^{\prime}}$ accounts for the local equilibrium temperature: $\bar{\mathsfit{n}}^{T}({\bm{x}},{\bm{q}})_{s}{=}\frac{d\bar{\mathsfit{N}}^{T}({\bm{q}})_{s}}{dT}(T(\bm{x}){-}T)$; the term $\vec{\mathsfit{n}}^{(1)}({\bm{q}})_{s,s^{\prime}}{\cdot}\nabla T(\bm{x})$ is, in general, non-diagonal in $s,s^{\prime}$ and contains the information concerning the deviation of the full solution from the local equilibrium solution — it is assumed to be of the order of the temperature gradient, as in previous workfugallo2013ab; cepellotti2016thermal. Considering only the first-order term in $\nabla T$, we obtain an equation that is decoupled in its diagonal and off-diagonal parts. The diagonal part turns out to be the usual LBTE for the populationsfugallo2013ab:

[TABLE]

which can be solved approximately (in the SMAPhysRevLett.106.045901) or exactlyfugallo2013ab; cepellotti2016thermal; omini1995iterative; carrete2017almabte; PhysRevLett.110.265506. The equation for the coherences ($s{\neq}s^{\prime}$) is:

[TABLE]

Equation (23) can be solved exactly and straightforwardly, yielding ($s{\neq}s^{\prime}$):

[TABLE]

We stress that the simplicity of the equation for the coherences (23) implies that the computational cost for finding its solution is negligible compared to the one for solving exactly the populations’ equation (22).

**SMA approximation and properties of $\kappa^{\alpha\beta}$

**In the kinetic regimecepellotti2015phonon, a good estimate of the Peierls’ conductivity $\kappa^{\alpha\beta}_{\rm P}$ appearing in equation (12) is given by the SMA approximation, i.e. $\kappa^{\alpha\beta}_{\rm P}{\simeq}\kappa^{\alpha\beta}_{{}_{\rm P,SMA}}$ where

[TABLE]

Within the SMA approximation, the total thermal conductivity $\kappa^{\alpha\beta}$ is approximated by $\kappa^{\alpha\beta}_{{}_{\rm SMA}}{=}\kappa^{\alpha\beta}_{{}_{\rm P,SMA}}{+}\kappa^{\alpha\beta}_{{}_{\rm C}}$, which can be written in the following compact form:

[TABLE]

where we stress that the summation over $s,s^{\prime}$ in equation (26) includes the diagonal terms $s{=}s^{\prime}$ (in contrast, these terms are excluded from the summation in equation (12)). We highlight how the coherences’ conductivity $\kappa^{\alpha\beta}_{{\rm C}}$ remains exact (i.e. as in equation (12)) also in the context of the SMA approximation. This can be understood from equation (11), where it is apparent that the coherences are affected only by scattering events that yield a decrease of the off-diagonal ($s{\neq}s^{\prime}$) Wigner distributioncohen2004atom $\mathsfit{N}(\bm{R},\bm{q},t)_{s,s^{\prime}}$ (depopulation). Therefore, when the scattering events that yield an increase of the populations are neglected in the SMA approximation, the coherence equation — affected by depopulation only — is left unchangedSimoncelli_full_2018 and thus yields the same value for $\kappa^{\alpha\beta}_{{\rm C}}$ reported in equation (12).

As a final remark for this section, it is worth mentioning that the equations for $\kappa^{xx}_{\rm P}$ in the work of Fugallo et al.fugallo2013ab are recovered by choosing the eigenstates of $\mathrm{H}^{\rm har}$ in the degenerate subspaces to be eigenstates of $\bm{\mathsfit{V}}^{x}(\bm{q})$ and neglecting the term $\kappa^{\alpha\beta}_{\rm C}$.

The limit of a harmonic crystal with infinite unit cell: Allen-Feldman equation. Here we show that in the limit of a harmonic crystal with an infinite unit cell (i.e. in the case of a harmonic glass, whose Brillouin zone is reduced to the $\Gamma$-point), the term $\kappa_{\rm C}^{\alpha\beta}$ mathematically reduces to the Allen-Feldman expression for the thermal conductivity of harmonic glassesallen1989thermal; allen1993thermal:

[TABLE]

where $i{=}1,{\dots},3{\cdot}N_{\rm at}^{\rm tot}$ with $N_{\rm at}^{\rm tot}$ the total number of atoms in the (ideally infinite) single unit cell of the glass and $V_{\rm t}$ its (ideally infinite) volume. The mode specific heat $C_{i}(T)$ is:

[TABLE]

where \bar{N}_{i}{=}\big{(}{\exp\big{[}\hbar\omega_{i}/(k_{B}T)\big{]}{-}1}\big{)}^{-1} is the Bose-Einstein distribution and the mode diffusivity $D_{i}(T)$ is:

[TABLE]

$\mathsfit{V}^{\beta}_{\text{\tiny{AA}}}(\bm{q})_{s,s^{\prime}}$ is the velocity operator defined by Auerbach and AllenPhysRevB.29.2884, which is related to the velocity operator (8) derived in this work:

[TABLE]

The Dirac delta in equation (29) implies that these two velocity operators in equation (30) are equivalent in the computation of the Allen-Feldman diffusivities (29).

Equation (27) is derived under the hypothesis of extreme disorder — i.e. the phonon wave-packets do not to propagate far enough to sample the periodicity of the medium, rendering it impossible to assign them a wave vector or a group velocity. In practice, this corresponds to neglecting the thermal conductivity contributions coming from the diagonal elements of the velocity operator (otherwise, at the harmonic order, these would yield a divergent thermal conductivityfeldman1993thermal). Within such extreme disorder assumption, we are left with the coherences’ term and we will show now that $\frac{1}{3}\kappa_{C}^{\alpha\alpha}$ is equivalent to the Allen-Feldman formula (27). Let us start by writing explicitly part of equation (28):

[TABLE]

where we have used $V_{\rm t}{=}\mathcal{V}N_{\rm c}$. Taking first the limit $\Gamma(\bm{q})_{s}{\to}\eta\;\forall\hskip 2.84526pt(s,\;\bm{q})$ and then $\eta{\to}0$, we have that the coherences’ term in equation (12) becomes:

[TABLE]

where we have obtained the delta function as the limit of a Lorentzian for ${\Gamma(\bm{q})_{s}{\to}\eta{\to}0}$. It follows that, defining $\kappa_{A}{=}{\kappa^{\alpha\alpha}_{\rm dis}}/{3}$ and in the case of a Brillouin zone reduced to the $\Gamma$-point, we recover the Allen-Feldman formula:

[TABLE]

where $\omega_{s}{=}\omega(\bm{q}{=}\bm{0})_{s}$, $\bar{N}_{s}{=}\bar{N}(\bm{q}{=}\bm{0})_{s}$ and $\vec{\bm{\mathsfit{V}}}{=}\vec{\bm{\mathsfit{V}}}(\bm{q}{=}\bm{0})$.

Phonon linewidths distribution. Figure (3) shows the distribution of the phonon linewidths for CsPbBr3 at different temperatures.

The distribution spans a range of about 4 orders of magnitude (from about $10^{-2}$ cm*-1* to about $10^{2}$ cm*-1*); consequently, couplings between vibrational eigenstates having very different frequencies are allowed, as shown in figure (2d,e). Increasing the temperature yields an increase of the phonon linewidths, with consequent higher $\kappa^{\alpha\beta}_{\rm C}$ — since more couplings are allowed — and lower $\kappa^{\alpha\beta}_{\rm B}$.

Computational details. The experimental structure of orthorhombic CsPbBr3exp_cell_param with Pnma space group (Crystallographic Open DatabaseCOD_database id: 4510745) is first converted into an equivalent structure with Pbnm space group using VESTAmomma2011vesta, to allow a direct comparison with other references on this materialwang2018cation; Miyatae1701217. The first-principles equilibrium lattice parameters are then determined using Quantum ESPRESSOgiannozzi2017advanced, the PBEsol exchange-correlation functional and GBRVgarrity2014pseudopotentials pseudopotentials. The PBEsol functional is selected on the basis of its good agreement between first-principles and experimental lattice parameterslee2017ultralow; exp_cell_param (see table (1)) and its capability to accurately describe the vibrational properties of inorganic halide perovskiteslee2017ultralow. Kinetic energy cutoffs of 50 and 400 Ry are used for the wave functions and the charge density. The Brillouin zone is integrated with a Monkhorst-Pack mesh of $5{\times}5{\times}4$ points, with a (1,1,1) shift. Second-order force constants are computed on a $4{\times}4{\times}4$ mesh using density-functional perturbation theoryRevModPhys.73.515 as implemented in Quantum ESPRESSO. In order to obtain the third-order force constants using the PBEsol functional, the finite-difference method implemented in ShengBTEli2014shengbte is used, together with the interconversion software from ShengBTE to Quantum ESPRESSO available in the D3Q packagepaulatto2015first. Third-order force constants are computed on a $2{\times}2{\times}2$ mesh and considering interactions up to the third-nearest neighbors. All the force constants are Fourier interpolated on $8{\times}7{\times}5$ grid for the thermal conductivity calculations. In the latter, natural-abundance isotopic scatteringPhysRevLett.106.045901; tamura1983isotope is considered in addition to third-order anharmonicitycarrete2017almabte; fugallo2013ab; phono3py; paulatto2013anharmonic; PhonTS; alamode. From the force constants, the phonon frequencies and the full non-diagonal anharmonic scattering operator are computed as in previous workfugallo2013ab. The anharmonic scattering operator is then modified to account for the coherences’ term reported in equation (11). In figure (1), a smearing $\sigma$ of $2\;{\rm cm^{-1}}$ is used to compute the delta function appearing in the scattering operator. Results are converged with respect to the choice of this parameter: changing $\sigma$ to $1\;{\rm cm^{-1}}$ and $5\;{\rm cm^{-1}}$ produces populations, coherences and total thermal conductivities with a maximum relative deviation of $2.6\%$ from the data of figure (1). These latter refer to the single-mode relaxation time approximation (SMA); the exact thermal conductivity is always compatible within $0.7\%$ to its SMA value. Performing the thermal conductivity calculations using force constants that are Fourier interpolated on a $4{\times}4{\times}3$ mesh produces data compatible with the ones reported in figure (1), which use a $8{\times}7{\times}5$ mesh. In particular, the maximum relative deviation are respectively 1.8$\%$ for the populations in the SMA and 0.9$\%$ for the coherences’ thermal conductivity. As anticipated in the main text, the data shown in figure (1) for $\kappa^{\alpha\beta}_{\rm P}$ are compatible with the theoretical results presented in Ref.lee2017ultralow (and related Supplementary Material). In particular, the direction $xx$ in this work corresponds to the $[010]$ direction in Ref.lee2017ultralow. Experimental data from Ref.lee2017ultralow do not allow to evaluate the effect of boundary scattering, as they refer to nanowires having roughly the same size, and are smaller than all the nanowires used in Ref.wang2018cation. Therefore, we compared the bulk thermal conductivity predictions of equation (12) with the experimental data from Ref.wang2018cation. In the calculation of the data reported in figure (1), temperature effects are accounted through the Bose-Einstein distributions entering in the scattering operator. Accounting for the renormalization of the force constants due to temperature goes beyond the conceptual scope of this work but it has been addressed recently by several studiesPhysRevB.96.014111; PhysRevLett.120.105901; PhysRevB.94.020303; PhysRevB.98.024301; PhysRevLett.112.058501; PhysRevB.98.085205; PhysRevB.91.214310; aseginolaza2018phonon; PhysRevB.84.180301; PhysRevB.88.144301; PhysRevB.87.104111; PhysRevLett.100.095901. Calculations for diamond and silicon are performed using a $26{\times}26{\times}26$ q-point grid and a smearing of $8\;{\rm cm}^{-1}$ and $4\;{\rm cm}^{-1}$ respectively. Details on the first-principles calculation of the force constants for these materials can be found in Refs.fugallo2013ab; cepellotti2016thermal. In figure (2d,e), contributions to $\kappa^{xx}_{\rm{C,DOS}(2)}$ smaller than $4{\cdot}10^{-6}\;\mathrm{\frac{W{\cdot}cm^{2}}{m\cdot K}}$ are not reported for graphical clarity.

Data availability. Raw data were generated using the SCITAS High Performance Computing facility at the École Polytechnique Fédérale de Lausanne. Derived data supporting the findings of this study are available at https://archive.materialscloud.org.

Code availability. Quantum ESPRESSO is available at www.quantum-espresso.org; the scripts related to the computation of the third order force constants using the finite-difference method are available at bitbucket.org/sousaw/thirdorder; the D3Q package for Quantum ESPRESSO is available at sourceforge.net/projects/d3q/. The custom code developed in this work will be made available in a next release of the D3Q package.