Thermal transport in crystals as a kinetic theory of relaxons

Thermal conductivity in dielectric crystals is the result of the relaxation of lattice vibrations described by the phonon Boltzmann transport equation. Remarkably, an exact microscopic definition of the heat carriers and their relaxation times is still missing: phonons, typically regarded as the relevant excitations for thermal transport, cannot be identified as the heat carriers when most scattering events conserve momentum and do not dissipate heat flux. This is the case for two-dimensional or layered materials at room temperature, or three-dimensional crystals at cryogenic temperatures. In this work we show that the eigenvectors of the scattering matrix in the Boltzmann equation define collective phonon excitations, termed here relaxons. These excitations have well defined relaxation times, directly related to heat flux dissipation, and provide an exact description of thermal transport as a kinetic theory of the relaxon gas. We show why Matthiessen's rule is violated, and construct a procedure for obtaining the mean free paths and relaxation times of the relaxons. These considerations are general, and would apply also to other semiclassical transport models, such as the electronic Boltzmann equation. For heat transport, they remain relevant even in conventional crystals like silicon, but are of the utmost importance in the case of two-dimensional materials, where they can revise by several orders of magnitude the relevant time- and length-scales for thermal transport in the hydrodynamic regime.

In particular, the linear BTE can nowadays be solved exactly, using empirical or first-principles interactions, with iterative [4,16,17], variational [7,[18][19][20] or direct diagonalisation algorithms [21][22][23], all of which do not need to simplify the scattering operator with the often-adopted single-mode relaxation time approximation (SMA). In the SMA, each phonon mode relaxes independently to equilibrium, and it has long been known that this is an incorrect assumption for solids at low temperatures [18,24]. Importantly, this approximation fails dramatically in lower dimensions, as first found in graphene [25], boron nitride and other two-dimensional (2D) materials [26][27][28], as well as in layered crystals [29]. The origin of such failure of the SMA has been up to now a matter of debate, with an emerging picture of collective phonon excitations being responsible for heat transfer [26,[29][30][31][32][33], while nevertheless lacking a definition of such excitations.
The microscopic interpretation of thermal transport in the BTE is based on the kinetic theory of gases, used in various contexts since its development in the 19 th cen- * nicola.marzari@epfl.ch tury, which relates the thermal conductivity to the velocities and relaxation times of the carriers, with phonons being usually identified as the relevant gas of excitations. However, as argued below, this identification is incorrect, since only the adoption of the SMA allows for the definition of a time interval (i.e. a lifetime) between heat flux dissipation events taken as phonon scatterings. Going beyond the SMA, the full, exact solution of the BTE provides the correct thermal conductivity (dramatically improved in 2D materials or at low temperatures), but adds complexity in its interpretation. In fact, solving exactly the BTE implies abandoning the concept of phonon relaxation time and the description of heat being carried by a gas of phonons. In other words, phonon lifetimes or phonon mean free paths are no longer relevant quantities to describe thermal transport, since phonons are not the heat carriers anymore. Yet, the beauty of the SMA description lies in the simplicity of its description of transport. A natural question then arises: can one define heat carriers, relaxation times, or mean free paths within the exact treatment of the BTE?
In this work, we provide an answer to these questions. First, we recall why phonon lifetimes are unrelated to heat flux dissipation. Then, we define a set of collective excitations, termed relaxons, that diagonalize the scattering matrix. The BTE is rewritten in the basis of these relaxons; in this representation, each eigenvector represents a collective excitation consisting of a linear combination of out-of-equilibrium phonon populations, and it describes the thermal relaxation of a collective excitation of out-of-equilibrium lattice vibrations. We show that each relaxon is characterised by a well defined relaxation time; in the case of a homogeneous system at the steady-state each relaxon has also a well defined velocity and mean free path, and the thermal conductivity can be interpreted exactly as a kinetic theory of the relaxon gas. As a practical example, we compare thermal conductivities in graphene and silicon contrasting the relaxon and the phonon representations, and highlight the profoundly different pictures that emerge.

II. APPROXIMATED RELAXATION TIMES
We start our derivation by recalling the microscopic description of heat transport given by the linearised phonon BTE [18]: (1) This equation describes the out-of-equilibrium dynamics of the phonon excitation number n µ at position x and time t, for all possible phonon states µ (in shorthand notation µ ≡ (q, s), where q varies over the Brillouinzone and s over the phonon branches). Furthermore, v µ is the phonon group velocity, V is a normalization volume, Ω µµ is the linear phonon scattering operator and ∆n µ = n µ −n µ is the deviation of the phonon distribution from thermal equilibrium, i.e. the Bose-Einstein distributionn µ (x, t) = (e ωµ/k B T (x,t) − 1) −1 , with ω µ being the phonon frequency and T (x, t) the local temperature. This linear approximation, commonly used in most studies of transport, allows to describe scattering as a linear operator represented by the action of the matrix Ω µµ on ∆n µ : this assumption holds for small deviations from thermal equilibrium and will always be used in the rest of the manuscript. The scattering matrix appearing in Eq. (1) is in its most general form and describes all possible mechanisms by which a phonon excitation can be transferred from a state µ to a state µ . For the sake of simplicity, we will limit this manuscript to the inclusion of three-phonon processes and isotopic scattering events [7], whose expressions are reported for completeness in Appendix A.
For later convenience, it is useful to write the left-hand side of Eq. (1) in terms of the unknown ∆n µ : where T is the reference temperature at which the BTE has been linearized. To obtain Eq. (2), we substituted n µ =n µ + ∆n µ in Eq. (1) and used the fact that the Bose-Einstein distribution depends on space and time only through the temperature T (x, t).
A closed-form solution of the above equation can be obtained in the SMA, which replaces the scattering operator with its diagonal terms To show that in this simplified diagonal form τ SMA µ represents indeed a relaxation time, let's consider a system at thermal equilibrium (T (x, t) = T ), so that the phonon distribution isn µ everywhere and thus in Eq.
(2) ∇(∆n µ ) = 0, ∂T ∂t = 0 and ∇T = 0. If we excite a single phonon at time t 0 , its population relaxes back to equilibrium as ∆n µ (t) = (n µ (t 0 ) −n µ )e −t/τ SMA µ , i.e. with a characteristic time τ SMA µ . The thermal conductivity tensor k ij (i and j are cartesian indices) is defined as the ratio between a heat flux Q i and a static gradient of temperature (∇T ) j . Two simplifications apply in this case. First, a steady-state condition allows to simplify the BTE by setting time derivatives to zero. Second, the spatial gradient can be simplified taking ∇(∆n µ ) = 0. This assumption, frequently adopted in literature, holds for a homogeneous perturbation of a bulk crystal (as in our case): if we apply a thermal gradient to a crystal at temperature T , the response ∆n µ should not depend on the particular position x inside the sample. Although we will not consider it further here, we note that this assumption cannot be applied when studying systems that break translational invariance, involving e.g. surfaces or pointlike heat sources. Under these conditions and the SMA, the resulting BTE can be solved analytically and, using the harmonic approximation for the heat flux Q = 1 V µ ω µ v µ ∆n µ [34] and the definition Q i = − j k ij (∇T ) j , the thermal conductivity is given by where (Λ j µ ) SMA is the component of the phonon mean free path in direction j. This expression can be interpreted as the thermal conductivity of a gas of phonons, each carrying a specific heat C µ = 1 k B T 2nµ (n µ + 1)( ω µ ) 2 = ∂nµ ∂T ω µ , traveling at velocity v i µ and with a mean free path (Λ j µ ) SMA = v j µ τ SMA µ before being thermalised by scattering. Crucially, the definition of phonon lifetime or mean free path cannot be extended beyond the SMA, since the off-diagonal terms of the scattering operator introduce couplings between phonons, and phonon thermalisation stops being governed by an exponential relaxation.

III. RELAXONS
An exact definition of relaxation times has been formally derived by Hardy [22], as an auxiliary result in his study of second sound. To recall it, let's first note that the left side of the BTE in Eq. (2) has a drifting operator diagonal in µ, whereas the right side has a scattering operator (determining scattering time scales) that is non-diagonal. To identify meaningful scattering times, we proceed with a change of basis that diagonalises the scattering operator while allowing the drifting term to become non-diagonal. To make more apparent the symmetries within the BTE, we perform the transformations [22,23,35,36]: These transformations are introduced to scale quantities appearing in the BTE in such a way thatΩ µµ =Ω µ µ (the matrix Ω does not obey this symmetry; see Appendix A for a detailed explanation). We note that sometimes these transformations appear in literature in the form of hyperbolic sines, by means of the identity sinh( . SinceΩ is a real symmetric matrix, it can be diagonalized, giving eigenvectors θ α µ and real eigenvalues 1 τα such that where α is the eigenvalue index. In passing, we define the scalar product α α ≡ 1 V µ θ α µ θ α µ that allows to define the orthonormalization condition for the eigenvectors ( α α = δ αα ) which will be helpful in the next algebraic operations. It can be shown [22,35] thatΩ is positive-semidefinite, i.e. 1 τα ≥ 0 ∀α, and that its eigenvectors are either even or odd, i.e. θ α µ = ±θ α −µ , where −µ = (−q, s) [22]. Little else is known on the eigenvalue spectrum of Eq. (7), which therefore has to be characterized numerically. In contrast with Refs. [21,22], we remark that the Bose-Einstein distribution is not an eigenvector with zero eigenvalue: the scattering operator acts only on the deviation from equilibrium ∆n µ , therefore thermal equilibrium (∆n µ = 0) is a stationary solution ( µ Ω µµ ∆n µ = 0) because it's an algebraically trivial solution. However, the Bose-Einstein distribution allows the introduction of a vector of unitary length where C = 1 V µ C µ , describing the increase of temperature. This vector is constructed as the linear deviation from equilibrium ofn µ (T + δT ), transformed using Eq. (6) and normalized to one. Note that θ 0 µ is not an eigenvector and doesn't have to be orthogonal to other eigenstates α.
Any response ∆ñ µ can be written as a linear combination of the θ α µ eigenvectors [22] ∆ñ µ ( and the BTE can be written in this θ α basis (to this aim, substitute Eq. (9) in (2) and take the scalar product of the equation with a generic eigenvector α ), becoming: where derives from the action of the diffusion operator on the deviation from equilibrium, while V α derives from the action of the diffusion operator on the equilibrium distribution.
The physical picture encoded in Eq. (10) underlines one of the key statements of this work: by diagonalising the scattering operator, the information about the characteristic relaxation time of the thermal excitations is now given by the eigenvalues 1 τα . The eigenvectors θ α µ for which this picture emerges represent collective excitations which we call here relaxons. Each relaxon represents a distribution of phonon excitation numbers (a wave packet), describing how the phonon distribution is relaxing to equilibrium. The coefficients f α are the relaxon occupation numbers, which are determined by the BTE in the out-of-equilibrium state, and that at equilibrium are all 0, so that the deviation from equilibrium ∆n µ vanishes.
The name relaxon is easily justified by considering a system at thermal equilibrium, so that ∂T ∂t = 0, ∇T = 0 and, in terms of relaxons, all states are empty everywhere, so that ∇f α = 0 ∀α in Eq. (10). If we excite a single relaxon α at time t 0 , its occupation will relax back to equilibrium as f α (t) = f α (t 0 )e −t/τα , therefore endowing τ α with the meaning of a relaxation time. Although the theory allows for zero eigenvalues, we will find in our examples only strictly positive relaxation times, so that all relaxons decay to zero for t → ∞. Using Eq. (9), one can show instead that phonon populations do not have well-defined relaxation times: since each phonon population decays as a linear combination of relaxon processes ∆n µ (t) = α f α (t 0 )θ α µ e −t/τα , the characteristic time for the decay depends on the initial conditions of the thermal excitation and can even display damped oscillations. Conversely, let's suppose to excite only one phonon mode µ at time t 0 . This initial state can be decomposed as the sum of different relaxons, each evolving with a different relaxation time. Therefore, at a subsequent time t one will also observe that new phonon modes µ = µ have been excited out of thermal equilibrium. In Fig. 1 we graphically illustrate our interpretation: each relaxon is a collective excitation of phonons, which interact through scattering events among themselves, but are decoupled from phonons belonging to different relaxons; owing to their positive relaxation time, relaxons disappear at long times, allowing the system to reestablish equilibrium.
Velocities appear in Eq. (10) with a matrix V αα coupling different relaxons, since it is the phonon basis that diagonalises the drifting operator; therefore, one cannot always identify a relaxon velocity. However, if we suppose to work in an infinite crystal at temperature T and apply a temperature perturbation homogeneously to the entire crystal, the response of the system is constant throughout the space and we can set ∇f α = 0, ∀α. Therefore, the BTE simplifies to (11) Both the drifting and the collision operator are now diagonal in α, and V α identifies a well-defined relaxon velocity.
Let's simplify the problem further and consider steady state. In this case, time derivatives are set to zero in Eq. (11) and one can, for small deviations from equilibrium, search for linear solutions of the form where i is a cartesian direction. The BTE reduces to whose solution for f i α is trivial. Using the relation between phonons and relaxon occupation numbers ∆n µ = n µ (n µ + 1)( iα ∇ i T f i α θ α µ ), we obtain the thermal conductivity where we introduced the relaxon mean free path Λ α (Λ j α is the component of Λ α in direction j). Therefore, the exact thermal conductivity in Eq. (13) is expressed in the framework of the kinetic theory of gases, and thermal transport can be thought as a flux of relaxons, each carrying a specific heat C, traveling at velocity V α for an average distance of Λ α before thermalisation occurs. At variance with the phonon picture, where each phonon participates to thermal conductivity with a mode specific heat C µ , all relaxons contribute with the same specific heat of the crystal C. Mathematically, the phonon-mode specific heat is moved in the vector θ 0 µ (note that (θ 0 µ ) 2 = Cµ C ) and thus is included into the relaxon velocity V α = 0 v α . To physically interpret this difference we recall that, from a thermodynamic point of view, the quantity CδT is the energy needed by the system to change temperature by δT . To observe such temperature change, all phonon modes must simultaneously change their occupation number according to the collective excitation θ 0 µ . The quantity C µ δT is the decomposition of such energy change in terms of each phonon mode. However, as explained before, one cannot suppose to excite a single phonon mode and bring it to a higher temperature without affecting the rest of the phonon ensemble: phonon scattering would redistribute the energy excess of such mode to the rest of the system. Only a collective excitation of phonons (θ 0 µ ) leads to a temperature change and the total energy cost for increasing temperature is necessarily associated with C; thus, from a thermodynamic point of view, one could state that the mode specific heat C µ does not have a well-defined meaning.
It's worth to point out the role played by the parity of relaxons. The quantity V α = 0 v α involves the odd function v µ (−v µ = v −µ ) and the even function θ 0 µ (owing to ω µ = ω −µ ). Therefore, relaxon velocities V α are different from zero only for odd relaxons α. Consequently, Eq. (12) implies also that only odd relaxons are excited in the steady state condition, and thus contribute to heat flux, while even relaxons have zero occupation number. The role of parity is reversed for determining the energy of the system, since the change from equilib-rium energy ∆E is In this case, even relaxons have a non-zero coefficient 0 α and contribute to an energy change, but for odd relaxons 0 α = 0 and thus do not change energy. We can thus deduce that at the steady state defined by Eq. (11), where only odd relaxons are excited, the energy of the system is conserved.

IV. GRAPHENE
As a first numerical example supporting these conclusions, we study relaxons in graphene, the material with the highest known thermal conductivity [37], and contrast the phonon and the relaxon pictures at 300K. Due to its symmetry, graphene's k ij tensor is diagonal and, since k xx = k yy , it has only one independent component (verified also numerically); therefore in the following we will drop cartesian indices and compute quantities numerically along the zig-zag direction. To proceed, we calculate harmonic and anharmonic force constants using density-functional perturbation theory [38][39][40][41][42][43][44] as implemented in the Quantum-ESPRESSO distribution [45] and construct the scattering matrix using 3-phonon and isotope-phonon interactions. The diagonalisation of Eq. (7) provides all the relaxon eigenvectors θ α µ ; each of them represents, at fixed relaxon index α, a difference in phonon populations with respect to thermal equilibrium (provided that is back-transformed using Eq. (6)). Notably, only few θ α µ have large relaxation times; for example, the longest-lived relaxon (α=1) is plotted in Fig.  2 as a function of the phonon index µ = (q, s). The first three (s = 1, 2, 3) branches are shown, corresponding to the out-of-plane, transverse or longitudinal acoustic phonons (ZA, TA and LA respectively). This particular relaxon induces a population difference for the ZA branch mainly located close to the Brillouin zone center, whereas TA and LA modes are altered throughout the Brillouin zone. The variations of optical modes (s = 4, 5, 6 not shown) are an order of magnitude smaller. The complex landscape drawn by these phonon distributions reflects the fact that out-of-equilibrium lattice properties cannot be described in terms of single phonon properties, as the action of scattering tightly couples phonons of any wavevector and branch.
We analyse the entire phonon and relaxon spectrum in Fig. 3, where the contributions to the SMA or the exact thermal conductivities are plotted as a function of the relaxon or phonon relaxation times. The thermal conductivities computed in the two pictures differ significantly in FIG. 2. Representation of the relaxon θ α µ with the longest relaxation time (α=1) in graphene at room temperature as a function of the phonon index µ = (q, s), where we choose s to be the out-of-plane/transverse/longitudinal acoustic mode (ZA, TA and LA respectively). We recall that the relaxon is a difference in phonon populations with respect to thermal equilibrium: overpopulated modes are colored in red, depopulated ones are in blue. The fine structure of the ridges is a numerical artifact due to discrete Brillouin zone sampling.
graphene (in this work we compute 3894 W/mK with the exact BTE against 495 W/mK with the SMA), hence, to have more comparable quantities, we plot the percentage contribution to thermal conductivity. We first note that the spectrum of phonon lifetimes (and phonon velocities and mean free paths) is continuous, with a divergence τ µ → ∞ for acoustic ZA phonons at the Γ point [46]. This divergence cannot be accurately described with a finite mesh of points sampling the Brillouin-zone (in our case, a full mesh of 128×128 points), resulting in a sparse tail of long-lived phonons on the right side of Fig. 3 3. Relaxation times and their contribution to the thermal conductivity of graphene at room temperature, considering relaxons or phonons as heat carriers. Relaxons tend to be longer lived than single phonon excitations, with large contributions to thermal conductivity coming from excitations with relaxation time larger than 10 3 ps, whereas phonons have lifetimes mainly in the range 10-100 ps. The shaded area is a guide to the eyes to stress that phonos form a continuous spectrum, while relaxons are discretized: thermal conductivity can be accurately described using a small number of relaxons. whose contribution to k SMA is negligible [46]. Instead, relaxation times for relaxons are discrete and sparse, in particular in the region of large values, so that only a small number of relaxons is sufficient to describe thermal transport with high accuracy. This observation is robust with respect to Brillouin Zone sampling: an improvement in the integration mesh has new phonon modes appear in the long-lifetime region; instead the longest relaxon relaxation times converge -from above -to the discretized values shown in figure. On average, relaxon relaxation times are larger by at least two orders of magnitude with respect to phonon lifetimes. The large difference between the time scales of phonons and relaxons appears because a single phonon scattering cannot thermalize the system [18], as instead implied by the SMA. Therefore, while phonons scatter at timescales of about 10-100 ps, heat flux is dissipated by relaxons within nano-and microsecond timescales.
Before analysing velocities, we note that the sign of V α is arbitrary, since both θ α µ and −θ α µ are relaxon eigenvectors. As a convention, we select the sign of odd eigenvectors such that V α is non-negative (and so also Λ α ), noting that in any case the contribution to k would be positive (as V 2 α ). Phonon velocities can also assume both signs: in figure we plot their absolute values. Fig. 4 reveals that the velocities of relaxons are much smaller than those of phonons: while the scale of phonon velocities is set by the speed of sound (the group velocity of the longitudinal acoustic phonon is about 20 km/s), relaxons are slower by two orders of magnitude, indicating that heat is transferred through the material at 0.1-1 km/s. Finally, we show the relaxon mean free paths in Fig. 5 (projected along the transport direction). As other firstprinciples studies reported [8], phonon mean free paths for graphene are distributed in the 0.1-1µm region [29]; this is confirmed here. For relaxons, most contributions to thermal conductivity come with mean free paths above 0.1 µm, the longest and most important contributions having mean free paths up to tens of µm. The contribution to k is roughly monotonic with the mean free path, and the large increase in τ α is partly compensated by the decreased V α . The saturation of relaxons' mean free paths at tens of µm appears to be in contrast with recent estimates for saturation lengths of 100 µm [29] or longer [47]; we will comment this discrepancy after discussing the next example. Analysis of phonon and relaxon mean free paths and their contributions to the thermal conductivity of graphene at room temperature. The spectrum of phonon mean free paths is peaked in the submicron region, whereas relaxon mean free paths are skewed to values larger than 1µm. The two largest relaxon mean free paths, illustrating the maximum distance that heat flux can propagate inside the material before decaying, lie between 10 and 100µm.

V. SILICON
As a further test for relaxons, let us now turn our attention to silicon and examine its thermal transport properties at room temperature. The thermal conductivity tensor in silicon is diagonal and the three cartesian directions are equivalent; we therefore only consider transport properties along the (100) direction. At variance with graphene, here the SMA introduces an appealingly small error: in our calculations we find 138 W/mK instead of 141 W/mK for the exact solution; these estimates are in line with previous first-principles studies [5,6,8]. The small difference between the two pictures is somehow replicated in their relaxation times, reported in Fig. 6. The time-scales covered by phonons and relaxons covers approximately the same range of values, except for one relaxon, not shown in the graph, that has relaxation time of 2 · 10 5 ps but negligible contribution to thermal conductivity (10 −9 %). However, one can note that the two distributions of values do not perfectly overlap: even if the scattering matrix is diagonally dominant, there are anyway small non-zero out-of-diagonal matrix elements that introduce deviations from the SMA.
It is enlightening to analyze the different velocity scales set by the two pictures, depicted in Fig. 7. Once again, most of the contributions to SMA thermal transport come from phonons with velocities close to the speed of sound, which in silicon is approximately 8 km/s. However, relaxon velocities are two orders of magnitude smaller than this limiting value, reaching merely 60 m/s. Despite the fact that the instantaneous velocity of a lattice vibration is determined by the phonon dispersion, the velocity at which the heat flux propagates can be much different: the scattering between phonons slows it down. The mean free paths for relaxons and phonons in silicon are compared in Fig. 8. The vast difference originating from the velocities is carried over, so that while mean free paths of phonons extend up to 100µm, in agreement with other first-principles studies [8,48,49], relaxons travel for a distance two orders of magnitude smaller than phonons. It seems therefore puzzling that the two pictures give such large differences in the estimate of velocities and mean free paths, despite the fact that thermal conductivities are essentially identical. To explain this discrepancy, let's compare Eqs. (4) and (13) for thermal conductivity and recall that specific heat is constant for relaxons and mode-specific for phonons. Also the SMA conductivity can be written in a form with constant spe- cific heat for each phonon, provided that we rescale velocities as v µ → Cµ VC v µ (consequently, also the mean free path is scaled Λ SMA . After this transformation, velocities and mean free paths of phonons and relaxons in silicon are again within the same order of magnitude, although residual discrepancies still persist (see Appendix B for the spectra after rescaling). Therefore, the differences observed in silicon between the two pictures arise mainly from the different interpretation of specific heat.
We note that other experimental and theoretical efforts have estimated heat mean free paths in silicon [8,[48][49][50][51][52], obtaining values that are comparable to those of the phonon mean free paths, and different from the relaxon mean free paths presented here. It is important to stress, though, that at least one of these two assumptions were used: first, that results can be interpreted considering phonons as the heat carriers; and second, that surface or grain boundary scattering can be exploited as a tool to estimate heat mean free paths. In the present work, we discussed already at length the limitations of the former assumption; as for the latter one, we cannot compare our results with studies that rely on surface scattering, since the data presented here pertain to a homogeneous bulk crystal. It is nevertheless possible to drop the bulk condition and solve the BTE in presence of surfaces, reconciling the different pictures; such discussion goes out of the scope of the present work and will be presented in an upcoming study.

VI. FURTHER PROPERTIES
A widely held assumption that is also violated by the exact BTE is the Matthiessen rule, which states that the total thermal resistivity (i.e. 1 k ) is the sum of the resistivities of each independent scattering mechanism; however, the Matthiessen rule is an approximation [18] relying on the possibility of exactly decoupling scattering mechanisms. To probe numerically this violation, we computed the resistivities of normal, Umklapp and isotopic processes, or any combination of these, and combined them according to Matthiessen rule. In Fig. 9 we show that, regardless of any particular decomposition, the conductivity obtained by imposing the Matthiessen sum deviates significantly from the exact conductivity.
. Due to the correlation between scattering events, there is no decomposition for which the Matthiessen rule is obeyed at all temperatures. The only curve that approaches the exact thermal conductivity is the decomposition (k N +U + k I ), just when the effect of isotopic scattering is negligible.
The only case in which a decomposition reproduces the exact result is when the effect of a separated mechanism is negligible; for the case shown in figure, one can sum separately the resistance due to isotopes only at high temperatures, when it's small. Finally, one can prove that the total thermal conductivity is always smaller or equal to the Matthiessen sum (see [18] or Appendix C); this is also verified in our calculations. Moreover, it is not possible to distinguish the contribution to a relaxon relaxation time due to each scattering mechanism, at variance with the case of a phonon lifetime in the SMA. This is because it would require that the eigenvalues of the sum of two matrices were the sum of eigenvalues of two matricesclearly not the case when the scattering matrix is not diagonal.
As an added benefit, the direct diagonalisation of the scattering matrix brings clear insight into the numerical stability of current methods used to solve the BTE. In particular, we show in Appendix D that the iterative method [4,16,17], often used to study 2D materials, is numerically unstable for graphene at room temperature, due to the dominant contribution of the out-of-diagonal terms in the scattering matrix (this is exactly the case when the relaxon picture differs significantly from the phonon picture).

VII. CONCLUSIONS
In summary, we have shown that by choosing the eigenvectors of the scattering matrix as a basis, the linear BTE can be greatly simplified. These eigenvectors are collective excitations of phonon populations, termed relaxons, that are characterised by well-defined relaxation times and, in the homogeneous case, also by proper velocities and mean free paths. Thermal transport can thus be described as a kinetic theory of a gas of relaxons. The characterisation of relaxon properties provides a description of the thermal transport in terms of proper time scales, and in the steady-state homogeneous case of velocity and length scales. For clarity, we report in Table I a summary of relaxon characteristics and how they compare with phonons. This theory is applied here first to graphene at room temperature where, as is typical of 2D materials or of 3D solids at low temperatures, the failure of the SMA and of its picture of phonons as heat carriers becomes dramatic; and to silicon at room temperature, where, although the SMA yields reasonable thermal conductivities, the theory brings new insight in the microscopic interpretation of heat flux and its typical velocities. Finally, we have shown that the Matthiessen rule is violated in the exact BTE, with significant consequences for all systems in which the SMA does not hold. As a final remark, the concept of relaxons has been applied in this work in the context of phonons; however, similar arguments will hold for the electron BTE or other semiclassical transport models.

First-principles simulations
Density-functional theory calculations have been performed with the Quantum-ESPRESSO distribution [45], using the local-density approximation and normconserving pseudopotentials from the PSLibrary [53]; for graphene a plane-wave cutoff of 90 Ry and a Methfessel-Paxton smearing of 0.02 Ry have been used and for silicon a plane-wave cutoff of 100 Ry. Graphene is simulated with a slab geometry, using an optimized lattice parameter a = 4.607 Bohr and a cell height c = 3a; for silicon we find an optimized lattice parameter of 10.18 Bohr. The Brillouin zone is integrated with a Gammacentered Monkhorst-Pack mesh of 24×24×1 points for graphene and 12×12×12 for silicon. Second and third order force constants are computed on meshes of 16×16×1 and 4×4×1 points respectively for graphene and 8×8×8 and 4×4×4 for silicon, and are later Fourier-interpolated on finer meshes.

Thermal conductivity simulations
The scattering matrixΩ includes 3-phonon interactions and harmonic isotopic scattering [6,7] at natural abundances [54] (98.93% 12 C, 1.07% 14 C for carbon, and 92.22% 28 Si, 4.67% 29 Si, 3.09% 30 Si for silicon). For graphene, the scattering matrix is constructed using the same computational parameters of Ref. [29] (a Gaussian smearing of 10 cm −1 and a mesh of 128×128×1 points for integrating the Brillouin zone), resulting in a matrix of order 98304, while for silicon we use a Gaussian smearing of 7 cm −1 and a mesh of 30×30×30, yielding a matrix of order 162000.Ω is diagonalised exactly using the routine PDSYEV of the Scalapack library [55]. The simulation cell of graphene is renormalized using the interlayer distance of bulk graphite (c/a = 1.367), in order to have a thermal conductivity comparable with the 3D counterpart. We verified the correctness of the software implementation ensuring that the thermal conductivity estimated with the diagonalization solver coincides with that computed with the variational method of Ref. [7] up to at least 4 significant digits. It's worth mentioning that these calculations are not prohibitively expensive and could be extended to other systems. The present software implementation computesΩ, diagonalizes it and computes the conductivity of graphene in about 5 hours using 256 CPUs on the Piz Daint supercomputer of the Swiss National Supercomputer Center (CSCS), for a total of 1300 CPU hours (1000 of which are spent in the diagonalization). For silicon the calculation completed in 8 hours on 576 CPUs, for a total of 4600 CPU hours. Calculations have been managed using the AiiDA materials' informatics platform [56] APPENDIX A: SCATTERING RATES In this appendix we report the expressions for building the scattering matrix using 3-phonon and isotope scattering events, which are discussed in detail in Ref. [7].
To make a connection with other studies, we note that most contemporary ones have largely preferred to solve the BTE using a phonon deviation from equilibrium of the form n µ =n µ +n µ (n µ + 1)F µ . Since the action of the collision operator must not change, we have the relation: where A is the scattering matrix when it acts on F , related with the scattering matrices used in our work by Ω µµ = 1 n µ (n µ + 1) A µµ 1 n µ (n µ + 1) .
The scattering rate for a phonon coalescence event is: where N is the number of q points, G is a reciprocal lattice vector and V (3) is the third derivative of the unit cell energy E cell with respect to atomic displacements with where b is an index running on the basis of atoms in the unit cell, R l is a Bravais lattice vector identifying the l th unit cell inside the crystal, α is a cartesian index, M b is the mass of atom b, z is the phonon polarization vector and u is the vector of atomic displacements. The scattering rate for a phonon-isotope scattering event is: where Combining these scattering rates, the scattering matrix A is: In the numerical implementation, the Dirac delta conserving the energy is replaced by a Gaussian smearing.
As the authors of Ref. [7] noted, the above expression guarantees that the scattering matrix A is symmetric and positive-definite also in presence of a Gaussian smearing (other expressions, which would be equivalent with a Dirac delta function, may introduce spurious negative eigenvalues). By virtue of Eqs. (16) and (17), it follows thatΩ is symmetric but not Ω, hence the necessity of the transformations (5) and (6). Finally, we recall that phonon lifetimes are related to the diagonal elements of the scattering matrix as:

APPENDIX B: SILICON THERMAL PROPERTIES
In this appendix we study the effect of the scaling of specific heat on velocities and mean free paths. While relaxon properties are as defined in the main text, phonon velocities and mean free paths are scaled as: With this choice of scaling, we can write the SMA thermal conductivity as: which treats specific heat in the same way as Eq. 13. In Figure 10 we report the comparison of these scaled phonon quantities with the relaxon properties in silicon.
One can readily see that the 2 orders of magnitude of difference that appeared in Figs. 7 and 8 have almost disappeared. Most of the discrepancy is thus due to the usage of specific heat. Nevertheless, the largest phonon velocities are still a factor of 3 smaller than those of relaxons, and the two pictures do not perfectly overlap.

APPENDIX C: MATTHIESSEN RULE
Here we recall a known result [18] that proves that the application of the Matthiessen rule results in an overestimation of the exact thermal conductivity. The BTE for a homogeneous system under a static gradient of temperature can be written in a matrix form (see for example Ref. [18] or more recently Ref. [7]): where A is related to the scattering matrix Ω via A µµ = Ω µµ n µ (n µ +1), b = − ∂nµ ∂T v µ and φ is the deviation from equilibrium defined as n µ =n µ +n µ (n µ + 1)∇T φ µ .
Another way of solving the BTE, besides the diagonalization approach discussed in the main article, is via the variational principle [18]. In particular, the solution of the BTE can be found from the minimization of the functional [18] Let φ be the function minimising F. The minimum of F is directly proportional to the thermal resistivity ρ [18]; therefore we write Now, let's separate the scattering matrix into two different components (for example 3-phonon and isotopic scattering) The exact resistivity is given by: The function φ that minimises the functional defined by A will not be, in general, the function that minimises the functionals F 1 and F 2 defined by A 1 or A 2 only. The functionals F 1 and F 2 instead will be minimised by the functions φ 1 and φ 2 respectively. By the variational Alternatively, this can be written as showing that the Matthiessen rule is a special case where the equalities holds exactly. More generally, its application leads to an overestimation of thermal conductivities.

APPENDIX D: ITERATIVE METHOD
In this appendix we examine the convergence properties of the iterative method for solving the BTE. Such method can be formalised as follows. The steady state homogeneous BTE in the phonon basis is: This is a linear algebra problem of the form AF = b, where b µ = −v µ ω µnµ (n µ + 1), F is defined by n µ = n µ +n µ (n µ + 1)∇T F µ and Ω µµ = A µµ n µ (n µ + 1). The iterative solution for F [4,16,17] can then be recast [7] as a geometric series F = where A d and A od are respectively the diagonal and the off-diagonal parts of A. This series is convergent if and only if all the eigenvalues λ of (A d ) −1 A od are |λ| < 1. In Fig. 11, we show that in graphene |λ| > 1 for more than half of the spectrum, proving that the iterative method is numerically unstable for graphene at room temperature. In general, one might expect convergence issues for the iterative method whenever the relaxon picture differs significantly from the phonon picture and the contribution of the off-diagonal part is large compared to the diagonal part. Eigenvalues λ of the matrix A −1 d A od (see main text for the definition), for graphene at room temperature, ordered by their magnitude. The red dots, roughly half of the eigenvalue spectrum, indicate eigenvalues |λ| > 1 that cause a divergence of the iterative solution for the Boltzmann transport equation. Most of the unstable eigenvalues are greater than 1, with only one eigenvalue lower than -1.