Wigner formulation of thermal transport in solids

Two different heat-transport mechanisms are discussed in solids: in crystals, heat carriers propagate and scatter particle-like as described by Peierls' formulation of the Boltzmann transport equation for phonon wavepackets. In glasses, instead, carriers behave wave-like, diffusing via a Zener-like tunneling between quasi-degenerate vibrational eigenstates, as described by the Allen-Feldman equation. Recently, it has been shown that these two conduction mechanisms emerge from a Wigner transport equation, which unifies and extends the Peierls-Boltzmann and Allen-Feldman formulations, allowing to describe also complex crystals where particle-like and wave-like conduction mechanisms coexist. Here, we discuss the theoretical foundations of such transport equation as is derived from the Wigner phase-space formulation of quantum mechanics, elucidating how the interplay between disorder, anharmonicity, and the quantum Bose-Einstein statistics of atomic vibrations determines thermal conductivity. This Wigner formulation argues for a preferential phase convention for the dynamical matrix in the reciprocal Bloch representation and related off-diagonal velocity operator's elements; such convention is the only one yielding a conductivity which is invariant with respect to the non-unique choice of the crystal's unit cell and is size-consistent. We rationalize the conditions determining the crossover from particle-like to wave-like heat conduction, showing that phonons below the Ioffe-Regel limit (i.e. with a mean free path shorter than the interatomic spacing) contribute to heat transport due to their wave-like capability to interfere and tunnel. Finally, we show that the present approach overcomes the failures of the Peierls-Boltzmann formulation for materials with ultralow or glass-like thermal conductivity.


I. INTRODUCTION
In 1929 Peierls [1] formulated the phonon Boltzmann transport equation (BTE) to explain heat conduction in crystalline solids, envisioning that in crystals the microscopic heat carriers are phonon wavepackets that diffuse and scatter as if they were particles of a gas.Recent computational advances have allowed to compute the parameters entering in the linearized form of the BTE (LBTE) from first principles [2][3][4][5][6], and to solve it either approximately, in the so-called single mode approximation [7], or exactly, using iterative [8][9][10], variational [11], or exact diagonalization [6,[12][13][14] methods.Several studies [7,9,[15][16][17][18] have highlighted the accuracy of the LBTE in "simple crystals", i.e. crystals characterized by phonon interband spacings much larger than the linewidths.Notably, in absence of disorder the LBTE predicts the relation between the temperature and the thermal conductivity (κ) to follow a universal asymptotic decay κ(T )∼T −m for T larger than the Debye temperature, where m=1 when anharmonic three-phonon interactions are the dominant source of thermal resistance [19].Even faster decays (m>1) are possible when higher-order phonon scattering processes become relevant [20].This is in marked contrast to the much milder decrease in conductivity, or even the increase, observed in disordered or glassy materials as a function of temperature [21][22][23].An attempt to explain the microscopic heat-conduction mechanisms in glasses was made by Kittel in 1949 [24], who introduced a phenomenological model where atomic vibrations have a constant mean free path that is determined by the disorder length scale; however, this model lacks rigorous validation and has aroused some controversy [25][26][27][28], since it describes nonperiodic glasses through Peierls' theory for periodic crystals.A key step forward was made by Allen and Feldman in 1989 [26]: they envisioned that in disordered systems heat diffuses in a wave-like fashion, specifically through a Zener-like tunneling between quasi-degenerate vibrational eigenstates, with a tunneling strength dependent on the off-diagonal elements of the velocity operator [29,30].Notably, this formulation reproduces the temperature-conductivity curves measured in several glasses [26][27][28]31], which increase up to saturation in the high-temperature limit.
Recent work has shown that the LBTE fails to describe materials with ultralow thermal conductivity [21][22][23][32][33][34][35], leading to the speculation that wave-like transport mechanisms might emerge and coexist with particle-like conductivity [32,33,35,36].Conversely, particle-like propagation mechanisms have been suggested also for glasses (propagons, albeit without a formal justification) in order to rationalize experimental results [28].At variance with these formulations, the Wigner transport equation we recently introduced [37] naturally encompasses the emergence and coexistence of particle-like and wave-like conduction mechanisms, providing a unified ap-proach to heat-transport phenomena in solids, including crystals (where particle-like propagation dominates), glasses (where wave-like tunneling dominates), and all intermediate cases.
Here, we discuss the theoretical foundations for the derivation of the Wigner thermal transport equation from the phase-space formulation of quantum mechanics.We introduce a general theoretical framework, grounded on Wigner's formalism [38,39], that allows to describe transport phenomena in solids in a very general and convenient way.We employ this to describe thermal transport, showing how this formulation leads to a well defined microscopic energy field and related microscopic heat flux, and we use the latter to derive the thermal conductivity from the solution of the Wigner transport equation.We highlight a subtle aspect of this framework, leading to a well defined "Wallace" [40] phase convention that provides a formulation that is size-consistent and invariant with respect to any possible choice of a crystal's unit cell.We elucidate the physics underlying the crossover from particle-like to wave-like heatconduction mechanisms and provide a quantitative criterion to distinguish the regimes where particle-like or wave-like conduction mechanisms are dominant.Specifically, we show that a "Wigner limit in time" naturally emerges from this formulation as the time scale determining when particle-like and wave-like conduction mechanisms coexist and are equally important.We discuss how such a time scale is related to the gaps between the energy levels of atomic vibrations, or of the phonon bands in a crystal.We rely on these findings to introduce a regime diagram for thermal transport which allows to determine the theoretical framework needed to describe transport in a solid just from the knowledge of its vibrations' frequencies and lifetimes.Moreover, we demonstrate that atomic vibrations with a lifetime comparable to the Wigner limit in time have a mean free path of the order of the typical interactomic spacings, a length scale which is known in the literature as the "Ioffe-Regel limit in space" [28,41], and can be used to assess phenomenologically the validity of particle-like transport formulations (i.e.kinetic theories in which microscopic carriers propagate particle-like as in the Peierls-Boltzmann transport equation) [24,28,35,42].
The paper is organized as follows.In Sec.II we briefly summarize the key quantities needed to model thermal transport in dielectric crystals using the well known semiclassical Peierls-Boltzmann formulation.In Sec.III we introduce the quantities needed to describe thermal transport beyond the semiclassical Peierls-Boltzmann approach, presenting a second-quantization formalism for atomic vibrations in real space.We use this formalism to describe a system driven out of equilibrium by a temperature gradient, where the one-body density matrix is space-dependent and undergoes a Markovian irreversible evolution.In Sec.IV we discuss the connection between the density matrix and Wigner phase-space formalisms.We show that the Wigner formalism is par-ticularly convenient to describe transport phenomena in solids, and we employ it to derive the Wigner transport equation, which generalizes the Peierls-Boltzmann equation accounting not only for the particle-like propagation of phonon wavepackets but also for their wave-like tunneling between different bands.In Sec.V we use the Wigner formalism to derive an expression for the microscopic harmonic energy field and the related microscopic harmonic heat flux.We discuss thoroughly the linearized form of the Wigner transport equation (LWTE), showing how its solution determines the microscopic heat flux and a general LWTE thermal conductivity expression that accounts for both particle-and wave-like conduction mechanisms.We show that particle-like conduction mechanisms dominate in "simple crystals" having phonon interband spacings much larger than the linewidths, implying that in this regime the LWTE conductivity becomes equivalent to the Peierls-Boltzmann conductivity for weakly anharmonic crystals [1,43].In contrast, we show analytically that in the limiting case of a disordered harmonic solid, wave-like mechanisms are dominant and the LWTE conductivity becomes equivalent to the Allen-Feldman formula for the conductivity of glasses [26].Most importantly, we show that the LWTE conductivity covers in all generality all intermediate cases, including "complex crystals" having phonon interband spacings comparable or smaller than the linewidths and ultralow or glass-like conductivity.In Sec.VI we show that a phase convention (discussed e.g. in the book of Wallace [40]) needs to be employed in the derivation of the Wigner transport equation (in particular, in the computation of the dynamical matrix in the reciprocal Bloch representation and of the related off-diagonal velocity operator's elements).We report numerical calculations demonstrating that only with such convention the thermal conductivity is size-consistent and invariant with respect to the multiple choices for a crystal's unit cell.In Sec.VII we show that the present approach allows to predict the ultralow thermal conductivity of complex crystals used for thermal barrier coatings and thermoelectrics, with applications to the zirconate La 2 Zr 2 O 7 and the perovskite CsPbBr 3 as materials representative of these classes.In Sec.VIII we discuss how the LWTE predicts the coexistence of particle-like and wave-like heat-conduction mechanisms, with a crossover between phonons that mainly propagate particle-like and phonons that mainly tunnel wave-like.We show how phonons at the center of such crossover have a lifetime approximately equal to the Wigner limit in time, and a mean free path approximately equal to the Ioffe-Regel limit in space (i.e. the typical interatomic spacing).Finally, Sec.IX discusses the implementation in computer programs of the general LWTE thermal conductivity formula, providing numerical recipes apt to modify LBTE solvers [4,6,10,11,44,45] with minimal effort.

II. PRELIMINARIES
We start by considering a 3-dimensional bulk crystal, i.e. an infinite lattice with a basis having primitive vectors a α (α=1, . .., 3 is a Cartesian index), Bravais lattice vectors R = n α a α (n α ∈ Z ∀α), a primitive cell of volume V containing N at atoms at positions τ b ∈ V (b=1, . .., N at is an atomic label) [19,40,46].Here with "primitive cell" we mean as usual a minimal-cardinality set of atoms whose periodic repetition allows to describe the crystal; with "unit cell" a set of atoms, not necessarily of minimal cardinality, whose periodic repetition allows to describe the crystal.We want to describe here the ionic contribution to thermal transport, focusing on the evolution of the atomic vibrational energy in the presence of a temperature gradient.The quantum operator describing this is the nuclear Hamiltonian, taken here in the Born-Oppenheimer approximation where p(R) bα and û(R) bα are the momentum and displacement-from-equilibrium operators [47] for atom b along the Cartesian direction α in the primitive cell labeled by the Bravais vector R. M b is the mass of the atom b and V ({û(R) bα }) is the lattice-periodic interatomic potential, which depends on all the displacement operators.The usual canonical commutation relations are satisfied: In a solid, atoms oscillate around their equilibrium positions, allowing to approximate the potential with a nth order Taylor expansion of the displacement operators {û(R) bα }.For now it is sufficient to consider explicitly the leading (second-order, or harmonic) term in such expansion: where the zeroth-order term can be taken as reference for the energy (i.e.equal to zero) without loss of generality, the first-order term vanishes when evaluated at equilibrium [19], and O represents the higher-order terms in the displacements.We consider the atomic displacements from equilibrium to be small and below Lindemann's threshold [48,49], as it is the case in solids.Using perturbation theory one can decompose the Hamiltonian (1) as the sum of a leading Hamiltonian and a perturbation, Ĥ= Ĥhar + Ĥper , where Ĥhar is the leading Hamiltonian given by Eq. ( 1) limiting V to the harmonic term, and Ĥper is the perturbation in the potential due to anharmonic third-or higher-order terms in Eq. (3).It is worth mentioning that also the presence of isotopic mass disorder can be treated as a kinetic-energy perturbation and taken into account in Ĥper [7,50].We first focus on the description of the leading Hamiltonian Ĥhar ; the perturbation introducing transitions between the eigenstates of Ĥhar will be discussed later.It is useful to recast the harmonic coefficients in Eq. ( 3) in terms of the tensor of mass-renormalized interatomic force constants commonly employed in the literature [40,51], which is symmetric and translation invariant: In fact, as extensively discussed in textbooks [40,46], at the leading harmonic order atomic vibrations in crystals can be decomposed in a linear combination of the normal modes of vibration (phonons), and the properties of these normal modes can be obtained computing the Fourier transform of the mass-renormalized interatomic force constants (4), hereafter referred to as "dynamical matrix in the reciprocal Bloch representation" where q is a wavevector that can be restricted to the first Brillouin zone of the crystal B, and the sum in Eq. (7) runs over a single lattice vector due to translation invariance (6) [40].The eigenvalues ω 2 (q) s , where s is a vibrational-mode index running from 1 to 3 N at , and eigenvectors E(q) s,bα of the dynamical matrix (7) b α D(q) bα,b α E(q) s,b α = ω 2 (q) s E(q) s,bα , provide the vibrational frequencies ω(q) s and displacement patterns of the normal modes (we recall that E(q) s,bα describes how atom b moves along the Cartesian direction α when the phonon with wavevector q and mode s is excited).In Eq. ( 7) the dynamical matrix is computed using the "smooth phase convention" employed e.g. in the book of Wallace [40], where the phases depend on the atomic positions in real space (R+τ b −τ b ).It is worth mentioning that a different "step-like" phase convention is often employed in the literature to compute the dynamical matrix, where the phases depend only on the Bravais lattice vector R and thus vary discontinuously (as explained in the following, see also the book of Ziman [19]), Hereafter we will refer to the matrix D(q) bα,b α (Eq.( 7)) as the "smooth dynamical matrix", and to the other matrix D(q) bα,b α (Eq.( 9)) as the "step-like dynamical matrix".It can be shown analytically that these two dynamical matrices are related by a unitary transformation as such, they have the same eigenvalues ω 2 (q) s and their eigenvectors are related by the unitary transformation While both these conventions allow to describe the equilibrium vibrational properties of materials (e.g. the phonon spectrum), we will show later that the smooth phase convention ("Wallace") has to be used to describe the out-of-equilibrium propagation of vibrations within the Wigner framework discussed here.
The derivative with respect to wavevector of the smooth dynamical matrix ( 7) is related to the propagation (group) velocity v β (q) s of the phonon wave packet centered at wavevector q and having mode s, which in absence of degeneracies (ω(q) s =ω(q) s for s =s ) is [11,52] In the presence of degeneracies (ω(q) s =ω(q) s for s =s ), the eigenvectors appearing in Eq. ( 12) should be chosen so that the matrix ∂q β E(q) s ,b α is diagonal in the indexes s, s corresponding to degenerate frequencies [11].We will show later that in the Wigner framework quantities related to the non-degenerate off-diagonal elements of the matrix m β (q) s,s become relevant.In this regard, it is possible to highlight already now that the smooth phase convention appears to be more suitable for the computation of these off-diagonal elements.For example, the smooth dynamical matrix (7) is not changed if using two primitive cells that differ just by one atom rigidly shifted by a Bravais-lattice vector R .This can be verified by replacing τ b → τ b +R in Eq. (7), and exploiting translation invariance.As a consequence, also the eigenvectors of the smooth dynamical matrix (8), and thus the off-diagonal elements of the matrix m β (q) s,s , do not vary.In contrast, the step-like dynamical matrix (9) and its eigenvectors (11) do not satisfy this invariance; applying the same shift to Eq. ( 9) leads to a different result.
The quantities described up to now enter into the well known semiclassical model for thermal transport in crystals developed by Peierls [1,43], in which the dynamics of out-of-equilibrium and space-dependent populations of phonon wavepackets (n(r, q, t) s , here r is a continuous Cartesian coordinate) is determined by the Boltzmann transport equation: where ∂n(r,q,t)s ∂t Ĥper is the phonon scattering superoperator that originates from Ĥper [11,19] and that will be discussed later.The Peierls-Boltzmann equation ( 13) is derived under the assumption that phonon wavepackets drift akin to the particles of a classical gas in the presence of a temperature gradient [43].Mathematically, this is apparent from the left-hand side of Eq. ( 13), which describes phonons' drift and has a mathematical form analogous to the drift term of the Boltzmann equation for the particles of a classical gas [53].In the next sections we discuss a formalism that allows to generalize Peierls' model beyond this particle-like behavior, and we discuss the conditions under which the Peierls-Boltzmann equation (13) fails and this generalization is relevant.

III. QUANTUM DESCRIPTION OF ATOMIC VIBRATIONS
In this section we discuss a second-quantization formalism for atomic vibrations in real space that is suitable to define the one-body density matrix of an outof-equilibrium system, where a temperature gradient enforces a space-dependent atomic vibrational energy.

A. Second quantization for atomic vibrations
The leading Hamiltonian Ĥhar is quadratic in the momentum and displacement operators, and thus represents a many-body system of harmonic oscillators.With the goal of exploiting the second-quantization formalism to simplify the description of this many-body problem, and tracking the space-dependent atomic vibrations that characterize the regime where a temperature gradient drives thermal transport (i.e. the regime in which atoms vibrate more in the warm region and less in the cold region), we introduce the class of bosonic operators in real space â(R) bα that are constructed so that â(R) bα annihilates a vibration centered around atom b in the unit cell R and along the Cartesian direction α.In Eq. ( 14), the matrices . These annihilation-of-localized-vibrations operators (14) are labeled by R, b, α, and satisfy the bosonic commutation relations together with their adjoints (the creation-of-localized-vibration operators â † (R ) b α ) [54] [ To better understand the physical action of the localized bosonic operators in ( 14) and ( 15), we proceed as follows: first, we create a vibrational excitation on the ground state |0 using the creation operator at position R+τ b , |ψ Rbα =â † (R) bα |0 ; then, we look at how the atom at the generic position R +τ b vibrates in such an excited state, computing the expectation value of the squared atomic displacement operator û2 (R ) b α .Specifically, from Eq. ( 14) it follows that the displacement-fromequilibrium operator for an atom can be written as: Then, the expectation value of the squared atomic displacement on the excited state is: where G −1 R b α ,Rbα goes to zero in the same limit [55].Therefore, Eq. (17) shows that the bosonic operator â † (R) bα creates atomic vibrations along direction α and centered around the position R+τ b .In other words, after the action of â † (R) bα on the ground state, atoms close to R+τ b have a larger mean square displacement than atoms far from R+τ b (Fig. 1).We note that an analogous calculation can be done for the expectation value of the square momentum operator; such a calculation is reported in Appendix A 1 and shows that considerations analogous to those for the mean square displacement reported here hold also for the mean square momentum.We have thus shown how these bosonic operators (14) can be used to describe space-dependent atomic vibrations.
The localized bosonic operators (14) allow also to write the leading harmonic Hamiltonian in second-quantized form as where the additive constant E 0 = 2 Rbα √ G Rbα,Rbα represents the zero-point energy, which can be taken as reference energy and thus will be omitted in the following.We see from Eq. ( 18) that the harmonic vibrational energy is a one-body operator; therefore, to compute its expectation value the knowledge of the complex manybody density operator ρ is not needed, but it is sufficient to know the much simpler one-body density matrix where Tr denotes the trace operation over the Fock space.In fact, it is possible to show that Up to now we have managed to reduce the complexity of the many-body problem by introducing the localized bosonic operators in real space (14) and the one-body density matrix (19).Now we show that the formulation can be further simplified exploiting the invariance under translation of the G tensor in (4).We start by defining the Fourier transform â(q) bα of the localized bosonic operator â(R) bα where are the Fourier transforms of the fourth root of the mass-rescaled harmonic interatomic force-constant matrix and of its inverse.It is possible to verify that, analogously to the matrix 4  √ G , also the fourth root of the smooth dynamical matrix in reciprocal Bloch representation satisfies , and an analogous calculation allows to obtain the smooth dynamical matrix (7) from its square root.Analogous considerations hold for the roots of the inverse smooth dynamical matrices in the reciprocal Bloch representation.
Combining Eq. ( 15) and Eq. ( 21), it is possible to show that also the operators in reciprocal space (21)  bosonic commutation relation where we have highlighted that the Dirac delta is obtained from the Kronecker delta in the bulk limit, i.e. when the crystal is made up of N c →∞ Bravais-lattice sites [57][58][59].
The key role of the bosonic operators in reciprocal space (21) is their capability to simplify to block-diagonal form the representation of translation-invariant quantities such as the harmonic Hamiltonian (18): So, the bosonic operators â(R) bα , â † (R) bα , localized in real space, allow to discuss local vibrational excitations, and their Fourier transforms â(q) bα , â † (q) bα allow to describe translation-invariant quantities in block-diagonal form in reciprocal space.
To compute the expectation value of the harmonic vibrational energy from Eq. ( 23), it is sufficient to know the one-body density matrix in reciprocal space, (q, q , t) bα,b α =Tr ρ(t)â † (q ) b α â(q) bα . ( To describe a weakly non-homogeneous out-ofequilibrium system, in which a small temperature gradient is present in real space, it is convenient to work in reciprocal space with the average and displacement coordinates (i.e.given a generic couple of wavevectors q 1 and q 2 , they are written in terms of the displacement q =q 1 −q 2 from their average q= q1+q2
In the context of Eq. (25) this implies that the one-body density matrix in reciprocal space q+ q 2 , q− q 2 , t bα,b α is close to the diagonal form, i.e. appreciably different from zero only for |q |→0.We will discuss later how this property can be exploited to simplify to solvable form the equation describing the system's evolution.
The transformation ( 25) is a reminder that the onebody density matrix in direct and reciprocal space contain the same information.Heat transport in solids is driven by a temperature gradient in real space, thus the formulation in real space is convenient to keep track of the vibration's locations and relate them to the spacedependent temperature driving transport.Conversely, we have seen that the formulation in reciprocal space simplifies to block-diagonal form the expression of key quantities for the description of thermal transport (e.g. the harmonic Hamiltonian ( 23)), and we will show later that using Wigner's phase-space framework it is possible to combine the advantages of these two formulations in real and in reciprocal space.
We note in passing that the standard phonon operators [19,40] where s is the phonon-band index, have a real-space representation that is not suitable to track where vibrations are centered.In fact, the eigenvectors of the dynamical matrix in the reciprocal Bloch representation E(q) s,bα (Eq.( 8)) have an undetermined phase, which allows to apply arbitrary unitary transformations to â(q) s and thus shift by an arbitrary Bravais lattice vector the corresponding real-space operator â(R) s = V (2π) 3 B â(q) s e +iq•R d 3 q.Details are discussed in Appendix A 2, where it is shown that the non-uniqueness of the real-space representation of the phonon operators â(q) s mirrors the electronic case, where the phase indeterminacy of the Bloch orbitals is reflected in the nonuniqueness of their transformation into Wannier functions.

B. Evolution equation for the density matrix
We recall that the goal of the present work is to describe the energy transfer across a solid material driven by a temperature gradient.More precisely, we consider the regime where temperature variations are appreciable over a macroscopic length scale L that is much larger than the interatomic spacing a.We assume that such a temperature gradient derives from a continuous energy exchange with some heat baths.This implies that the system undergoes an irreversible evolution [61], which can be modeled as Markovian and described by the master equation [62][63][64] where the scattering superoperator ∂ ρ(t) ∂t Ĥper is determined by the perturbative Hamiltonian Ĥper (which accounts e.g. for the presence of anharmonicity and isotopes) and describes transitions between the eigenstates of the leading Hamiltonian Ĥhar .
Recalling Eq. ( 23), we have that Ĥhar is a one-body operator; thus, applying the operator â † (q− q 2 ) b α â(q+ q 2 ) bα to both sides of Eq. ( 27) and taking the trace as in Eq. ( 24), gives us the evolution equation for the one-body density matrix We also note that Eq. ( 28) can be derived starting from real space [55], applying the operator â † (R ) b α â(R) bα to both sides of Eq. ( 27), then taking the trace as in Eq. ( 19), and finally performing the Fourier transform (25).
We recall that Eq. ( 25) implies that at equilibrium, where the system is homogeneous (translation-invariant), the one-body density matrix is diagonal in reciprocal space, i.e. q+ q 2 , q− q 2 , t bα,b α = q, q, t bα,b α V δ(q ).As a consequence, Eq. ( 28) yields a non-trivial evolution only out of equilibrium, where a perturbation (temperature gradient) in real space implies that the one-body density matrix has non-zero off-diagonal elements.In the next section we will discuss a theoretical framework that makes the description of this non-homogeneous regime particularly intuitive and manageable.

IV. PHASE-SPACE FORMALISM
The ideal framework to describe heat transport in solids should keep track of: (i) the equilibrium eigenstates towards which the system relaxes, for which the quasimomentum q is a good quantum number; (ii) the spatial dependence of the perturbation driving transport, whose localization in real space would be conveniently described by an Hilbert-space basis having as quantum number the Bravais lattice vector R.However, the quasi-momentum and Bravais-lattice operators ( q and R, respectively) are a pair of canonically conjugate operators [65], and thus their eigenvalues ( q and R, respectively) cannot be used to label simultaneously a quantum-mechanical representation within the usual Dirac formalism.The Wigner phase-space formalism [38,39,[66][67][68] allows to describe quantum mechanics in terms of distributions having as arguments eigenvalues of canonically conjugate operators, and in this section we discuss the application of such formalism to solids.We show that the central quantity appearing in this formulation generalizes the semiclassical distribution appearing in the Peierls-Boltzmann equation [43], describing not only intraband propagation of particle-like phonon wavepackets, but also wave-like interband (Zener) tunneling of phonons.

A. Wigner transform
As anticipated, the central objects of the Wigner formalism are phase-space distributions having as arguments eigenvalues of non-commuting operators [39,68,69].In order to show how such a framework allows to describe conveniently transport in non-homogeneous, out-of-equilibrium systems, it is useful to start from the Wigner transformation that maps the matrix elements of a one-body operator O q+ q 2 ,q− q 2 bα,b α into a "phasespace" distribution that depends on a wavevector q (belonging to a Brillouin zone B) and a Bravais lattice vector R.More generally, we will discuss soon that the spatial dependence of a Wigner distribution can be extended from the discrete Bravais lattice vectors R to a continuous Cartesian position r; for the sake of generality we use this definition then for the Wigner transformation W [O] of a one-body operator O: When the Wigner transformation ( 29) is applied to the one-body density matrix (24), it becomes apparent that the Wigner distribution w (r,q, t) bα,b α =W [ ] (r,q, t) bα,b α does not depend on space if the system is homogeneous (i.e.translation-invariant; as discussed before, Eq. ( 25) implies that for a translation-invariant system the onebody density matrix is nonzero only for q =0).Conversely, for an out-of-equilibrium, non-homogeneous system the one-body density matrix is non-diagonal in q and thus its Wigner representation depends on space.
In general, the position r appearing in Eq. ( 29) can be any (continuous) position in real space.Nevertheless, restricting such a position to the discrete Bravais-lattice vectors is sufficient to obtain a phase-space distribution that contains the same information of the Dirac matrix element of the corresponding operator.In fact, knowing the phase-space distribution at the discrete Bravais vectors (r=R) allows to determine the matrix element of the corresponding operator (proof in Appendix B) We note in passing that the time dependence for O and W [O] in Eqs.(29,30) has been reported for generality; in the present paper the only operator (phase-space distribution) that can depend on time is the density matrix (Wigner distribution w (R, q, t) bα,b α ).
In order to shed light on why the Wigner framework is particularly convenient to describe transport, it is useful to discuss how expectation values are computed in this framework, and compare it the usual Dirac formalism.We start by recalling that in the Dirac formalism the expectation value of a quantum mechanical observable Â, on a state represented by the density matrix ρ(t) is determined from the trace operation Tr ρ(t) Â = Â .In the Wigner framework, instead, expectation values are computed as phase-space-like integrals of the product between the Wigner representations of the density matrix and that of the observable [68].For the one-body operators and phase-space distributions in focus here and related via Eq.( 29) these two equivalent methods to compute expectation values are summarized in the following equation (proof in Appendix B): The third line in Eq. ( 31) contains a sum over all lattice vectors in real space (R) and an integral over wavevectors in reciprocal space (q); therefore, carrying out only the integration in reciprocal space yields an expression for the expectation value as a spatial integration of a spacedependent quantity (here with "spatial integration" we mean the sum over R).This shows that the Wigner framework provides prescriptions to resolve in space expectation values [68,70].Explicit calculations will be reported later, when we will exploit Eq. ( 31) to define a space-dependent vibrational energy field, compute the related heat flux, and thus determine the thermal conductivity.

B. Wigner transport equation
The first key step in obtaining the equation describing the temporal evolution of the Wigner distribution at the Bravais site R and wavevector q, w (R, q, t) bα,b α is applying the transformation (29) to the evolution equation for the one-body density matrix (28).Then, we note that we are interested in the close-to-equilibrium regime characterized by weak non-homogeneities in real space, i.e. in the regime where temperature varies slowly and causes spatial variations of the one-body density matrix (19) to be appreciable only at "mesoscopic scales" l that are much larger than the lengths a at which atomic positions can be resolved (l a).We start by recalling that as the spatial variations of the one-body density matrix in real space become smaller, the one-body density matrix in reciprocal space becomes closer to being diagonal in the arguments (see Eq. ( 25) and related discussion).This implies that, close to equilibrium, the one-body density matrix q+ q 2 ,q− q 2 bα,b α is sharply peaked around q, and significantly different from zero only for |q | 2π/|a i | (a i is the i-th direct lattice vector).This property allows us to simplify the evolution equation (28), performing a Taylor expansion for |q |→0 of the coefficients appearing in that equation: Inserting such approximation into Eq.( 28), one obtains a simplified evolution equation: Multiplying all terms in Eq. ( 33) by V (2π) 3 e iq •R and integrating q over the Brillouin zone yields the evolution equation for the Wigner distribution where we have introduced the notation [ D(q), w (R, q, t)] bα,b α to denote the matrix element of the commutator, i.e.
the quantity b ,α ; an analogous notation is employed for the anticommutator {, }.Moreover, the derivative of the Wigner distribution with respect to the position R is obtained taking first the derivative with respect to the continuous position r (see Sec. IV A), and then evaluating such derivative at the Bravais lattice vector R.
After having performed the first-order Taylor expansion (32), for which localization and differentiation properties are crucial, one can apply any unitary transformation to the simplified evolution equation (34), leaving the physics unchanged.Therefore, we can apply to Eq. ( 34) the unitary transformation (8) that diagonalizes the square root of the dynamical matrix in the reciprocal Bloch representation, obtaining distributions that depend on position R, wavevector q, and phonon band indexes (s, s ), and thus can be compared directly with the distribution appearing in the Peierls-Boltzmann equation.Specifically, we introduce the quantities (36) which generalize the concepts of phonon populations and group velocities beyond the particle-like interpretation provided by the semiclassical Peierls-Boltzmann equation (13).In fact, in the absence of degeneracies (i.e. for s =s , ω(q) s =ω(q) s ) their diagonal elements (s=s ) coincide with the phonon populations (n(R, q, t) s,s = n(R, q, t) s ) and group velocities (12) (v β (q) s,s =v β (q) s ) appearing in the Peierls-Boltzmann equation ( 13), respectively, while their off-diagonal elements emerge from the wave-like nature of atomic vibrational eigenstates and are related to the phase coherence between pairs of vibrational eigenstates s and s [62,64,71,72].
From now on we will use the textbook nomenclature [73] and refer to the diagonal elements (s=s ) of the Wigner distribution (35) as "populations" and to the off-diagonal ones (s =s ) as "coherences".As discussed in the literature [61, 62,64,71,72], populations have a well defined energy (frequency) and therefore can be interpreted as particle-like excitations with a well defined wavevector q and mode index s.In contrast, coherences do not have an absolute energy and cannot be related to a single eigenstate; they describe oscillations between pairs of eigenstates and correspond to an evolution which does not preserve the nature of the excitation (e.g.we will show later that coherences describe couplings between two different vibrational modes s, s at the same wavevector q).Moreover, the velocity operator (36) has the following properties: it is Hermitian (since it is the incremental ratio of an Hermitian operator), its different Cartesian components do not commute in general, s (v α (q) s,s v β (q) s ,s −v β (q) s,s v α (q) s ,s ) =0 if α =β, and from the time-reversal symmetry of the dynamical matrix in reciprocal Bloch representation [51] it follows that v β (q) s,s =−v β (−q) s ,s .
Representing Eq. ( 34) in this phonon basis, and denoting with ∂n(R,q,t) s,s ∂t Ĥper the scattering superoperator in the basis of the phonon eigenmodes (which will be discussed later), we obtain the Wigner transport equation where the terms in the square bracket are reminiscent of the commutator appearing in the quantum evolution equation ( 28); the time derivative and the terms in curly brackets generalize the drift term of the LBTE (13), since the diagonal elements (s=s ) of these terms are equal to the product between the phonon group velocities and phonon populations appearing in the LBTE (this will be discussed more in detail later in Sec.V).The Wigner scattering superoperator in Eq. ( 37) is derived from the many-body scattering superoperator of Eq. ( 27) following a textbook approach [64].Here, we account for scattering due to third-order anharmonicity and perturbative mass disorder [50]; we neglect energy-renormalization effects [62], and we follow the standard procedure of treating scattering at linear order in the deviation from equilibrium [19].The result is Here, N(q Bose-Einstein distribution at the equilibrium temperature T , and Ω in (q, q ) s,s is the scattering matrix that describes the repumping of the populations due to the incoming scattered phonons (Ω in (q, q ) s,s =−A in qs,q s / N(q ) s [N(q ) s +1] , where A in qs,q s is defined by the anharmonic and mass-disorder contributions in Eq. ( 15) of Ref. [11]).
Γ(q) s is the linewidth of the phonon with wavevector q and mode s (related to the phonon lifetime τ (q) s via [Γ(q) s ] −1 =τ (q) s , see e.g.Eq. ( 14) of Ref. [11] for the definition of the phonon lifetime τ (q) s ), and accounts for all the scattering events that decrease the populations or coherences involving phonons with wavevector q and mode s.The conservation of energy in anharmonic and mass-disorder scattering [74] implies the following relation between the "depumping rate" Γ(q) s and the "repumping matrix" Ω in (q, q ) s,s From Eq. ( 38) it is evident that populations scatter only with populations, and coherences scatter only with coherences.Importantly, Eq. (38) shows that scattering does not create but only destroys off-diagonal coherences s =s (since the scattering term for coherences in the third line of Eq. ( 38) has negative sign), while it can both create (term with plus sign in the second line of Eq. ( 38)) or destroy (term with minus sign in the first line of Eq. ( 38)) populations.It follows that scattering drives the Wigner distribution towards a distribution that is diagonal in the phonon-mode indexes, which coincides with the local-equilibrium state towards which the Peierls-Boltzmann transport equation evolves to [74,75].We conclude by noting that the scattering superoperator has been linearized for convenience in the deviation from equilibrium, to report equations containing quantities available in LBTE solvers [4,6,10,11,44,76], but it can be extended to the non-linear order to describe thermal transport in the far-from-equilibrium regime [19,64].Moreover, the description of scattering employed here relies on the assumption that phonons are well defined excitations, i.e. their linewidth is small compared to their energy, γ(q) s < ω(q) s [77].When this assumption breaks down, it is no longer possible to describe scattering in terms of the phonon wavevector q and mode s as in Eq. ( 38), and one has to consider phonon spectral functions as discussed in Refs.[78,79].

V. THERMAL CONDUCTIVITY A. Vibrational energy field and flux
We have seen that the Wigner framework provides prescriptions to resolve in space expectations values (Eq.( 31) and related discussion).This property allows to compute the expectation value of the harmonic vibrational energy of a crystal integrating the space-dependent Wigner energy field.To see this, we first recast the expectation value of the harmonic Hamiltonian (20) in reciprocal space using the inverse of Eq. ( 25): Combining Eq. ( 40) with Eq. ( 31) allows to compute the expectation value of the harmonic vibrational energy as spatial integration of Wigner's energy field E(R) [70], , where The first line of Eq. (41) shows that the Wigner representation of the translation-invariant square root massrenormalized force-constant matrix √ G Rbα,R b α does not depend on space and is equal to the square root of the dynamical matrix in reciprocal Bloch representation D(q) bα,b α (see Appendix B for details).The second line of Eq. ( 41) has been obtained recasting the first line in the phonon eigenmodes basis (see Eq. ( 8) for the definition of the phonon eigenmodes basis).
The time derivative of the Wigner energy field E(R), which is intensive, can be related directly to the gradient of the intensive heat flux J (R, t).Using Eq. ( 37) to evaluate the time derivative of Eq. ( 41) one obtains a continuity equation that relates the time derivative of the energy field to the heat flux J (R, t): is the heat flux that originates from the leading harmonic Hamiltonian Ĥhar , while σ(R, t) originates from the perturbation Ĥper .In the following σ(R, t) will be neglected, since with respect to ∇ • J (R, t) it is of higher order in the perturbative treatment of anharmonicity.Importantly, the spatial average of the harmonic Wigner heat-flux (42) differs from the well-known harmonic heat-flux derived by Hardy [29].Such a difference is discussed in detail in Ref. [78], and originates from a difference in the expressions for the harmonic microscopic energy field used as starting point for the computation of the heat flux here (Eq.( 41)) and in Hardy's work [29].

B. Steady-state solution of the LWTE
The thermal conductivity is defined as the tensor κ αβ that relates a temperature gradient to the heat flux generated in response to it: J α =− β κ αβ ∇T β .Consequently, one can determine the thermal conductivity as follows: (i) fix a temperature gradient constant in time; (ii) solve Eq. ( 37) with boundary conditions corresponding to the temperature gradient in (i); (iii) determine the heat flux inserting the solution obtained in (ii) into Eq.( 42); (iv) determine the thermal conductivity as a tensor which relates the temperature gradient (i) and the heat flux (iii).In this section we will determine the thermal conductivity following the protocol above.We consider a system where a steady and spacedependent temperature T (r) T +∇T •r is enforced along a certain direction.In addition, the temperature gradient is assumed to be small and to be related to temperature variations appreciable over a length scale L much larger than the interatomic spacing a, i.e. ∇T − δT L , where δT T .These considerations allow us to look for a steady-state solution for Eq. ( 37) as a perturbation of order ∇T in the local-equilibrium distribution corresponding to the local temperature T (r) (see e.g.Ref. [80] for the mathematical details).Therefore, we look for a solu-tion of the form n(R, q) s,s =N(q) s δ s,s +n(R, q) s δ s,s +n (1) (q) s,s •∇T, (43) where N(q) s is the equilibrium Bose-Einstein distribution at temperature T , n(R, q) s δ s,s = d N(q)s d T (T (R)− T )δ s,s is the distribution that describes the local-equilibrium state corresponding to a space-dependent temperature (such a local-equilibrium distribution has been discussed e.g. in Refs.[81,82]), and n (1) (q) s,s •∇T is in general nondiagonal in s, s 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 work [11,13]).We recall that, even though the phase-space distributions (29) are defined everywhere in space, only the values of these distributions at the Bravais lattice sites appear in Eq. ( 37) and Eq. ( 43), since the knowledge of the Wigner distributions at these points is sufficient to fully describe the problem, as explained in Sec.IV A.
Inserting the expansion (43) into the LWTE (37) at steady state, and considering only terms linear in the temperature gradient, we obtain a matrix equation that is decoupled in its diagonal and off-diagonal parts.The (diagonal) equation for the populations is the usual steady-state LBTE [11]: where on the left-hand side of this equation we have used the property that the local-equilibrium distribution n(R, q) s δ s,s is an eigenvector with zero eigenvalue of the populations' scattering integral (the sum of the first two terms in Eq. (38), Ω tot (q, q ) s,s =Ω in (q, q ) s,s −(2π) 3 Γ(q) s δ(q−q )δ s,s ) [14,74,81]; thus, s B Ω tot (q, q ) s,s n(R, q ) s d 3 q =0.Eq. ( 44) can be solved exactly using iterative [6,10,44], variational [11] or exact diagonalization methods [12,13].
The off-diagonal equation for the coherences (s =s ) is and its solution reads (s =s ): We note that solving exactly the equation for the populations (44) requires accounting for all the interactions between populations (those described by the scattering matrix Ω tot (q, q ) s,s ); in contrast, the exact evolution equation for the coherences (45) contains only the coupling of a coherence n(R, q) s,s (with s =s ) to itself.From a practical point of view, the computational cost for solving exactly the coherences' equation ( 45) is negligible compared to that for solving exactly the LBTE (i.e. the populations' equation ( 44)) [83].
After having determined the steady-state solution of the LWTE (43) (i.e. the sum of the solutions of Eq. ( 44) and of Eq. ( 45)), one can insert it into Eq.( 42) and thus evaluate the heat flux.At this point it is possible to show that only the deviation-from-equilibrium n (1) (q) s,s •∇T of the solution (43) gives a non-zero contribution to the heat flux, since the Bose-Einstein term (N(q) s δ s,s ) and the local-equilibrium term (n(R, q) s δ s,s ) are both even functions of the wavevector; thus, when multiplied by the diagonal elements of the velocity operator -which have odd parity in the wavevector -they yield an oddparity function whose integral over the symmetric Brillouin zone yields zero [14].Therefore, one obtains a linear relation between heat flux and temperature gradient, J α =− β κ αβ ∇ β T , with the proportionality tensor being the thermal conductivity: where is the specific heat of a phonon with wavevector q and mode s, κ αβ P is the "populations conductivity" obtained from the exact (i.e.iterative [8,10], variational [11], or exactdiagonalization [6,12,13]) solution of the diagonal (population) part of the LWTE (44) -which is exactly the Peierls-Boltzmann equation for phonon wavepackets [43] -and the additional tensor derives from the equation for the coherences (45), thus it is called here "coherences conductivity" (κ αβ C ). Eq. ( 47) can be explicitly calculated by evaluating the integral over the Brillouin zone on a uniform mesh of wavevectors, i.e. replacing 1 (2π) 3 B d 3 q → 1 VNc q , where N c is the number of qpoints that sample the Brillouin zone and over which the sum runs, and V is the volume of the primitive cell [59].
We stress that the LWTE conductivity (47) derives from an exact solution of the LWTE (37) and thus is also valid in the hydrodynamic regime, where the scattering between phonons that conserve the crystal momentum ("normal processes") dominate and yield a Peierls-Boltzmann thermal conductivity that is very high (κ P (300 K) > ∼ 10 3 W/mK) and several order of magnitude larger than the coherences conductivity [14].It is also worth mentioning that in the case of low-thermalconductivity solids, where "Umklapp" scattering between phonon wavepackets that dissipate the crystal momentum are dominant, Eq. ( 38) contains depumping (negative sign) terms that are dominant with respect to the repumping (positive sign) terms.Therefore, one can apply the single-mode relaxation-time approximation (SMA) [11,84], which consists in discarding the repumping terms in Eq. (38).In practice, this approximation implies discarding the second line (the repumping term) in the populations' equation (44), and leaves the coherences' equation (45) unchanged.It follows that within this approximation nothing changes for the coherences, while the populations equation assumes a simpler form that can be solved at a greatly reduced computational cost compared to the full solution.Therefore, solving the populations' equation simplified with the SMA, one obtains: When the approximated SMA solution for the populations ( 48) is inserted into the expression for the heat flux (42), the populations' conductivity (term κ αβ P in Eq. ( 47)) becomes exactly equal to the SMA expression for the Peierls-Boltzmann conductivity [59]: Eq. (49) shows that the populations' conductivity can be interpreted in terms of particle-like phonon wavepackets labeled by (q) s , which carry the specific heat C(q) s and propagate between collisions over a length Λ α (q) s =v α (q) s,s [Γ(q) s ] −1 (as mentioned before, v α (q) s,s is the propagation (group) velocity in direction α, and [Γ(q) s ] −1 is the typical inter-collision time).Importantly, we note that κ αβ P can be interpreted in terms of microscopic carriers propagating particle-like without relying on the SMA approximation, since it has been shown that the conductivity deriving from the exact solution of the population's equation ( 44) can be determined exactly and in a closed form as a sum over relaxons [13,14], i.e. particle-like collective phonon excitations that are the eigenvectors of (a symmetrized version of) the full LBTE scattering matrix Ω tot (q, q ) s,s discussed after Eq. (44).
The microscopic conduction mechanisms that determine the term κ C in the conductivity expression (47) is that of phonons tunnelling between two different bands at the same wavevector q.Such an interband tunneling mechanism originates from the wave-like nature of phonons, and in general becomes stronger as the co-herence between two phonons (q) s and (q) s becomes stronger (i.e.their frequency difference ω(q) s −ω(q) s becomes smaller).We note in passing that this mechanism has some analogies with the Zener interband tunnelling of electrons discussed e.g. in Refs.[85][86][87][88][89] relying on the density-matrix formalism.
It has been shown in Ref. [14] that in "simple" crystals, i.e. those characterized by phonon interband spacings much larger that the linewidths (e.g.silicon, diamond), particle-like mechanisms dominate and thus κ αβ P κ αβ C .We will discuss in Sec.VII how in the "complex" crystal regime, characterized by phonon interband spacings smaller that the linewidths, particle-like and wavelike mechanisms coexist and are both relevant, implying Finally, we note that the Wigner formulation can be applied also to amorphous solids, describing them as limiting cases of disordered but periodic crystals in the limit of infinitely large primitive cells (N at →∞ and consequently V→∞).Importantly, it is possible to show analytically that in the ordered limit describing a harmonic glass (i.e.first letting the primitive cell volume go to infinity V→∞, and then letting each linewidth go to the same infinitesimal broadening η, Γ(q) s → η→0, ∀ q, s [26,27]) the normalized trace of the LWTE conductivity (47) becomes equivalent to the Allen-Feldman formula for the conductivity of glasses [37].More generally, the Wigner formulation extends Allen-Feldman theory accounting also for the effects of anharmonicity on thermal transport in disordered solids, since Eq. ( 47) accounts for anharmonicity thorough the linewidths.In practice, evaluating the LWTE conductivity (47) in disordered solids accounting for anharmonicity is a challenging task, since it requires computing all the quantities appearing in Eq. ( 47) from atomistic models that are sufficiently large to realistically describe disorder and thus computationally expensive.Moreover, we recall that, analogously to Allen-Feldman theory, the Wigner formulation relies on the hypothesis that atoms vibrate around equilibrium positions; thus, it can be applied only to disordered solids and glasses with a structural stability for which this hypothesis is realistic (see e.g.Refs.[90][91][92][93][94][95][96][97] for a discussion of structural stability and related phenomena in various glasses).Because of the multiple challenges highlighted above, the application of the Wigner formulation to glasses will be subject of future work.A summary of the results obtained here relying on the Wigner framework, and their comparison with the standard Peierls-Boltzmann framework is reported in Table I.
It is worth mentioning that Eq. ( 47) -with the populations conductivity in the SMA approximation (49)can be obtained from the Green-Kubo formalism considering the Lorentzian spectral-function limit [78].We conclude by noting that the thermal conductivity formula (47) has been derived employing two approximations that might be improved: (i) We considered a time-independent and spatially uniform temperature gradient; (ii) We employed the standard perturbative treatment of anharmonic effects [19,43,62], which considers only the lowest third-order anharmonic terms and treats them as a perturbation to the harmonic terms (we recall that in this standard perturbative treatment temperature effects are accounted for through the Bose-Einstein distributions entering in the scattering operator (38)).Approximation (i) can potentially be improved following the procedure discussed in Refs.[82,98] to account for time-and space-dependent effects on the conductivity.Improving approximation (ii) requires accounting for anharmonic terms in the heat flux expression (i.e.considering the term σ in Sec.V A and thus refining Eq. ( 42)) as well as for higher-than-third-order phonon scattering processes [20,99] in the collision operator (38).Another related possible improvement concerns accounting for temperature-renormalization effects on the phonon frequencies and collision operator with accuracy beyond perturbation theory, a challenge that can potentially be tackled relying on the stochastic self-consistent harmonic approximation (SSCHA) [100][101][102][103]. Rigorously accounting for all these anharmonic effects in the Wigner formulation requires additional work, to ensure a consistent treatment of the approximations that have to be performed.More details on the capability of the approximated treatment of anharmonicity employed here to accurately describe real solids will be provided later in Sec.VII and Appendix C.

VI. PHASE CONVENTION AND SIZE CONSISTENCY
As anticipated in Sec.II, different phase conventions, smooth or step-like, are employed in the literature for the Fourier transform.Up to now all the derivations have been performed using the smooth phase convention, since in Sec.II we discussed how quantities computed employing such a convention are invariant under different possible choices of the crystal's primitive cell, while quantities computed using the step-like phase convention can be ill-defined (i.e. they can depend on the non-univocal possible choice for the primitive cell).From an intuitive point of view, the choice of using the smooth phase convention in the derivation (via the Taylor expansion (32)) of the LWTE (37), might be motivated noting that the derivative of the square root of the smooth dynamical matrix (7) has variations over its elements that are smoother than those of the derivative of the square root of the step-like dynamical matrix (9).Therefore, the use of the smooth phase convention intuitively yields to a more appropriate description of the smooth variations typical of the close-to-equilibrium regime in focus here.In the following we report quantitative evidence that confirms this intuitive argument.
We start by discussing how the LWTE (37) and related thermal conductivity expression (47) depend on the phase convention adopted in the derivation.
⌦ in (q, q 00 ) s,s 00 h n(R,q 00 ,t) s 00 N(q 00 ) s 00 i d 3 q 00 < l a t e x i t s h a 1 _ b a s e 6 4 = " R C K 9 S Q R 1 / g N g W g r n F N j j o 7 r p

Regime of validity
Vibrations are well defined quasiparticles excitations, ω(q) s > Γ(q) s ∀ q, s (when this condition is not verified, vibrations are overdamped, neither the Wigner nor the Peierls-Boltzmann approach can be applied and spectral-function approaches are needed [78]).
• Simple crystals, where κ αβ P κ αβ C ; • complex crystals, where ∆ω avg ∼Γ(q) s for most phonons (q) s , and Simple crystals characterized by where ω max is the maximum phonon frequency, and 3N at is the number of phonon bands.
As shown by Eqs.(10,11), quantities in Fourier space obtained using the step-like phase convention are related by a unitary transformation to those obtained using the smooth phase convention.In practice, using the step-like phase convention implies removing the atomic positions τ b , τ b from the complex exponentials in the derivations of Secs.III,IV; from Eq. ( 21) to Eq. ( 28) this produces no differences in the physical description.However, the phase convention adopted affects the LWTE (37), since the velocity operator appearing in such an equation depends on the phase convention; it is possible to demonstrate instead that all other coefficients in Eq. ( 37), which are not obtained from a differentiation procedure, are unaffected by the phase convention.Specifically, using the relation (10) between the dynamical matrices in the smooth (D(q) bα,b α , Eq. ( 7)) or step-like (D(q) bα,b α , Eq. ( 9)) phase conventions, it is possible to show that the velocity operator in the step-like phase convention (where E(q) s ,b α are the eigenvectors of the step-like dynamical matrix defined in Eq. ( 11)) is related to that obtained using the smooth phase convention (v β (q) s,s , defined in Eq. ( 36)) by: From Eq. ( 50) it is evident that the elements v β (q) s,s of the velocity operator, either on the diagonal (s=s ) or coupling degenerate vibrational modes (i.e.s =s with ω(q) s =ω(q) s ), do not depend on the phase convention.
In contrast, off-diagonal elements of the velocity operator that couple non-degenerate vibrational modes depend on the phase convention.We note that in any subspace spanned by degenerate vibrational modes with wavevector q and with frequency ω d [105], one can exploit the freedom of rotating the degenerate eigenstates of the dynamical matrix (Eq.( 8)) to render the representation of the velocity operator in direction β diagonal in the mode indexes [11].As mentioned in Sec.IV B, the different Cartesian components of the velocity operator in general do not commute, thus they can not be all diagonalized simultaneously.Nevertheless, since one has the freedom to diagonalize in the degenerate subspace at least one Cartesian component of the velocity operator -thus to recover along such direction a populations' conductivity exactly equivalent to the Peierls-Boltzmann conductivity discussed in Ref. [11] -we consider the diagonal or degenerate velocity-operator elements as contributing exclusively to the Peierls-Boltzmann conductivity, implying that coherences emerge exclusively from off-diagonal and non-degenerate velocity-operator elements.
We investigate the differences between the smooth and step-like phase conventions on the conductivity, focusing on the general property that the conductivity has to respect -namely, size consistency.To this aim, we compute the conductivity (47) of a silicon crystal described in two mathematically different but physically equivalent ways: (i) repeating N ×N ×N times a primitive cell; (ii) Size consistency of the LWTE conductivity in the smooth phase convention.Size consistency is tested computing the conductivity (47) of a silicon crystal at T =300 K and using a fictitious linewidth of Γ(q)s=30 cm −1 ∀ q, s (see text).We sketch in panels a,b) two equivalent crystal representations: a) repeating a 2-atom primitive cell the same number of times along each Cartesian direction; b) repeating a 16-atom unit cell that is larger than the primitive.Panel c) shows that the smooth phase convention (used in Eq. ( 7) yields size-consistent numerical results, since the conductivity is unchanged when the silicon crystal with fictitious linewidths is modeled as: (i) a 16×16×16 repetition of the 2-atom primitive cell; (ii) a (physically equivalent) 2×2×16 repetition of a unit cell that is a supercell 8×8×1 of the primitive.Panel d) shows instead that the steplike phase convention (used in Eq. ( 9)) is not size-consistent, since it yields different conductivities for the two physically equivalent representations of the crystal.The black and red lines serve as a guide to the eye, to show that the Peierls-Boltzmann (populations) conductivity does not depend on the phase convention adopted (black), while the coherences conductivity depends on the phase convention adopted (red).
repeating N 8 × N 8 ×N times a unit cell that is a supercell 8×8×1 of the primitive cell used at point (i).A schematic representation of these two cases is shown in Fig. 2a,b for N = 8, with panel a) representing case (i) and panel b) representing case (ii).Since the differences between the conductivities in smooth and step-like phase conventions can emerge only from the off-diagonal elements of the velocity operator (Eq.( 50)), we perform tests where we use fictitious large linewidths Γ(q) s =30 cm −1 ∀ q, s to enhance the importance the off-diagonal elements of the velocity operator in the conductivity formula (47); therefore, the conductivities obtained in these tests are informative about the size-consistency of Eq. ( 47) in the smooth or step-like phase convention, but their val-ues have no physical meaning.Results show that the primitive-cell and supercell descriptions, which are equivalent from the physical point of view, result having the same thermal conductivity when the smooth phase convention is adopted (Fig. 2c), while they lead to different conductivities when the step-like phase convention is adopted (Fig. 2d).This reiterates that the smooth phase convention (7) has to be used in the derivation of the LWTE (37) (as opposed to the step-like one used in Ref. [37], albeit with negligible numerical differences for that case study), in agreement the intuitive arguments made earlier.

VII. CASE STUDIES: La2Zr2O7 AND CsPbBr3
In this section we consider complex crystals, whose thermal conductivity is very low and not correctly described by the LBTE, as case studies for this more general Wigner framework (Eq.( 47)).In particular, we analyze the La 2 Zr 2 O 7 lanthanum zirconate and the CsPbBr 3 perovskite, as materials used for thermal barrier coatings and thermoelectric devices, respectively.Experimental measurements in these materials [106][107][108][109][110][111] have highlighted a high-temperature trend for the temperaturethermal conductivity relation that has a decay much slower than the T −1 predicted by Peierls-Boltzmann (we recall that the LBTE predicts this T −1 trend to be universal for all crystals [19]).The LWTE conductivity (47) encompasses the T −1 -decaying Peierls-Boltzmann conductivity κ P , but also adds the coherences conductivity κ C , which can increase with T (e.g. when it reduces to the Allen-Feldman conductivity in the limiting case of a harmonic glass [26,112]).These analytical considerations suggest that the LWTE conductivity allows for the emergence of much milder decays or glass-like trends in the thermal conductivity, and further motivates the study of these materials with such more general Wigner framework.

A. Lanthanum zirconate
La 2 Zr 2 O 7 is a material characterized by wide thermal stability and ultralow thermal conductivity; thus, it is widely used for thermal barrier coatings [113].This is a good test case for the Wigner formulation discussed here, since at high temperatures La 2 Zr 2 O 7 falls in the category we defined of complex crystals, with many overlapping phonon bands and a temperature-conductivity curve that is not correctly described by the LBTE.Here, we calculate the thermal conductivity of La 2 Zr 2 O 7 as a function of temperature, computing from first principles all the quantities needed to evaluate Eq. ( 47) (see Appendix F 1 for details).
La 2 Zr 2 O 7 is characterized by a cubic structure (spacegroup: Fd-3m, number 227) and by an isotropic thermal conductivity tensor; in the following only the single in-dependent component of the conductivity tensor will be reported (κ=κ xx =κ yy =κ zz ).Moreover, we rely on the SMA approximation to reduce the computational cost for solving the populations' part of the LWTE (i.e. the LBTE, see Sec.V B), since it is known that for materials with ultralow conductivity the solution of the LBTE determined within the SMA approximation is practically indistinguishable from the exact solution of the LBTE [11,37,42,84].Therefore, in this work the conductivity is computed from Eq. ( 47) with the population term evaluated in the SMA (49).

Thermal conductivity
The conductivity as a function of temperature for La 2 Zr 2 O 7 is shown in Fig. 3.In La 2 Zr 2 O 7 at low temperature, the populations conductivity κ P dominates over the coherences conductivity κ C .
Increasing the temperature, κ P decreases following the universal 1/T decay of the Peierls-Boltzmann equation [19,21,23,116] -a trend that is in broad disagreement with experiments [106][107][108][109][110]. The contribution of the coherences to the conductivity instead increases with temperature, and becomes dominant at high temperature in the complex-crystal regime, where it offsets the incorrect 1/T decay of Peierls-Boltzmann conductivity, leading to a total conductivity (κ TOT =κ P +κ C ) that is in much closer agreement with experiments.The experimental values of the thermal conductivity of La 2 Zr 2 O 7 reported in Fig. 3   show the phonon spectrum of La2Zr2O7 (gray lines) with the phonon linewidths (the full amplitude of the shaded gray areas correspond to the full phonon linewidth Γ(q)s) at 200 K and 1300 K, respectively.Colored circles represent the phonon eigenstates (q)s sampled in the calculation; the area of each circle is related to its contribution to κ (area ∝ √ contribution to κ [115]) and the color identifies the origin of the contribution: green for populations' propagation and blue for coherences' anharmonic couplings (red corresponds to 50% of each, see Eq. ( 54) and related discussion for details).In the coherences' couplings between two modes (q)s and (q) s the contribution of the single mode s is determined as C(q)s/[C(q)s+C(q) s ], where C(q)s is the specific heat of a given phonon (see Appendix E for the details).Panels a(ii), b(ii) show the thermal conductivity density of states κω of populations ("P", green) and coherences ("C", blue) at 200 K and 1300 K, respectively.The gray line shows the vibrational density of states (VDOS).Panels a(iii), b(iii) show the cumulative total thermal conductivity (in black) as a sum of the populations' contribution ("P") in green and coherences' one ("C") in blue, at 200 K and 1300 K, respectively.At 200 K (panel a) La2Zr2O7 is a simple crystal with interband spacings larger than the linewidths and κP kC , while at 1300 K it is a complex crystal with interband spacings smaller than the linewidths and κP < ∼ kC .
mined by the lattice thermal conductivity, with radiative effects having negligible impact on these values.In particular, the data from Yang et al. [110] refer to experimental measurements of the lattice (phonon) contribution to the conductivity of bulk La 2 Zr 2 O 7 .Importantly, Yang et al. [110] discuss how radiative effects are significantly reduced in sintered composites, and all the other experimental data reported in Fig. 3 (Refs.[106][107][108][109]) refer to sintered samples for which radiative effects are suppressed.Finally, given that the smooth phase convention is the only one yielding a size-consistent conductivity, the predictions reported in Fig. 3 have been computed using the smooth phase convention.For completeness, we discuss in Appendix D how the conductivity would change if the step-like phase convention was adopted, showing that the coherences' conductivity at high temperature becomes overestimated, while the populations' conductivity is not affected (this is clear from Eq. ( 50), which shows that the diagonal elements of the velocity operator that enter in the populations' equation (44) are the same in both the smooth and step-like phase conventions).

Transport mechanisms: from simple to complex crystal
In Fig. 4a(i) we show that in La 2 Zr 2 O 7 at low temperature (200 K) the predominance of the populations conductivity over the coherences conductivity is related to linewidths smaller than the phonon interband spacings, demonstrating that at 200 K La 2 Zr 2 O 7 behaves as a simple crystal.More precisely, we show in Figs.4a(ii,iii) that such a dominant populations conductivity is mainly determined by low-frequency phonon bands with large group velocities and weak anharmonicity.Fig. 4b(i) investigates instead the behavior of La 2 Zr 2 O 7 at high temperature, showing that at 1300 K the predominance of the coherences conductivity over the populations conductivity is related to interband spacings smaller than the linewidths, i.e. to a complex-crystal behavior.Figs.4a(ii,iii) show that the dominant coherences conductivity originates from highly anharmonic flat bands.
The density of states of the coherences' conductivity (blue histogram in Fig. 4a,b(ii)) can be also resolved in terms of the frequencies ω 1 , ω 2 of the two modes coupled, as shown in Fig. 5.At 200 K, such a 2-dimensional density of states for the coherences' thermal conductivity (κ C,ω1ω2 , Fig. 5, left panel) shows couplings between quasi-degenerate vibrational modes, similarly to the case of harmonic glasses (Γ(q) s →η→0 ∀ q, s).At 1300 K (Fig. 5, right panel) κ C,ω1ω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.We conclude by noting that the results presented in this section were obtained employing the approximated treatment of anharmonicity discussed in Sec.V B. We show in Appendix C that the frequencies computed with this approximated scheme yield predictions for the specific heat in agreement with experiments in the temperature range analyzed.Moreover, we also show that simulating the Raman spectra using the frequencies and linewidths obtained from this approximated scheme yields predictions in agreement with experiments over a temperature range in which coherences are not negligible.This shows that the standard perturbative treatment for anharmonicity adopted is accurate enough for the illustrative scope of this analysis.

B. Halide perovskite CsPbBr3
The perovskite CsPbBr 3 belongs to a family of ultralow-thermal-conductivity materials that are promising for thermoelectric energy conversion [22,111], and it has been used in Ref. [37] to showcase the predictive capabilities of the Wigner formulation.The temperature-conductivity curve for CsPbBr 3 reported in our earlier work [37] was computed using the step-like phase convention, which has been shown in Sec.VI to yield size-inconsistent results.Therefore, in this section we update the LWTE predictions for the temperatureconductivity curve of CsPbBr 3 using the correct sizeconsistent smooth phase convention discussed before.We show in Fig. 6 a comparison between the conductivity predicted with the LWTE conductivity formula (47) in the smooth phase convention and the experimental measurements [111], finding results which are not significantly different from those reported in Ref. [37] (see Appendix D for a detailed analysis on how the x component of the conductivity tensor of CsPbBr 3 changes between the size-consistent smooth phase convention and the size-inconsistent step-like phase convention).Specifically, Fig. 6 shows that at 300 K the populations con- FIG. 6. Bulk thermal conductivity of CsPbBr3 in the smooth phase convention."Exp(L)", "Exp(M)" and "Exp(S)" refer to measurements [111] of κ xx in nanowires having respectively sections of 800×380 nm 2 , 320×390 nm 2 and 300×160 nm 2 : their broad agreement supports the hypothesis of negligible finite-size boundary scattering [111].Lines are theoretical predictions of the conductivity from equation ( 47) using the smooth phase convention, with solid lines referring to the tensor component κ xx (to be compared with experiments), and dotted lines referring to the component κ yy (this last is reported here for completeness and is practically indistinguishable from the not-reported component κ zz ).Green, Peierls' LBTE conductivity (κP), which displays the universal T −1 asymptotics (dashed line) [19,21,23,116].Blue, coherences' thermal conductivity (κC).Black, total conductivity from equation ( 47): κTOT = κP + κC.
ductivity contributes just ∼ 1 3 of the total conductivity, while the coherences' term provides an additional ∼ 2  3 , leading to a much closer agreement with experiments that becomes even more relevant in the high-temperature asymptotics.Conversely, at low temperature the populations' conductivity becomes dominant, and at 50 K it provides ∼ 3 4 of the total conductivity.It can be seen that in the high-temperature limit κ xx P ∝ T −1 , as predicted by Peierls' theory [19,21,23,116]; this in disagreement with experiments.Instead, the LWTE formulation predicts a decay of κ xx much milder than T −1 , as shown here for CsPbBr 3 , and as present in many other complex crystals [21, 23, 33-35, 117, 118].Finally, in Appendix C we check the accuracy of the standard approximated treatment of anharmonicity we employed, showing that also for CsPbBr 3 the frequencies and linewidths we computed yield predictions for the specific heat and Raman spectra in agreement with experiments.

VIII. PARTICLE-WAVE DUALITY FOR HEAT TRANSPORT
We recall that the central object of the Wigner formulation is a matrix distribution n(R, q, t) s,s , where position (R), wavevector (q) and time (t) are the arguments, and the two matrix indexes are the phonon band labels s, s .Importantly, in the out-of-equilibrium steady-state regime considered, the Wigner matrix distribution (35) is not diagonal in the basis which diagonalizes the leading phonon Hamiltonian (8) (from Eq. ( 37) it is clear that at fixed R, q, t, the matrices n(R, q, t) s,s , ω(q) s δ s,s , and v (q) s,s , do not commute in general), and the evolution of the diagonal elements (s=s ) of the Wigner distribution is decoupled from the evolution of the off-diagonal elements (s =s ), see Eq. ( 44) and Eq. ( 45).As discussed by Rossi [62], the semiclassical limit corresponds to neglecting the off-diagonal Wigner distribution elements (s =s ), thus considering only the diagonal Wigner-distribution elements n(R, q, t) s,s and interpreting these as number of phonon wavepackets that behave as if they were particles of a classical gas, since these have well defined energies ω(q) s , group velocities v (q) s,s = ∂ω(q)s ∂q , and lifetimes τ (q) s =1/Γ(q) s [119].As mentioned before, such a semiclassical approximation turns out to be accurate in simple crystals characterized by well-separated phonon bands (i.e.interbranch spacings much larger than the linewidths), where the exact solution of the LBTE yields a populations' thermal conductivity that is order of magnitude larger than the coherences conductivity [14] and in good agreement with that measured in experiments [7,9,[15][16][17].In contrast, when off-diagonal coherences s =s cannot be neglected, such a particle-like (semiclassical) picture is no longer valid, since off-diagonal elements are related to the wave-like coherences between quantum vibrational eigenstates (phonons) [62].Specifically, the presence of the difference between frequencies ω(q) s −ω(q) s in the coherences evolution equation (45) demonstrates that: (i) coherences between phonons are possible thanks to their wave-like nature and capability to interfere (such interference is possible also in the presence of damping); (ii) coherences do not have an absolute energy (frequency) akin to that of a particle-like phonon wave-packet.In addition, when there is wave-like coherence between two phonon eigenstates, heat is transfered between these via an interband (Zener-like) tunneling mechanism (see Eq. ( 47)); such a wave-like tunneling between phonon bands is intrinsically different from the particle-like propagation of a phonon wavepacket, since in the latter the "identity" of a phonon (defined by the wavevector and index of the band to which a phonon belongs) is preserved while in the former is not.In the following we derive a quantitative criterion that allows to assess when the particle-like (semiclassical) picture breaks down and wave-like coherences have to be considered.

A. Wigner and Ioffe-Regel limits in time
We start by analyzing the strength of particle-like and wave-like conduction mechanisms, focusing on the conditions under which they become comparable.We start by resolving the contributions of a phonon (q) s to the particle-like and to the wave-like conductivities, recasting the normalized trace of the conductivity tensor (47) in the SMA approximation as [120] where K avg P (q) s and K avg C (q) s quantify the contributions of the phonon (q) s to the particle-like and to the wave-like (coherences) conductivities of Eq. ( 47), respectively (see Appendix E for details).The ratio between these wave-like and particle-like conductivity contributions quantifies the relative strength of the corresponding microscopic conduction mechanisms.To discuss the relative strength of particle-like and wave-like conduction mechanisms, it is useful to introduce the average phonon interband spacing, which is defined as the ratio between the maximum phonon frequency (ω max ) and the number of phonon bands (3N at ).In fact, the ratio between the wave-like and particle-like contributions is approximately equivalent to the ratio between the phonon linewidth and the average phonon interband spacing, The approximated equality in Eq. ( 53) is demonstrated with order-of-magnitude accuracy in Appendix E, relying on the assumption that at fixed s and for variable s , the ratio between the off-diagonal and diagonal velocityoperator elements has an order of magnitude equal to one, 1.Such an assumption can be understood and justified by noting that in Eq. ( 47) the trend of the velocity-operator elements α |v α (q) s,s | 2 as a function of the frequency difference ω(q) s −ω(q) s determines how the coherences conductivity behaves at high temperature.In particular, an approximatively saturating trend for the temperature-coherences conductivity curve -which appears in complex crystals with glasslike conductivity such as La 2 Zr 2 O 7 (Fig. 3) and CsPbBr 3 (Fig. 6) -is related to off-diagonal velocity operator elements |v α (q) s,s | 2 that for fixed s and variable s are approximately independent from s (i.e. they are independent from the frequency difference ω(q) s −ω(q) s ).
We highlight how Eq. ( 53) predicts the crossover from particle-like to wave-like conduction for phonons to be not sharp, with the possibility to have phonons that contribute simultaneously and comparably to populations' particle-like conductivity and to coherences' wave-like conductivity.Importantly, phonons with a lifetime τ (q) s equal to the inverse of the average interband spacing [∆ω avg ] −1 are predicted to be at the center of such nonsharp crossover.Since such prediction is obtained via the Wigner framework, hereafter we will refer to the time scale [∆ω avg ] −1 as "Wigner limit in time".
It is worth mentioning that another condition on the phonon lifetime is discussed in the literature, in relation to the validity of theoretical frameworks that describe transport relying on the concept of phonon as quasiparticle excitation.Specifically, such condition requires phonons to be above the "Ioffe-Regel limit in time", i.e. to have a lifetime longer than their reciprocal angular frequency, τ (q) s =[Γ(q) s ] −1 > 1 ω(q)s [41,121], implying that the atomic vibrations are not overdamped [122].We note that also this limit is not sharp, and definitions of this limit differing by a factor of 2π can be found in the literature [123].The definition employed here follows the book of Landau et al. [77], where the comparison between the quasiparticle energy ( ω(q) s ) and its "degree of broadening" ( Γ(q) s = [τ (q) s ] −1 ) is used to assess if the concept of quasiparticle excitations is well defined.Therefore, having phonons with a lifetime longer the Ioffe-Regel limit in time is a necessary condition for the validity of the Wigner formulation discussed here, since the particlelike and wave-like conduction mechanisms discussed here are derived relying on the assumption that phonons are well defined quasiparticle excitations (i.e. one can label them with a wavevector q and with a mode index s).We also note that having phonons with a lifetime Regime diagram for thermal transport.Three regimes can be distinguished in a diagram reporting the phonon lifetime (τ ) as a function of the phonon frequency (ω).The horizontal black like is the Wigner limit in time, τ = 1/∆ωavg (i.e. the inverse average interband spacing, see Eq. ( 52), and denotes the center of the non-sharp particle-wave crossover for phonons.The green region well above the line represents phonons that propagate particlelike and mainly contribute to the populations conductivity; the blue region below the line represents phonons that tunnel wave-like and mainly contribute to the coherences conductivity; the red points around the line are phonons that contribute simultaneously and comparably to the particle-like and wave-like conductivities.The purple line τ =1/ω is the center of the non-sharp Ioffe-Regel limit in time [41,121].The region above it represents well defined phonons; the Wigner transport equation ( 37) describes all these phonons (blue, red, green), while the semiclassical Peierls-Boltzmann equation accounts only for phonons that propagate particlelike (green).The yellow region below the Ioffe-Regel limit represents overdamped phonons, which require full spectralfunction approaches [78,79] to be described correctly.
longer than the Ioffe-Regel limit in time is a necessary but not sufficient condition for the validity of the semiclassical particle-like description, since having well defined phonons does not imply that wave-like conduction mechanisms are negligible.When, instead, the phonon lifetimes become shorter than the Ioffe-Regel limit in time (τ (q) s < 1 ω(q)s ), i.e. in the so-called overdamped regime, phonons are no longer well defined quasiparticle excitations.This implies that it is no longer possible to describe scattering in terms of the phonon wavevector q and mode s, thus Eq. ( 38) is no longer valid and one has to address this regime considering phonon spectral functions, as discussed in Refs.[78,79].All these regimes are summarized in Fig. 7.
With the goal of validating the analytical considerations above (Eq.( 53)), we plot in Fig. 8 the phonon lifetime τ (q) s as a function of the phonon energy ω(q) s and at different temperatures for La 2 Zr 2 O 7 (panels a, b, c) and CsPbBr 3 (panels d, e, f ).The particle-like and wavelike conductivity contributions of Eq. ( 51) (K avg P (q) s and K avg C (q) s , respectively) are used in Fig. 8 to quantify the magnitude and type of the contribution to the conductivity: the area of each circle is proportional to the total contribution to the conductivity, K avg P (q) s +K avg C (q) s , and the color of each circle is set according to the value of c= Specifically, green corresponds to c=1 (dominant particle-like conductivity contribution), blue corresponds to c=−1 (dominant wave-like conductivity contribution), and intermediate colors represent coexisting particle-and wave-like conduction (e.g.red corresponds to 50% of each, c=0).Results in Fig. 8 are in broad agreement with the predictions of Eq. ( 53), with phonons above the Wigner limit in time (i.e. with long lifetime) contributing mainly to the populations' particle-like conductivity, and phonons below the Wigner limit in time (i.e. with extremely short lifetime) contributing mainly to the coherences' wave-like conductivity.These results can be rationalized intuitively: phonons with an extremely short lifetime give a negligible contribution to the particle-like propagation mechanisms because they are suppressed too quickly (i.e.before having the time to propagate significantly), but they can give a significant contribution to the wave-like conductivity since they are still allowed to interfere and tunnel (here "interference" is used in the sense of the frequency difference appearing in the coherences equation ( 45)).Finally, we note that in both La 2 Zr 2 O 7 and CsPbBr 3 all the phonons are above the Ioffe-Regel limit in time (i.e.not overdamped, for orthorhombic CsPbBr 3 below 350 K this is in broad agreement with recent experiments [124]).This confirms that phonons in these materials are well defined quasiparticle excitations, validating a key hypothesis employed in the derivation of the LWTE (37).
FIG. 8. Phonon lifetimes in La2Zr2O7 and CsPbBr3.We show the phonon lifetimes τ (q)=[Γ(q)s] −1 as a function of the energy ω(q)s for La2Zr2O7 (at a temperature of 200 K (a), 800 K (b), and 1300 K (c)), and for CsPbBr3 (at a temperature of 50 K (d), 225 K (e), and 350 K (f )).The area of each circle is proportional to the contribution to the total conductivity and colored according to the origin of the contribution: green for particle-like propagation of populations; blue for wave-like tunneling of coherences; intermediate colors represents phonons contributing to both mechanisms, with red corresponding to 50% of each (see Eq. ( 54)).The dashed-black line is the Wigner limit in time, corresponding to a phonon lifetime equal to the inverse of the average interband spacing (τ = [∆ωavg] −1 see Eq. ( 52)).We highlight how the non-sharp crossover from dominant particle-like conduction to dominant wave-like conduction occurs around this value, in agreement with the theoretical predictions from the Wigner framework (Eq.( 53)).The dashed-purple hyperbola is the Ioffe-Regel limit in time, τ = 1/ω, the Wigner formulation is applicable when phonons have lifetime longer than this value, a condition which is always verified here.The pie charts have an area proportional to the total conductivity, and the slices resolve the populations conductivity (particle-like, green) and the coherences conductivity (wave-like, blue).

B. Ioffe-Regel limit in space
The analysis of the previous section has been performed in the frequency/lifetime domain because frequencies and lifetimes appear explicitly in the LWTE (37) and related conductivity formula (47).Such an analysis can be recast also in the the frequency/mean free path domain, a representation that is commonly employed in the literature [24,27,28,35,42,125] and thus worth to be discussed also here.
We plot in Fig. 9 the phonon mean free path (MFP, Λ(q) s =v avg (q) s,s τ (q) s , where is the spatiallyaveraged modulus of the group velocity, more details are provided in Eq. (E1)) as a function of the phonon energy ω(q) s at different temperatures for La 2 Zr 2 O 7 (panels a, b, c) and CsPbBr 3 (panels d, e, f ).The non-sharp crossover from dominant particle-like propagation mech-anisms to dominant wave-like tunneling mechanisms is evident also in the frequency/MFP domain, and its center is approximatively identified by having MFPs equal to the typical interatomic spacing a, i.e. by phonons which are at the so-called "Ioffe-Regel limit in space" [27,28].More precisely, phonons with Λ(q) s > a (above the Ioffe-Regel limit in space) contribute mainly to the populations' particle-like conductivity, while phonons with Λ(q) s < a (below the Ioffe-Regel limit in space) contribute mainly to the coherences wave-like conductivity.
The appearance of the interatomic spacing as length scale determining the center of the particle-wave crossover can be rationalized with an order-of-magnitude analysis (similar to that leading to Eq. ( 53)), which yields FIG. 9.
Phonon mean free paths in La2Zr2O7 and CsPbBr3.We plot the phonon mean free paths (MFP) Λ(q)s=v avg (q)s,sτ (q)s as a function of the energy ω(q)s for La2Zr2O7 (at a temperature of 200 K (a), 800 K (b), and 1300 K (c)), and for CsPbBr3 (at a temperature of 50 K (d), 225 K (e), and 350 K (f )).The area of each circle is proportional to the contribution to the total conductivity and colored according to the origin of the contribution: green for particle-like propagation of populations; blue for wave-like tunneling of coherences; intermediate colors represents phonons contributing to both mechanisms, with red corresponding to 50% of each (see Eq. ( 54)).The dashed-black line in panels a,b,c (d,e,f ) is the average bond length in La2Zr2O7 (CsPbBr3).We highlight how phonons at the center of the non-sharp crossover from dominant particle-like conduction to dominant wave-like conduction have a mean free path approximately equal to the average bond length (i.e. they are at the Ioffe-Regel limit in space [28,41]), a behavior which is explained by Eq. ( 55).The pie charts have an area proportional to the total conductivity, and the slices resolve the populations conductivity (particle-like, green) and the coherences conductivity (wave-like, blue).
Eq. ( 55) follows from Eq. ( 53), recasting the lifetime in terms of the mean free path, τ (q) s =Λ(q) s /v avg (q) s,s , and using the following estimate for the group velocity: (56) This first approximated equality in Eq. ( 56) derives from the following considerations: (i) each phonon band spans throughout half the Brillouin zone a frequency interval that is approximatively equal to the average interband spacing ∆ω avg (we consider half the Brillouin zone because of the time reversal symmetry ω(q) s = ω(−q) s ; for example, in a representative one-dimensional case, all the variations taking place in the positive half of the Brillouin zone are mirrored in the negative part); (ii) the typical magnitude of a reciprocal lattice vector can be estimated as 2π a3 √ Nat (a is the typical interatomic spacing, N at the number of atoms in the primitive cell, and thus a 3 N at is an estimate of the primitive-cell volume and a 3 √ N at is an estimate of a direct lattice vectors' length).These considerations allow to estimate the typical group velocity as ratio between the typical frequency interval spanned by a phonon band over half the Brillouin zone and half the average length of the reciprocal lattice vector.We note that in the last approximated equality in Eq. ( 56) we have approximated the quantity a π 3 √ Nat a. Within the order-of-magnitude analysis performed here, this means that multiplying the typical bond length a by the factor π 3 √ Nat yields a quantity that can still be considered as typical bond length.Such an approximation is realistic for typical complex crystals, since these usually have a primitive cell containing tens of atoms and for N at ranging from 10 to 100 the ratio π has some freedom -compatible with a factor of π 3 √ Nat -in the definition of the typical interatomic spacing.This freedom justifies the last approximation performed in Eq. ( 56).
In summary, we have shown the existence of a nonsharp crossover from the regime where phonons mainly propagate particle-like to the regime where phonons mainly tunnel wave-like.We have shown that the Wigner limit in time (corresponding to phonons with lifetime equal to the inverse average interband spacing, τ (q) s = [∆ω avg ] −1 ), or equivalently the Ioffe-Regel limit in space (corresponding to phonons with MFP equal to the typical interatomic spacing, Λ(q) s = a), can be used to identify the center of the non-sharp particle-wave crossover.The analysis and data reported in this section demonstrate that the Winger formulation allows to describe thermal transport beyond the Ioffe-Regel limit in space, solving rigorously a long-standing problem that has hitherto been tackled with phenomenological models [35,42].The results reported in Figs.8,9, and their explanation with Eqs.(53,55), shed light on the strengths and weaknesses of the phenomenological heat-conduction models that conjecture the existence of two heat-conduction mechanisms sharply divided by the Ioffe-Regel limit in space [35,42].Specifically, we have shown that: (i) phonons mainly propagate particle-like (tunnel wave-like) above (below) the Ioffe-Regel limit in space, and this supports where the separation between the particle-like and wavelike channels is sharply performed in the phenomenological models of Refs.[35,42]; (ii) the particle-wave crossover for phonons is not sharp, with phonons contributing simultaneously to both populations' particlelike conductivity and coherences' wave-like conductivity, and this is in contrast with the phenomenological sharp separation between these two mechanisms done in Refs.[35,42].

IX. SOFTWARE IMPLEMENTATIONS
Here we report the numerical protocols to implement the general LWTE thermal conductivity expression (47) in commonly used computer programs.Modern LBTE solvers contain almost all the quantities needed to evaluate the general LWTE thermal conductivity; the only missing quantity is the off-diagonal part of the velocity operator, and this can be obtained with minimal effort.
The starting point to compute the velocity operator (36) is the mass-renormalized interatomic-forceconstant tensor G Rbα,R b α (Eq.( 4)).Once the tensor G Rbα,R b α is known, one has to compute its Fourier transform in the smooth phase convention using Eq. ( 7), to obtain the smooth dynamical matrix in reciprocal Bloch representation D(q) bα,b α .We remark that some LBTE solvers (e.g.Phono3py [6]) employ the smooth phase convention as default, others instead employ the step-like phase convention (e.g.D3Q [11,76]).As discussed before, it is fundamental to use the smooth phase convention to implement correctly the generalized LWTE conductivity (47), and we refer the reader to Sec.II for details on how to transform quantities from the steplike phase convention to the smooth phase convention needed here.Then, one needs to compute the square root of the smooth dynamical matrix in reciprocal Bloch representation, D(q) bα,b α .To compute this quantity, we rely on the property that the square-root matrix D(q) bα,b α has the same eigenvectors E(q) s,bα of the matrix D(q) bα,b α , and the eigenvalues of D(q) bα,b α are the square roots of the eigenvalues of D(q) bα,b α .Eigenvectors and eigenvalues of D(q) bα,b α have been introduced in Eq. ( 8) (and are available in modern LBTE solvers), and allow to compute the square root of the smooth dynamical matrix in reciprocal Bloch representation as Once the square root of the dynamical matrix is known, its derivative -needed to evaluate the velocity operator (36) -is computed as: and the limit here can be evaluated numerically employing standard finite-difference methods.Once Eq. ( 58) is evaluated numerically, inserting its result into Eq.( 36) allows to evaluate numerically the velocity operator.
After having computed the velocity operator, one can use the vibrational frequencies, phonon specific heats, and phonon linewidths -which are available in modern LBTE solvers [3,4,6,10,11,44,76] -to evaluate numerically the coherence conductivity (κ αβ C ) appearing in generalized LWTE conductivity formula (47).The total conductivity is finally obtained adding the Peierls-Boltzmann conductivity (κ αβ P , already computed by LBTE solvers) and the coherences conductivity, Finally, we stress that the total LWTE conductivity κ αβ TOT is the physically relevant quantity.To have a Peierls-Boltzmann particle-like conductivity κ αβ P equivalent to that discussed in previous work [11], one has to consider the diagonal or degenerate velocity-operator elements (i.e.v β (q) s,s such that ω(q) s = ω(q) s ) as contributing exclusively to the Peierls-Boltzmann conductivity.A possible way to account for this prescription consists in exploiting the freedom of rotating the degenerate eigenstates of the dynamical matrix in reciprocal Bloch representation to obtain a velocity operator that is diagonal in the degenerate subspaces, as done in Ref. [11].

X. CONCLUSIONS
In summary, we have derived from the Wigner phase space formulation of quantum mechanics a Wigner transport equation (37) able to account rigorously and simultaneously for the interplay between the quantum Bose-Einstein statistics of atomic vibrations, anharmonicity, and disorder in determining the thermal properties of solids.
We have discussed how the Wigner phase-space framework determines an expression for the microscopic harmonic energy field and the related microscopic harmonic heat flux; such an expression allows to compute the thermal conductivity from the solution of the linearized Wigner transport equation (LWTE).At leading order in the temperature gradient, the solution of the LWTE (37) yields a generalized LWTE thermal conductivity (47) that is the sum of two terms: a populations conductivity that describes the particle-like intraband propagation of phonons exactly as in the Peierls-Boltzmann LBTE, and a coherences conductivity describing the wave-like interband tunnelling of phonons.We note in passing that the coherences thermal conductivity has some analogies with the Zener interband electrical conductivity [85,88,89,127,128]; however, thermal transport is fundamentally different from electronic transport since the full phonon spectrum contributes to the heat current, while only the bands close to the Fermi level contribute to the electronic Zener current.The LWTE conductivity (47) reduces to the LBTE conductivity in the limit of a simple crystal [14] (i.e.characterized by interband spacings much larger than the linewidths) and to the Allen-Feldman conductivity in the limit of a harmonic glass; most importantly, the LWTE conductivity ( 47) is more general and encompasses all intermediate cases, applying, for example, also to complex crystals characterized by interband spacings smaller than the linewidths.Moreover, we have discussed how the Wigner formulation extends the Allen-Feldman theory accounting for the interplay between disorder, anharmonicity, and quantum Bose-Einstein statistics of atomic vibrations; however, evaluating the Wigner formulation in disordered solids is a computationally challenging task.Fittingly, recent work [129] has discussed techniques that allow to account for anharmonicity in disordered systems at a reduced computational cost, and future work will aim at understanding if these techniques can be combined with the Wigner formulation to better understand the anomalous thermal properties of glasses [130][131][132][133][134][135][136][137][138].
We have discussed how the coefficients appearing in the LWTE (37), and in the LWTE conductivity (47) that follows from it, depend on the phase convention adopted for the Fourier transform.Specifically, we have discussed the necessity to compute the dynamical matrix in reciprocal Bloch representation and related off-diagonal velocity operator's elements using the smooth phase convention (7), employed e.g. by Wallace [40], in order to obtain a thermal conductivity that does not depend on the non-univocal possible choice of the crystal's unit cell and is size-consistent.Therefore, the formulation presented here improves and corrects our previous work [37], where the commonly employed [19,43] step-like convention (9) was adopted.We have also shown that these improvements/corrections become increasingly more relevant as the size and anisotropy of the unit cell used to describe a crystal increase.
We have shown that the present approach overcomes the limitations of the Peierls-Boltzmann formulation, predicting accurately the temperature-conductivity curves of complex crystals such as La 2 Zr 2 O 7 and CsPbBr 3 , both displaying a decay much milder than the universal T −1 trend predicted by Peierls' theory (Figs.3,6).We have shown that the LWTE formulation predicts a non-sharp crossover from dominant particle-like heat-conduction mechanisms to dominant wave-like heatconduction mechanisms, which is centered at the Wigner limit in time (corresponding to τ (q) s ∼[∆ω avg ] −1 , where τ (q) s is the phonon lifetime and ∆ω avg is the average interband spacing).We have discussed the equivalence between the Wigner limit in time and the Ioffe-Regel limit in space, this latter corresponding to a phonon mean free path (Λ(q) s ) equal to the the typical interatomic spacing (a), Λ(q) s a. Specifically we have shown that phonons with τ (q) s [∆ω avg ] −1 (or equivalently with Λ(q) s a) mainly propagate particle-like and contribute to the populations conductivity; phonons with τ (q) s [∆ω avg ] −1 (or equivalently with Λ(q) s a) mainly tunnel wave-like and contribute to the coherences conductivity; phonons with τ (q) s [∆ω avg ] −1 (or equivalently with Λ(q) s a) behave simultaneously particle-and wave-like and contribute comparably to the populations and coherences conductivities.We have exploited these findings to propose a classification, based on the heat conduction mechanisms, of the various regimes of thermal transport, also identifying the regime of applicability of the semiclassical Peierls-Boltzmann equation (which accounts for particlelike conduction only) and of the Wigner equation (37) (which accounts for both particle-and wave-like conduction).We have provided numerical recipes to implement with minimal effort the general LWTE thermal conductivity formula in LBTE solvers [4,6,10,11,44].
It is also worth mentioning that other thermal conductivity expressions valid in both crystals and glasses have been proposed relying on the Green-Kubo theory of linear response [79,[139][140][141], which differ under some aspects from Wigner expression (47) discussed here.On the one hand, the Green-Kubo formulations reported to date are restricted to the so-called "kinetic regime" of thermal transport, where Umklapp scattering processes dominate and the single-mode relaxation-time approximation is accurate [11,84,142].Such a restriction is not present in the Wigner transport equation (37), which can be used also in the hydrodynamic regime of thermal transport (i.e. when normal scattering processes dominate and the single-mode relaxation-time approximation breaks down [11,14,142]).On the other hand, the Green-Kubo framework allows to describe transport in the overdamped regime, i.e. below the Ioffe-Regel limit in time, where phonons are no longer well defined quasiparticle excitations and they have to be described using spectral functions [78].We highlight that the complex crystals studied in this work are not in the overdamped regime (i.e. they have phonons above the Ioffe-Regel limit in time), but conditions under which the overdamped regime can be reached have been discussed recently [124].The approaches discussed up to now account for anharmonicity approximatively, truncating at third or higher order the expansion of the interatomic potential around equilibrium (Eq.( 3)).We note that approaches based on molecular dynamics (MD), either in the equilibrium [143][144][145][146][147][148][149], non-equilibrium [112,[150][151][152][153], or approach-to-equilibrium [154,155] flavors, allow to account for anharmonicity exactly, although are computationally more expensive than the Wigner approach and do not account for the quantum Bose-Einstein statistics of atomic vibrations [156], which affects significantly the thermal properties around and below the Debye temperature.

√
G −1 R b α ,Rbα , see main text).Therefore, Eq. (A2) shows that the bosonic operator â † (R) bα creates atomic vibrations in which atoms close to R+τ b have a larger mean square momentum than atoms far from R+τ b .

Localization problems of the phonon-mode basis
In this Appendix we show that the standard phonon operators â(q) s defined in Eq. ( 26), have in general nonunique real-space representation and as such are not suitable to track vibrations in real space.
As anticipated in Sec.III, the phase of the eigenvector E(q) s,bα appearing in Eq. ( 26) is undetermined at all qpoints; i.e., there exists a gauge freedom in the definition of the eigenvectors E(q) s,bα that allows to apply them unitary transformations where U(q) s,s is a unitary matrix of size 3N at ×3N at periodic in q and with non-zero matrix elements exclusively within a subspace of phonon bands that is invariant with respect to the space group of the crystal [176].The freedom inherent in Eq. (A3) implies that the real-space (Wannier) representation of the phonon operator ( 26) is non-unique and its localization properties depend on the particular choice made for the arbitrary unitary matrix U(q) s,s .To illustrate this, we consider the case in which U(q) s,s =e −iq•R U δ s,s , where R U is an arbitrary Bravais-lattice vector, and we use the notation âU (q) s = â(q) s e iq•R U to emphasize that the phonon operator (26) depends on the arbitrary unitary transformation U(q) s,s (i.e. on the choice of R U in this example).The real-space (Wannier) representation of the phonon operator ( 26) is: and clearly âU (R) s inherits from âU (q) s a dependence from the aforementioned arbitrary unitary transformation.Now we investigate the localization properties of âU (R) s performing an analysis analogous to that dis-cussed in Sec.III.We create on the ground state |0 the excitation â † U (R) s |0 , and we compute the mean-square atomic displacement at a generic position R +τ b in such an excited state, finding where 0| û2 (R ) b α |0 is the mean-square atomic displacement due to the zero-point motion discussed in Sec.III.From Eq. (A5) it is apparent that the excitation created by the operator â † U (R) s on the ground state |0 is centered at a Bravais lattice position R −R−R U , which can be arbitrarily shifted since the Bravais-lattice vector R U appearing in the unitary transformation e iq•R U can be chosen arbitrarily.
In conclusion, we have shown that the standard phonon bosonic operators (26) have non-unique real-space representations, each of which creates atomic vibrations with different center.This mirrors the electronic case, where the phase indeterminacy of the Bloch orbitals is reflected in the non-uniqueness of the transformation into Wannier functions [176].In contrast, the localized bosonic operators in real space (14) have by construction well-defined localization properties (see Eq. ( 17) and Fig. 1), which allow to track vibrations and thus perform the derivation detailed in Sec.IV.

Appendix B: Properties of the Wigner matrix distribution
Here we discuss the properties of the Wigner matrix defined in Eq. (29).Specifically, the Wigner representation of a one-body operator Ô is a distribution that depends on a direct lattice vector R and on a wavevector q, it can depend on time t, and carries two matrix indexes bα, b α that denote positions inside the primitive cell of the crystal; in this Appendix, such a distribution will be denoted with W [O] (R, q) bα,b α =O(R, q) bα,b α (the time dependence will be reported only if needed).Despite this phase-space distribution is discussed for the particular case of vibrational excitations (phonons) here, in principle it can be employed also to describe other types of quasiparticles in a periodic potential.Specifically, the one-body matrix element O(q+ q 2 , q− q 2 ) bα,b α appearing in Eq. ( 29) is a matrix element in the Zak basis [177][178][179], which has as quantum numbers eigenvalues of translation operators in direct and reciprocal spaces (that is a position inside a unit cell τ , which is a continuous position for electrons and the discrete atomic position for phonons, and a wavevector q belonging to a Brillouin zone).We also note that, in order to adapt the formulation discussed here to other quasiparticles, the phase convention for the Zak eigenfunctions has to be chosen according to the convention (21), i.e. it must yield a transformation that relates the matrix elements in the Zak basis to the matrix elements in the standard position basis analogous to Eq. (25).
We want now to discuss and demonstrate the mathematical properties satisfied by the direct and inverse Wigner transformations (Eq.( 29) and Eq. ( 30)).
a. Proof of the inverse Wigner transform (30).To prove that the transformation (30) is the inverse of the transformation (29), it is sufficient to show that inserting the first into the second (or vice versa) one obtains the identity: In Eq. (B1) we have used the exponential representation of the Dirac delta [58,59] where N c → ∞ is the number of Bravais lattice vectors over which the sum R runs, G is a reciprocal lattice vector, δ q −q ,G is the Kronecker delta, and δ(q −q −G) is the Dirac delta.In the last line of Eq. (B2) we have exploited the property that only the reciprocal lattice vector G = 0 has to be considered when both q and q are restricted to the first Brillouin zone B [59,180] (we stress that this is the case in Eq. (B1)).
b.The Wigner transform of a translation-invariant operator does not depend on space.Eq. ( 29) implies that the Wigner transform of a translation-invariant operator, i.e. diagonal in reciprocal space, is independent from the spatial position r.To show this, we start by considering a translation-invariant one-body operator O(R, R ) bα,b α =O(R−R ) bα,b α .Performing the Fourier transform (25), and defining R =R−R , yields where we have used the exponential representation of the Dirac delta (B2) in the case where the wavevector q ∈ B. It is now evident that inserting the translationinvariant operator (B3) in the Wigner transform (29) yields a distribution that does not depend on space.This property has been exploited to transform Eq. ( 40) into Eq.( 41).Specifically, following the procedure above one can verify that the Wigner representation of the square root force-constant matrix √ G Rbα,R b α -which is translation-invariant, as a consequence of Eq. ( 6)does not depend on space and is equal to the square root of the smooth dynamical matrix in reciprocal Bloch representation D(q) bα,b α .c.The Wigner matrix distribution is Hermitian.From the hermiticity of a generic observable, O(q, q ) bα,b α =O * (q , q) b α ,bα , it is straightforward to prove that the Wigner matrix distributions defined by Eq. ( 29) are Hermitian, i.e. they satisfy d. Marginal of the Wigner matrix distribution.Integrating the Wigner matrix distribution over all the Bravais-lattice vectors yields the diagonal elements of the one-body density matrix in reciprocal space: where we have used the exponential representation of the Dirac delta (B2) in the case where the wavevector q ∈B.e. Inner product.As discussed in Eq. ( 31) of the main text, the trace of the product of two one-body operators in the Dirac framework can be rewritten as phasespace integral involving the Wigner representations of these operators.Here we prove such a property.We start using the properties of the Dirac delta to recast the first line of Eq. ( 31) as follows ×A q− q 2 ,q+ q 2 b α ,bα d 3 qd 3 q d 3 q . (B5) Since in Eq. (B5) both q and q belong to the first Brillouin zone B, we can rewrite the Dirac delta using the exponential representation (B2), obtaining: q+ q 2 b α ,bα d 3 qd 3 q d 3 q . (B6) Now we exchange the order of sum and integrals, thus we can recognize that the integrals over q and q yield the phase-space representation of the one-body density matrix and of the one-body operator A(q, q ) bα,b α , i.e.Eq. (31).
Appendix C: Accuracy of the perturbative treatment of anharmonicity in La2Zr2O7 and CsPbBr3 As discussed in Sec.V B, the standard perturbative treatment of anharmonicity employed in this work considers the vibrational frequencies as temperature independent (bare phonon frequencies), and accounts for the effect on temperature via the Bose-Einstein distributions appearing in the scattering operator (38).
To check how accurate it is to consider the bare phonon frequencies (and related atomic displacement patterns defined in Eq. ( 8)) as a representation of the actual vibrational frequencies (and related atomic displacements) over the entire temperature range analyzed, we use them to compute the specific heat and compare our calculations with experiments.Specifically, the specific heat at constant volume (since thermal expansion is neglected in the standard perturbative treatment we adopted) is computed as 2 ω 2 (q)s k B T 2 N(q) s N(q) s +1 d 3 q, where D is the volumetric mass density and N(q) s the Bose-Einstein distribution at temperature T .Fig. 10 shows that the vibrational frequencies used in our calculations yield for La 2 Zr 2 O 7 a specific heat in good agreement with experiments up to 1300 K. Above 1300 K, Specific heat of La2Zr2O7 and CsPbBr3.Dashed-blue lines are theoretical calculations at constant volume and using the bare phonon frequencies.Scatter points are measurements performed at constant pressure.For La2Zr2O7 these are taken from Vassen [107], Bolech [181], Zhang [114] and Sedmidubskỳ [182].For CsPbBr3 these are taken from Evarestov [183] and Elbaz [184].[185] for La2Zr2O7, and from Du et al. [186] for CsPbBr3.The experimental measurements on La2Zr2O7 in panel a) have been shifted rigidly of 10 cm −1 towards left at all temperatures, to ease the comparison with our calculations.The vertical lines in panel a) highlight that in La2Zr2O7 the experimental Raman peaks do not have an appreciable temperature-dependent shift, thus the aforementioned rigid shift is not compensating for the temperature-dependent renormalization of the vibrational frequencies not accounted for in our calculations (see text for details).Overall, the temperature-dependent broadening of the experimental Raman peaks is in remarkably good agreement with our calculations, thus the perturbative treatment of anharmonicity employed is accurate enough for our illustrative purposes.
experiments start deviating from theory, suggesting that at these very high temperatures a more refined treatment of anharmonicity is needed.For this reason, calculations in Sec.VII have been limited to the upper-bound temperature of 1300 K.For CsPbBr 3 , the theoretical and experimental specific heats agree reasonably well over the temperature range analyzed (which corresponds to the temperature interval in which the orthorhombic phase is stable).
To assess the degree of accuracy of our perturbative treatment of anharmonicity in La 2 Zr 2 O 7 and CsPbBr 3 , we use the theoretical anharmonic linewidths to predict the temperature dependence of the Raman spectrum of these complex crystals, and compare our predictions with experiments [185,186].Specifically, the Raman tensor of La 2 Zr 2 O 7 and CsPbBr 3 is computed following the approach detailed in Ref. [187].The Raman tensor and bare phonon frequencies at q = 0 are used to compute the modal Raman intensity I s for La 2 Zr 2 O 7 and CsPbBr 3 (we use the formula for the powder-averaged intensity discussed in Eq. ( 5) of Ref. [188], in such for-mula we consider laser frequencies equal to those used in the experiments [185,186] against which we compare our calculations).Finally, from the the Raman intensity I s the Raman spectrum is computed as [189]: 2 , where ω(0) s and Γ(0) s are the bare phonon frequencies and anharmonic linewidths at q=0, respectively, and the linewidth Γ ins =2 cm −1 accounts for the instrumental broadening (see e.g.Ref. [190]).Therefore, the Raman spectrum is simulated under the same approximations employed in our thermalconductivity calculations, i.e. it does not account for the temperature-dependent shift of frequencies, and considers anharmonicity up to the lowest third order.Fig. 11 shows that the simulated Raman spectrum agrees reasonably well with experiments, both for La 2 Zr 2 O 7 [185] and CsPbBr 3 [186] over the temperature range from 100 to 300 K. We note that in Fig. 11a) the experimental measurements in La 2 Zr 2 O 7 are shifted rigidly of 10 cm −1 at all temperatures, to ease the comparison with our calculations.The vertical lines in Fig. 11a) highlight that experiments in La 2 Zr 2 O 7 do not show an appreciable temperature-dependent shift of the Raman peaks in the temperature range analyzed, thus the shift we applied is not compensating for the temperature-dependent renormalization of the vibrational frequencies that is not accounted for in our calculations.Assessing the origin of such a systematic shift goes beyond the scope of the present study, since Fig. 10 shows that the bare phonon frequencies calculated for La 2 Zr 2 O 7 yield a specific heat in satisfying agreement with experiments and are thus accurate enough for our scope.Although Raman experiments in La 2 Zr 2 O 7 do not cover the entire temperature range over which the thermal conductivity is computed, we note that in the range from 100 to 300 K the temperature-dependent broadening of the experimental Raman peaks is in remarkably good agreement with our calculations, and at 300 K the coherences conductivity in La 2 Zr 2 O 7 is a non-negligible fraction of the total thermal conductivity.These considerations suggest that the perturbative description of anharmonicity employed is accurate enough to describe the main features of the temperature-thermal conductivity curve of La 2 Zr 2 O 7 .Finally, we note that often in the literature the inverse of the analysis discussed here is performed, i.e. the experimental values for the linewidths are extracted, via a fitting procedure, from the measured Raman spectrum and then compared with their theoretical counterparts.Such a procedure is error-prone in complex crystals such as La 2 Zr 2 O 7 and CsPbBr 3 .In fact, these complex crystals feature many Raman-active vibrational modes that are very close in frequency and have very different intensity, therefore in these cases very often the linewidths of individual vibrational modes cannot be resolved from the Raman measurements.For these reasons such a procedure has not been followed here, and instead the simulated Raman spectra has been compared directly with experiments.
In summary, for both La 2 Zr 2 O 7 and CsPbBr 3 , the bare phonon frequencies computed yield a specific heat in good agreement with that measured in experiments.Moreover, the anharmonic linewidths computed yield a temperature-dependent broadening of the simulated Raman peaks in good agreement with experiments over a temperature range in which the coherences conductivity is not negligible.This shows that the perturbative treatment for anharmonicity adopted is accurate enough for the illustrative scope of the present study.In this section we discuss the effects of the phase convention adopted for the Fourier transform on the thermal conductivity, comparing the thermal conductivity obtained using the velocity operator in the smooth phase convention (36) and that obtained using the velocity operator in the step-like phase convention (50).This com-parison is reported only for illustrative purposes.We remark that the velocity operator in the smooth phase convention (36) is the only one yielding a size-consistent conductivity -thus the convention to be used -and the results reported in Fig. 3 and Fig. 6 have been computed using the smooth phases.We show in Fig. 12 and Fig. 13 how the predictions of the conductivity of La 2 Zr 2 O 7 and CsPbBr 3 , respectively, change if the sizeinconsistent step-like phase convention is used in place of the correct (size-consistent) smooth phase convention.We note that in both La 2 Zr 2 O 7 and CsPbBr 3 the populations' terms of the conductivity κ P is not affected by the phase convention adopted, since the diagonal elements of the velocity operator do not depend on the phase convention (see Eq. ( 50)).The differences in the conductivities obtained in the two phase conventions are increasingly more relevant as temperature increases and arise exclusively from the coherences' conductivity κ C , with the step-like phase convention yielding a larger κ C .12. Effect of phase convention on the thermal conductivity of La2Zr2O7.Solid lines are the conductivities computed using the correct size-consistent smooth phase convention, dashed lines are the conductivities computed using the size-inconsistent step-like phase convention.The populations' conductivities are green (κP), and they are equal in the two phase conventions (see text).The coherences' conductivities (κC) are blue.The total conductivity (black) is given by the sum κTOT = κP + κC.Using the step-like phase convention yields a total conductivity larger than that obtained using the smooth phase convention.We recall that the populations and coherences conductivity tensors of La2Zr2O7 are proportional to the identity (see Sec VII A), thus populations and coherences conductivities are represented here as scalars.Effect of the phase convention on the thermal conductivity of CsPbBr3.Solid lines are the conductivities computed using the correct size-consistent smooth phase convention, dashed lines are the conductivities computed using the size-inconsistent step-like phase convention (also reported in Fig. 1 of Ref. [37]).The populations' conductivities are green (κ xx P ), and they are equal in the two phase conventions (see text).The coherences' conductivities (κ xx C ) are blue.The total conductivity (black) is given by the sum κ xx TOT = κ xx P + κ xx C .Using the step-like phase convention yields a total conductivity larger than that obtained using the smooth phase convention.
Appendix E: Ioffe-Regel limit in space and center of the non-sharp particle-wave crossover for phonons In this Appendix we report the details of the derivation of Eq. (53).We start by showing how to obtain the single-phonon contribution to the particle-like and wave-like (coherences) conductivities appearing in such an equation.The expression for the particle-like contribution K avg P (q) s is trivially obtained from the integrand of the populations (Peierls-Boltzmann) conductivity in the SMA (49) (we recall that the SMA is accurate for the ultralow-thermal-conductivity materials in focus here), K avg P (q) s =C(q) s α |v α (q) s,s | 2 3 1 Γ(q) s =C(q) s α v avg (q) s,s v avg (q) s,s Γ(q) s , where we have defined the spatially-averaged group velocity v avg (q) s,s = 1 3 3 α=1 |v α (q) s,s | 2 .The expression for the single-phonon contribution to the wave-like conductivity, K avg C (q) s , is obtained requiring that in the coherences' coupling between two phonons (q) s and (q) s (which is quantified by the normalized trace of the integrand of the coherence term in Eq. ( 47)), each phonon contributes to the coupling with a weight equal to the relative specific heat (e.g. for phonon (q) s the weight is C(q)s C(q)s+C(q) s , and correspondingly for phonon (q) s the weight is C(q) s C(q)s+C(q) s .In practice, denoting with k(q) αβ s,s the integrand of the coherences term in Eq. ( 47), and noting that for α=β this term is symmetric in the mode indexes, k(q) αα s,s = k(q) αα s ,s , we can rewrite the average trace of the coherences conductivity tensor as s,s =s B k(q) αα s,s d 3 q = 1 (2π) 3 s B s =s C(q) s C(q) s +C(q) s 1 3 α 2k(q) αα s,s d 3 q. (E2) The term in square brackets in Eq. (E2) has the following explicit form C(q) s C(q) s +C(q) s ω(q) s +ω(q) s 2 C(q) s ω(q) s + C(q) s ω(q) s × 1 3 α |v α (q) s,s | 2 1 2 Γ(q) s +Γ(q) s [ω(q) s −ω(q) s ] 2 + 1 4 [Γ(q) s +Γ(q) s ] 2 , (E3) and represents the contribution of the phonon (q) s to the wave-like (coherences) conductivity.

(E4)
We now proceed with performing an order-of-magnitude analysis of Eq. (E4).First, we note that the Lorentzian in Eq. (E4) is peaked around ω(q) s =ω(q) s , implying the following order-of-magnitude estimate for the term involving frequencies and specific heats: ω(q)s+ω(q) s C(q)s+C(q) s 1 2 C(q)s ω(q)s + C(q) s ω(q) s ∼1.Second, we rely on the considerations reported in Sec.VIII to approximate the ratio between the velocity-operator elements to one, α |v α (q) s,s | 2 α |v α (q)s,s| 2 ∼ 1.These considerations imply that the relation reported in Eq. (E4), with order-of-magnitude accuracy, can be simplified to: s =s 1 2 Γ(q) s +Γ(q) s [ω(q) s −ω(q) s ] 2 + 1 4 [Γ(q) s +Γ(q) s ] 2 .(E5) Now we distinguish what happens in Eq. (E5) in the two regimes of simple and complex crystals.In the simple-crystal regime, where the typical phonon interband spacings are much larger than the linewidths (∆ω avg Γ(q) s ∀q, s), Eq. (E5) assumes values much smaller than one on both sides and is trivially verified.In the complex-crystal regime, where the typical phonon interband spacings are comparable or smaller than the linewidths (∆ω avg < ∼ Γ(q) s ∀q, s), we need to perform two approximations to proceed with our order-of-magnitude analysis: (i) we consider the 3N at phonon bands as uniformly distributed over the frequency range (0, ω max ], (i.e. for a generic mode index s , we have ω(q) s n • ∆ω avg , where n ∈ [1, . . ., 3N at ] [191]); (ii) we approximate the average linewidth as Γ(q)s+Γ(q) s 2 =Γ(q) s .These approximations allow us to rewrite the sum in Eq. (E5) as an integral: [ω−ω(q)s] 2 +Γ 2 (q)s dω.Now we note that the integral of the Lorentzian distribution yields a number smaller than π and of the order of one.In fact, such an integral would yield π if the integration domain was running from −∞ to +∞, in the case considered here the linewidths are a fraction of the finite integration domain (see e.g.Fig. 4b) and thus a considerable portion of the Lorentzian is always integrated, yielding a value of the order of unity.These considerations imply that the two members in Eq. (E5) are approximately equal also in the case of a complex crystal, proving with order-of-magnitude accuracy Eq. ( 53).The crystal structure of La 2 Zr 2 O 7 is cubic with spacegroup Fd 3m [227] and has been taken from the materials project database [192] (id = mp-4974 [193]).The equilibrium crystal structure has been computed using density functional theory (DFT), specifically performing a "vc-relax" calculation with the Quantum espresso [194] software.The following functionals have been tested: LDA (with ultrasoft pseudopotentials from the psilibrary [195]), PBEsol (with GBRV [196] pseudopotentials), and PBE with Grimme D2 correction [197] (with pseudopotentials from the SSSP precision library [198,199]).A summary of the computational parameters and equilibrium lattice parameter is reported in Table II.The PBEsol functional has been selected on the basis of the best agreement between first-principles and experimental lattice parameters [200] (see Table II).Second-order force constants have been computed on a 4×4×4 mesh using  [200] 10.805 density-functional perturbation theory [2] and accounting for the non-analytic term correction due to the dielectric tensor and Born effective charges.Third-order force constants have been computed using the finite-difference method implemented in ShengBTE [4], on a 2×2×2 supercell, integrating the Brillouin zone with a 2×2×2 Monkhorst-Pack mesh, and considering interactions up to 7.8 Å.Then, harmonic force constants have been converted from Quantum espresso format to hdf5 format using phonopy [201], and third order force constants have been converted to hdf5 format using hiphive [5].
The LWTE conductivity formula (47) has been implemented in the phono3py software [6].Results reported in Fig. 3 have been computed using a mesh 19×19×19, and evaluating the scattering operator accounting for natural-abundance isotopic scattering [7,50] and thirdorder anharmonicity [3,6,10,11,44,45].The conductivity has been computed using the single-mode relaxation time approximation (SMA), since for La 2 Zr 2 O 7 the exact Peierls-Boltzmann thermal conductivity is known to be practically indistinguishable from the SMA value [42] (in agreement with expectations for materials with ultralow thermal conductivity [11,84]).The tetrahedron method has been used for the computation of the phonon linewidths.The data shown in Fig. 3 for κ P are compatible with the Peierls-Boltzmann conductivity of La 2 Zr 2 O 7 presented in Ref. [42].The data needed to reproduce the calculations reported here (DFT-relaxed primitive cell, second-order and third-order force constants) are available on the Open Science Platform Materials Cloud [202,203].

CsPbBr3
The crystal structure of orthorhombic CsPbBr 3 has been taken from Ref. [204] (Crystallographic Open Database [205] id 4510745) and converted to Pnmb spacegroup using VESTA [206].The equilibrium crystal structure has been obtained performing a "vc-relax" calculation with Quantum espresso [194], employing the PBEsol exchange-correlation functional and GBRV [196] pseudopotentials.The PBEsol functional is selected on the basis of its good agreement between first-principles and experimental lattice parameters [22,204] (see Table III) and its capability to accurately describe the vibrational properties of inorganic halide perovskites [22].Kinetic energy cutoffs of 50 and 400 Ry have been used for the wave functions and the charge density; the Brillouin zone has been integrated with a Monkhorst-Pack mesh of 5×5×4 points, with a (1,1,1) shift, and second-order force constants have been computed on a 4×4×4 mesh using density-functional perturbation theory [2] and accounting for the non-analytic term correction due to the dielectric tensor and Born effective charges.Third-order force constants have been computed using the finite-difference method implemented in ShengBTE [4], on a 2×2×2 supercell, integrating the Brillouin zone with a 2×2×2 Monkhorst-Pack mesh, and considering interactions up to 6 Å.Then, anharmonic force constants have been converted from ShengBTE format to mat3R format using the conversion software d3_import_shengbte.xprovided with the D3Q package [3,76], including interactions up to the third neighbouring cell for the re-centering of the third-order force constants (parameter NFAR=3).Thermal conductivity calculations have been performed with the software D3Q, using a mesh 8×7×5.The scattering operator has been computed using a Gaussian smearing of 2 cm −1 , considering third-order anharmonicity [3,11] and scattering due to isotopes at natural abundance [7,50].The conductivity is computed using the single-mode relaxation time approximation (SMA), since for this system the exact thermal conductivity is practically indistinguishable from the SMA value (in agreement with expectations for materials with ultralow thermal conductivity [11,84]).The data shown in Fig. 6 for κ αβ P are compatible with the theoretical results presented in Ref. [22] and related Supplementary Material.In particular, the direction xx in this work corresponds to the [010] direction in Ref. [22].Experimental data from Ref. [22] 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. [111].Therefore, we compared the bulk thermal conductivity predictions of equation (47) with the experimental data from Ref. [111].The results reported in Fig. 6 have been reproduced using the software phono3py [6] and the tetrahedron method (the softwares phonopy [201] and hiphive [5] have been used to convert force constants from Quantum espresso or ShengBTE format to the hdf5 format used by phono3py).An implementation of Eq. ( 47) will be made available as an update for the packages phono3py [6] and D3Q [3,11,76].The data needed to reproduce the calculations reported here (DFT-relaxed primitive cell, second-order and third-order force constants) are available on the Open Science Platform Materials Cloud [202,203].should not contain the term s = s .In the limit of a large number of bands, including in the sum also the term s=s does not alter the order-of-magnitude analysis performed here, and we will see soon that allows to approximate the sum with an integral.For these reasons the term s=s is included

4 √ G R b α ,Rbα and 4 √ G − 1
Rbα,R b α are the fourth roots of the mass-rescaled harmonic interatomic force-constant matrix (4) and its inverse; they satisfy the relations R ,b ,α 4 0b α is the mean square displacement due to the zero-point motion.The force constant matrix G Rbα,R b α goes to zero for |R+τ b −R −τ b |→∞ [2], and this implies that also 4 √ FIG. 1. â † (R) bα applied to the ground state |0 creates an atomic vibration centered on the atom at position R+τ b .The gray rectangles represent primitive cells of a lattice with basis, orange and gray circles symbolize atoms.The horizontal red lines show the mean square displacements of the atoms.The vertical lines have the same length of the horizontal lines below them, they show that in the excited state |ψ Rbα =â † (R) bα |0 (where |0 is the ground state), atoms at position R +τ b have a mean square displacement (MSD) that becomes smaller as the distance |R+τ b −R −τ b | increases.The envelope of the atomic MSD (vertical red lines) is decomposed into pink and blue, to emphasize that the total MSD in the excited state is the result of the sum of two terms: a space-dependent MSD due to the excitation created by â † (R) bα (pink), and a spaceindependent MSD due to the zero-point motion (blue rectangle).
d p n V 0 5 j p O 1 a i c h d o D K + F P 4 G D 6 D T + E N O 0 m 3 7 b a o n S j S a C 5 n Z s 6 M g 4 x R I f v 9 P z O z j b n m t e v z N 1 o 3 b 9 2 + c 3 d h 8 d 6 e S I s c k 1 2 c s g 6 u 2 0 x M l Y m p G t o y o k 7 z 3 W r u q i j 7 B 9 S F s l f O r p 7 q i 4 E w z W r l r M K P L w 3 X b s C g 4 Z J R T K U Z K i S 5 s d 6 s 5 A K S J n H g g R 3 J s o A 7 V V j 2 m N Y h I w Q + c x K g m g d o x q t a r p i q 4 c o Q J 9 u W D l 4 H H w 9 d 5 k / x L i Z j U PF 3 R B + F w v X S Z G s f s r G j X m 9 7 U c u m r W T u 3 N U v o t N H 2 t G Z u 5 c q 7 7 O r R Q q e / 2 i 8 F n F e 8 W u k 4 t e y M F m d + w z D F B S e J x A w J M f D 6 m f S V v U z M i G 7 B Q p A M 4 U M U k 4 F R E 8 S J 8 F X 5 k D R 4 b C w h i N L c / I k E p f V 0 h k J c i C M e m M h y q 9 M + a 7 z Q F / C p y j J 6 4 S u a Z I U h A 1 e F o 4 I B m Q L 7 S k F I c 4 I l O z I K w j k 1 v Q M 8 R o Z q a d 5 y C 7 4 h Z r a c v D N 1 X r F s j A I i V X 1 n V G o F D W Y a 0 i Q O S Y Q K Z i 0 i m u h c K 8 F M O 5 + J / B 9 A k L L w M p j g R 4 2 T k O 9 m m xw l Y X V l e u D 5 5 V I n g K r j a X 1 B p D 2 C K v p M p F m 6 N 7 3 i 8 8 r e 2 q q 3 s f r s 4 0 Z n c 6 t e / 7 z z w H n k u I 7 n P H c 2 n W 1 n x 9 l 1 c G O u 8 a S x 3 t i Y + 9 t 8 2 F x q u l X o 7 E y d c 9 8 5 I 0 3 v H 5 v K 3 o k = < / l a t e x i t > W 8 I / o r P 4 D v 4 A e w k 3 V 5 R c W T p a C 5 n 5 n g m Q S a 4 0 r 3 e 7 4 n G h e b F S 5 c n r 7 S u X r t + 4 + b U 9 K 1 N l R Y 5 Z R s 0 F W m + H R D F B E / Y h u Z a s O 0 s Z 0 Q G g m 0 F u 8 + d f + s T y x V P k / d 6 L 2 M D S e K E R 5 w S b U 3 D 6 Y k / y J 7 7 O M o J B c A Z y T U n A h m s W a J o K j N I T B e 3 A Q c S 3 p o 5 3 H b g o w P a z A 5 B G Y T b 5 m g i Q t p e h A M e x / t D w CO i Y d X s A M 4 l y l h u j K W w H z w 1 8 N D g l 0 R K 0 q 0 4 S z q 8 w u P + f x d 3 F K m V 5 9 T D Y d Z r Y 7 p Q x R 6 y D j B u l V r h g a n k + g a 6 C z j j s z u L r i t V S C y 4 5 F o N A d Q M b s 9 U z S L M E z 3 2 Y E n 0 y K b v w k q t xR l U B H h N s p j U S r l r v G q 2 a q S i c 8 3 M j b n P l 1 o G H s g 9 y D t X 9 L j W 0 U o D F O 4 s l i 7 c H k 5 1 e v O 9 8 q D T w K 9 B x 6 v P u l 2 T X z h M a S F Z o q k g S v X 9 X q Y H 4I Z O B T M t X C i W E b p L Y t a 3 M C G S q Q G U 6 2 n Q P W s J U Z T m 9 i Y a l d a j G U C k U n s y s J H l W 5 7 0 O e O Z v k C e q K y j J w P g S V b Y V 6 F V 4 a g Q S K f I 7 T 4 K e c 6 o F n s W E J p z 2 z u i I 2 J X Q d s / p I V f M K s t Z 6 9 s n W c i G 5 G Aa a i n y 7 U B b D n T k C d x y C J S C G d R 0 R h L A 0 r Y d t 4 x / S + C I B X h e T T B l 5 o n Y Z / t W C V J w m r G p u 8 P y v G O C a H j G 3 N G p N u G K v p Y p B 2 6 f 3 L E p 8 H m w r y / N P / o z V J n e a U e / 6 R 3 x 7 v r d T 3 f e + w t e 6 v e u r f h 0 c a H x r f G j 8 b P 5 l q z a O 4 3 v 1 a h j Y k 6 5 7 Z 3 7 D S / / w U y H K T + < / l a t e x i t > FIG. 2.Size consistency of the LWTE conductivity in the smooth phase convention.Size consistency is tested computing the conductivity (47) of a silicon crystal at T =300 K and using a fictitious linewidth of Γ(q)s=30 cm −1 ∀ q, s (see text).We sketch in panels a,b) two equivalent crystal representations: a) repeating a 2-atom primitive cell the same number of times along each Cartesian direction; b) repeating a 16-atom unit cell that is larger than the primitive.Panel c) shows that the smooth phase convention (used in Eq. (7) yields size-consistent numerical results, since the conductivity is unchanged when the silicon crystal with fictitious linewidths is modeled as: (i) a 16×16×16 repetition of the 2-atom primitive cell; (ii) a (physically equivalent) 2×2×16 repetition of a unit cell that is a supercell 8×8×1 of the primitive.Panel d) shows instead that the steplike phase convention (used in Eq. (9)) is not size-consistent, since it yields different conductivities for the two physically equivalent representations of the crystal.The black and red lines serve as a guide to the eye, to show that the Peierls-Boltzmann (populations) conductivity does not depend on the phase convention adopted (black), while the coherences conductivity depends on the phase convention adopted (red).

FIG. 4 .
FIG. 4.Vibrational properties and heat conduction mechanisms in La2Zr2O7.Panels a(i), b(i) show the phonon spectrum of La2Zr2O7 (gray lines) with the phonon linewidths (the full amplitude of the shaded gray areas correspond to the full phonon linewidth Γ(q)s) at 200 K and 1300 K, respectively.Colored circles represent the phonon eigenstates (q)s sampled in the calculation; the area of each circle is related to its contribution to κ (area ∝ √ contribution to κ[115]) and the color identifies the origin of the contribution: green for populations' propagation and blue for coherences' anharmonic couplings (red corresponds to 50% of each, see Eq. (54) and related discussion for details).In the coherences' couplings between two modes (q)s and (q) s the contribution of the single mode s is determined as C(q)s/[C(q)s+C(q) s ], where C(q)s is the specific heat of a given phonon (see Appendix E for the details).Panels a(ii), b(ii) show the thermal conductivity density of states κω of populations ("P", green) and coherences ("C", blue) at 200 K and 1300 K, respectively.The gray line shows the vibrational density of states (VDOS).Panels a(iii), b(iii) show the cumulative total thermal conductivity (in black) as a sum of the populations' contribution ("P") in green and coherences' one ("C") in blue, at 200 K and 1300 K, respectively.At 200 K (panel a) La2Zr2O7 is a simple crystal with interband spacings larger than the linewidths and κP kC , while at 1300 K it is a complex crystal with interband spacings smaller than the linewidths and κP < ∼ kC .
FIG. 7.Regime diagram for thermal transport.Three regimes can be distinguished in a diagram reporting the phonon lifetime (τ ) as a function of the phonon frequency (ω).The horizontal black like is the Wigner limit in time, τ = 1/∆ωavg (i.e. the inverse average interband spacing, see Eq. (52), and denotes the center of the non-sharp particle-wave crossover for phonons.The green region well above the line represents phonons that propagate particlelike and mainly contribute to the populations conductivity; the blue region below the line represents phonons that tunnel wave-like and mainly contribute to the coherences conductivity; the red points around the line are phonons that contribute simultaneously and comparably to the particle-like and wave-like conductivities.The purple line τ =1/ω is the center of the non-sharp Ioffe-Regel limit in time[41,121].The region above it represents well defined phonons; the Wigner transport equation (37) describes all these phonons (blue, red, green), while the semiclassical Peierls-Boltzmann equation accounts only for phonons that propagate particlelike (green).The yellow region below the Ioffe-Regel limit represents overdamped phonons, which require full spectralfunction approaches[78,79] to be described correctly.

La 2 FIG. 11 .
FIG. 11.Raman spectra of La2Zr2O7 and CsPbBr3.Solid lines are temperature-dependent simulations performed relying on the same perturbative treatment of anharmonicity that we employed in out thermal-conductivity calculations (see text for details).Dashed-black lines are experiments, taken from Paul et al.[185] for La2Zr2O7, and from Du et al.[186] for CsPbBr3.The experimental measurements on La2Zr2O7 in panel a) have been shifted rigidly of 10 cm −1 towards left at all temperatures, to ease the comparison with our calculations.The vertical lines in panel a) highlight that in La2Zr2O7 the experimental Raman peaks do not have an appreciable temperature-dependent shift, thus the aforementioned rigid shift is not compensating for the temperature-dependent renormalization of the vibrational frequencies not accounted for in our calculations (see text for details).Overall, the temperature-dependent broadening of the experimental Raman peaks is in remarkably good agreement with our calculations, thus the perturbative treatment of anharmonicity employed is accurate enough for our illustrative purposes.
Appendix D: Effects of the phase convention on the thermal conductivity FIG.12.Effect of phase convention on the thermal conductivity of La2Zr2O7.Solid lines are the conductivities computed using the correct size-consistent smooth phase convention, dashed lines are the conductivities computed using the size-inconsistent step-like phase convention.The populations' conductivities are green (κP), and they are equal in the two phase conventions (see text).The coherences' conductivities (κC) are blue.The total conductivity (black) is given by the sum κTOT = κP + κC.Using the step-like phase convention yields a total conductivity larger than that obtained using the smooth phase convention.We recall that the populations and coherences conductivity tensors of La2Zr2O7 are proportional to the identity (see Sec VII A), thus populations and coherences conductivities are represented here as scalars.
FIG.13.Effect of the phase convention on the thermal conductivity of CsPbBr3.Solid lines are the conductivities computed using the correct size-consistent smooth phase convention, dashed lines are the conductivities computed using the size-inconsistent step-like phase convention (also reported in Fig.1of Ref.[37]).The populations' conductivities are green (κ xx P ), and they are equal in the two phase conventions (see text).The coherences' conductivities (κ xx C ) are blue.The total conductivity (black) is given by the sum κ xx TOT = κ xx P + κ xx C .Using the step-like phase convention yields a total conductivity larger than that obtained using the smooth phase convention.

TABLE I .
Comparison between the Peierls-Boltzmann and Wigner frameworks.