Green function, quasi-classical Langevin and Kubo–Greenwood methods in quantum thermal transport

With the advances in fabrication of materials with feature sizes at the order of nanometers, it has been possible to alter their thermal transport properties dramatically. Miniaturization of device size increases the power density in general, hence faster electronics require better thermal transport, whereas better thermoelectric applications require the opposite. Such diverse needs bring new challenges for material design. Shrinkage of length scales has also changed the experimental and theoretical methods to study thermal transport. Unsurprisingly, novel approaches have emerged to control phonon flow. Besides, ever increasing computational power is another driving force for developing new computational methods. In this review, we discuss three methods developed for computing vibrational thermal transport properties of nano-structured systems, namely Green function, quasi-classical Langevin, and Kubo–Green methods. The Green function methods are explained using both nonequilibrium expressions and the Landauer-type formula. The partitioning scheme, decimation techniques and surface Green functions are reviewed, and a simple model for reservoir Green functions is shown. The expressions for the Kubo–Greenwood method are derived, and Lanczos tridiagonalization, continued fraction and Chebyshev polynomial expansion methods are discussed. Additionally, the quasi-classical Langevin approach, which is useful for incorporating phonon–phonon and other scatterings is summarized.


List of symbols
Symbol Explanation

Introduction
Thermal transport plays a central role in most of the state-ofthe-art technological applications. In electronics, one requires good thermal conduction in order to get rid of the excess heat because the device temperature limits the operation speed.
On the other hand, there are also cases where a poor thermal conduction is preferred. Thermoelectrics is a typical example, where one uses temperature difference to drive an electric current with minimal exchange of heat between the reservoirs or vice versa. At the macroscopic scale, Fourier's Law gives an accurate description of thermal transport, which states that the thermal current is opposite and proportional to the temperature gradient and the proportionality constant is thermal conductivity, i.e. J = −σ∇T . At the nano-scale, where the device lengths are less than or comparable to the phonon mean-free-path, ballistic transport regime becomes relevant and the diffusion picture fails. In fact, the entire phonon spectrum is involved in thermal transport, in general, therefore a wide range of length scales are relevant for phonon transport. This could be at the order of µm for acoustic modes, and as short as a few nm for high energy optical modes with small group velocities. The relation between length scales and thermal transport is important not only for understanding the effects of nano-structuring on phonon transport, but also because nano-structuring has the potential to bring novel concepts and applications involving phonons. One such concept of fundamental nature is the quantum of thermal conductance, which is observed in suspended insulating nano-structures at very low temperatures [1][2][3][4]. Moreover, there have been considerable efforts to build phononic analogs of electronic devices. The emerging field of phononics aims to use phonons to carry and process information with thermal diodes, thermal transistors, thermal logic gates and thermal memories [5].
Parallel to these developments, theory and modeling of phonon transport has been an active field of research in the recent years [6][7][8][9]. Classical molecular dynamics (MD) simulations, Boltzmann formalism, equilibrium and nonequilibrium Green function (GF) methods, generalized quasi-classical Langevin (QCL) approach, master equation formalism and the realspace Kubo-Greenwood (KG) methodology are among the widely used methodologies to compute phonon transport at the nano-scale. MD methods are based on the Verlet algorithm, that is integrating Newton's equations of motion. This is a disadvantage if combined with Boltzmann statistics neglecting quantum mechanics entirely. In graphene, for example, Debye temperature is 2300 K for in-plane, and 1300 K for out-ofplane modes [10]. Methods which are based on quantum statistical mechanics, e.g. GF, QCL, and KG do not suffer from this problem. On the other hand, MD methods include anharmonicity intrinsically and to all orders, whereas for GF methods one has to include anharmonic interactions separately. At distances shorther than the anharmonic mean-free-path, GF is quite accurate in predicting scatterings due to defects, interfaces, grain boundaries, and nano-structuring schemes. A way of include phonon-phonon scatterings in those systems is to use generalized QCL approach. Apart from these, accurate treatments of anharmonic scatterings has been implemented based on density functional theory and Boltzmann formalism [11]. With this approach, one can obtain excellent agreement between the calculated and measured intrinsic lattice thermal conductivities of semiconductors. Yet, the method is more suited for periodic systems without disorder because of the computational cost.
Here, a review of the GF, QCL and KG methodologies is presented so as to introduce the reader with the fundamentals and applications of these approaches. It is neither our intention nor it is possible to cover these topics in full detail in a single review paper. There are excellent reviews on the very details of the GF methodology [12][13][14][15] and the QCL approach [16] in the literature. Our goal in this review is to guide the reader directly to the transport equations. Therefore some of the important details of the GF methodology, like the equation of motion on contour, Feynmann diagrammatic perturbation theory are excluded. The refer the reader to the aforementioned reviews for further details.
The paper is organized as follows. In section 2, the GF methodology is reviewed, where thermal current expression is derived using nonequilibrium Green functions (NEGF) and the equilibrium expression is obtained, which is of Landauer type. Decimation algorithms and surface GFs for reservoirs are also discussed in detail. The examples of the use of GF formalism include defected carbon nanotubes, graphene with molecular antiresonances, a hybrid nanostructuring scheme for reducing thermal conductance, graphene grain boundaries, anistropy of thermal transport in puckered structures and strain Gt Anti-time-ordered Green's function g µ,r Free GF of the reservoir µ. µ ∈ {L, R} Σ L/R Self energy due to coupling to the left/right reservoir Σ L/R = Φ CL/CR g L/R Φ LC/RC Γ L/R Broadening due to coupling to the left/right reservoir Transmission spectrum f B (ω, T) Bose-Einstein distribution function κ (T) Conductance J Thermal current σ Thermal conductivity elastic Elastic mean-free-path ξ Driving force ξ µ (t) = Φ Cµ g µ (t)u µ (0) S µ Gaussian noise S µ = dt e iωt ξ µ (t)ξ µ (0) engineering of thermal transport properties. A simple model for the reservoirs and a force constant model for obtaining the dynamical matrices for graphene based systems are also included, which enable an easier starting with the reviewed methods. In section 3 the QCL approach is discussed, which bridges GF and MD approaches. Section 4 is reserved for the real-space KG approach, which is a linearly scaling method for disordered systems. Disordered carbon nanotubes, edge disordered graphene nanoribbons and graphene/BN heterostructures are the examples reviewed using this method.

Green function method
Green functions are being used in almost all branches of physics. It finds its roots in quantum field theory and has been used extensively in many-particle physics [17]. The NEGF formalism was formulated in the pioneering works of Schwinger [18], Keldysh [19], Kadanoff and Baym [20]. NEGF was used by Meir and Wingreen [21] to study electronic transport, and later used to address phonon transport [22][23][24].

The Hamiltonian
The total energy of the system E tot can be expanded in terms of the displacements of the constituent atoms from their equilibrium positions. Denoting m i and x iα as the mass of ith atom and its displacement along α-direction (α = x, y, z), one defines mass-normalized displacements as u iα = √ m i x iα and writes the Hamiltonian as where Φ ij,αβ = ∂ 2 E tot /∂u iα ∂u jβ are the elements of the dynamical matrix. Ψ ijk,αβγ = ∂ 3 E tot /∂u iα ∂u jβ ∂u kγ and Ψ ijkl,αβγθ = ∂ 4 E tot /∂u iα ∂u jβ ∂u kγ ∂u lθ are the mass-normalized anharmonic force constants up to fourth order. The third and the higher order terms belong to the non-linear part of the Hamiltonain, H int , whereas the harmonic part ( H harm ) defines an eigenvalue problem.
One can define |u as the normalized column vector of displacements u iα , u| as its transpose such that u|u = 1, and denote the eigenvalue problem defined by the dynamical matrix, Φ, as with u q being the eigenvector and ω 2 q is the corre sponding eigenvalue with index q. The equivalent form without mass-normalization is a generalized eigenvalue problem, ij,αβ = ∂ 2 E tot /∂x iα ∂x jβ are the harmonic force constants and M is the diagonal matrix of atomic masses.
The annihilation and creation operators are defined as as which satisfy the commutation relation [a q , a + q ] = δ qq , and diagonalize the harmonic Hamiltonian, H harm = q ω q (a + q a q + 1 2 ). Also note that [u iα ,u jβ ] = i δ ij δ αβ .

Six phonon Green functions
In equilibrium, it is enough to define one Green function, which is generally chosen as the retarded GF. In nonequilibrium cases, it is common practice to define six Green functions, which are not linearly independent. Their definitions are [17], which will be referred to as the greater, lesser, time-ordered, anti-time-ordered, retarded, and advanced Green functions, respectively. The following relations can be derived from the above definitions [17].  The contour C, on which t 1 is before t 2 , is deformed as 2.2.1. Contour ordering. In a non-equilibrium situation, the system does not return to its initial state at asymptotically large times, in general. In such cases, contour ordering is the suitable method for calculating correlations. Contour ordered GFs are non-equilibrium analogs of time-ordered GFs. Generally, the initial correlations are neglected and the contour is chosen like in figure 1.
In GF calculations, one evaluates terms like In order to calculate C < , one chooses a contour C such that t 1 is before t 2 (see figure 2(a)). Then the contour is deformed as C → C 1 + C 2 (see figure 2(b)), so that t 1 is on C 1 and t 2 is on C 2 . Using these, one obtains the expression for C < as Here, one has B < in the first term because all τ on C 1 are before t 2 , which is on C 2 . Similarly, one has A < in the second term because t 1 is always before τ , which is on C 2 . Now, we do the the first integral. On C 1 , τ can be before or after t 1 , which is expressed as Likewise, the C 2 integral takes the form (11) and (12)), one obtains It is possible to use a matrix representation for the GFs and the self energies as [25] Using the matrix notation, Ã =BC, one can easily write (9) and (10)), one again obtains A < = B r C < + B < C a . It will be necessary to calculate such correlations during the course of phonon transport and equation (21) will be used in deriving the expression for the heat current. Adiabatic switching is employed to study the non-equilibrium and interacting systems. It is assumed that the subsystems (reservoirs and the device region) are decoupled and the nonlinear interactions are turned off at t → −∞. The decoupled subsystems are at temperature T α and the lesser and greater GFs satisfy ). (24) The couplings between the subsystems (Φ LC/CL and Φ RC/CR , see the next section) are turned on slowly so that steady state is established at t 0. Later the nonlinear interactions are turned on at t = 0.
Applying perturbation expansion, Dyson's equation with contour ordering can be obtained as and using the Langreth rules one obtains with Σ ≶ originating from couplings between subsystems.

Partioning scheme
In the GF and QCL approaches the system is partitioned (figure 3). The most common case is to include two reservoirs, which are labeled as the left (L) and the right (R) reservoirs and assumed to be kept at different temperatures. In this review T L > T R , if the opposite is not stated. The reservoirs are assumed to be perfect, namely no scatterings take place inside them. The central region (C) connects the reservoirs, which are disconnected otherwise. The boundaries of the central region are chosen such that all scattering events, whether elastic or inelastic, are to be included inside the central region.
In other words, convenient parts of the reservoirs are included in the central region so as to make the force constant matrices of the reservoirs perfectly periodic in one direction. The partitioning scheme applies to the Hamiltonian as The retarded Green function, also referred to as the resolvent, can be written as where δ is an infinitesimal positive number and is the unit matrix. The elements of the resolvent read G r qq = δ qq /(ω 2 − ω 2 q + iδ) in the diagonal basis, δ qq representing Kronecker-δ function.
Using the partitioning scheme, the dynamical matrix and the resolvent can be written as 3 × 3 matrices, and equation (30) reads (31) The second column of the product matrix consists of the following equalities.
where g L/R,r = [(ω + iδ) 2 − Φ LL/RR ] −1 are the retarded Green functions for uncoupled reservoirs. Using equations (32)-(34), Green function for the central region is obtained as with Σ L/R = Φ CL/CR g L/C Φ LC/RC being the self energies due to coupling to the reservoirs. The broadening matrices are defined in terms of the self energies as

Expression for the heat current and Landauer-type formula
At steady state, the heat current can be defined in terms of the rate of energy change at one of the reservoirs. That is, Here, the indices l 1 , c 1 etc include both the atom index and the Cartesian components (see equation (6)). Performing contour ordering within the partition scheme, one has which leads to the expression Applying Fourier transformation and symmetrizing the real valued function as J = (J + J * )/2, one finds (43) Using equation (23) for the uncoupled GFs, one can write Similarly, using equation (24), one gets Employing Dyson's equation (see equation (26)) for G CC,< and G CC,> as This equality reproduces the Landauer-type formula given Rego and Kirczenow [2]. (49) where T L > T R are reservoir temperatures and ζ(ω) is the transmission spectrum. The expression for the transmission function was obtained by Yamamoto and Watanabe [24] as Equivalent formulas for the transmission function were previously derived by Özpineci et al [22] and Mingo et al [23] In the limit of small temperature difference, namely The function p(x) is equal to unity at x = 0, it decays rapidly for x>1 and ∞ 0 dx p(x) = π 2 /3. Therefore at the low temperature limit (T→0), the thermal conductance is quantized as is the temperature dependent quantum of thermal conductance and M is the number of acoustic modes [1-3].

Generalization to multiple reservoirs
The Landauer-type expression can be generalized to multiterminal configurations, in which case the central region is connected to more than two reservoirs. Each reservoir µ is assumed to be kept at temperature T µ . Transmission between reservoirs µ and ν is where Γ µ is the mode broadening due to coupling to reservoir µ, Γ µ = Im [Σ µ,a − Σ µ,r ] (see equation (36)). G CC collects self energy contributions from all reservoirs, The multi-terminal systems are of particular interest for applications such as thermal rectification [26][27][28][29][30][31].

Self-consistent reservoir approach
A method which is based on the multiple reservoirs is the selfconsistent reservoir (SCR) approach. The method is originally developed for a classical harmonic chain consisting of N atoms [32][33][34]. Each atom is coupled to an independent reservoir (see figure 4), which acts like a Büttiker probe [35,36]. Temperatures at the end points are set to T 1 and T N , which correspond to the temper atures of the left and right reservoirs, T 1 = T L and T N = T R . The rest of the temperatures from T 2 to T N−1 are to be determined in a self-consistent manner by setting the net heat current equal to zero for each reservoir at steady state. The classical model is able to reproduce diffusion dynamics and Fourier's law [34]. The model is also extended to the quantum domain [37,38], which has been solved at the linear response regime [39][40][41], approximately in an analytical way [42,43], and numerically exactly [44,45]. It was shown by Bandyopadhyay et al that, at high temperatures quantum mechanical self-consistent model reproduces the classical results, whereas at low temperatures classical simulations overestimate the currents [45]. At large bias and for low temperatures, the quantum model can also reveal thermal rectification for asymmetric chains, which is absent in the classical model [45]. The SCR approach is also applied to realistic systems beyond toy model calculations. Sääskilahti et al computed the temperature profile in a grpahene nano-constriction using the SCR approach, where the temperatures of the left and right reservoirs are 300 K and 280 K [46]. The computed temperature profiles are found to be in agreement when quantum exact and classical approximate approaches are used, which is interpreted as transport being mainly due to the low frequency acoustic modes.

Decimation algorithms for efficient computation of Green functions
Decimation is a numerically exact method for reducing the sizes of block tridiagonal matrices. Consider a tripartite dynamical dynamical matrix of the form Self-consistent reservoir approach. Each atom is coupled to a self-consistent reservoir. Temperatures of the end reservoirs are kept fixed while the others' are determined self-consistently with the demand that no net current flows from them.
The retarded GF is defined as Let us obtain an effective dynamical matrix by decimating the middle layer. To do, we first need to define the free GF for the middle layer from which the self energy terms arise as and the 2×2 effective dynamical matrix is obtained as where Φ eff is a function of ω.

Recursive algorithm.
It is straightforward to generalize the above scheme for an extended tridiagonal matrix, where i, j of the submatirces Φ ij are the layer indices and the layer thickness is to be chosen such that Φ ij = 0 if |i−j|>1. One can start decimating from the second layer, where the effective dynamical matrix parts and the corresponding G (n) of the (n+1)st recursion step are Repeating the decimation N−2 times, one arrives at the 2×2 effective GF for the central region consisting of G 1,1 , G 1,N , G N,1 and G N,N .

Renormalization-decimation algorithm.
In GF computations, the reservoir regions are considered to be perfect, in the sense that they are free of scatterings. Therefore the reservoirs are made up of identical layers. That is, the dynamical matrix consists of identical blocks. For such periodic systems, one has Efficient algorithms have been devised to compute the bulk and surface Green functions of such systems [47,48]. Decimating each second layer, the effective bulk dynamical matrix becomes Repeating this procedure by using which converge very fast compared to the standard decimation method and incorporation of a thousands of layers is possible within a few steps.

Surface Green functions (SGF)
For semi-infinite reservoirs, dynamical matrices have the following form.
and the surface terms are obtained as The convergence of the iteration is checked by comparing the difference between dynamical matrices at two consecutive steps.

Generic reservoir as a simple model for SGFs.
Atomistic quantum transport calculations are quite expensive computationally. One of the most time consuming parts is the computation of SGFs. A widely used approximation in calculating electronic transport across nano-junctions is the wide band limit. The wide band limit produces fairly good results, especially when the reservoir density of states (DOS) is approximately constant around the Fermi energy. (see e.g. [49]) Recalling that the DOS is proportional to the imaginary part of the Green function, and that the real part is its Hilbert transform, the wide band limit introduces ' structureless' reservoirs. Namely, during its injection into the central region, the information carried by an electron regarding the electronic structure of the reservoir is the Fermi velocity. Such a wide band limit approximation is not as sensible for phonons as it is for electrons, because at ambient temperatures phonon transport involves a substantial range of the phonon spectrum. The phonon modes have wide range of group velocities. For acoustic modes, the modes at the zone center have the speed of sound while at the zone boundaries their velocities are zero. Still it is possible to define structureless reservoirs, which are physical in the sense that they sample the group velocities that are correctly matched with the wave vectors inside the Brillouin zone, and the acoustic phonons having the desired sound velocity. Keeping in mind that it is the acoustic phonons which carry most of the heat, such structureless reservoirs should result in fairly good agreement with more realistic SGF computations. Indeed, when phonon transport is computed across molecular junctions consisting of carbon nanotubes and azobenzene molecules, thermal conductance values have been found to be in fairly good agreement with those obtained from the 'generic reservoir model', which is explained below [50].
Consider a semi-infinite portion of a chain of identical atoms having mass m and harmonic spring constant It was shown by Müller et al [51] that the elements of the GF are given as where θ = arccos(1 − ω 2 /2ω 2 0 ), and {m, n} ∈ [1, ∞ ) with m = n = 1 corresponding to the on-site term for the end atom. Therefore the SGF is obtained as with η = ω/ω 0 . The real and imaginary parts of the SGF are plotted in figure 5. We emphasize that the dynamical matrix in equation (79) is not the DM of a semi-infinite system but a semi-infinite portion of a periodic system. The difference is in the Φ 11 element of the DM, which would be equal to ω 2 0 for a semi-infinite system. We take it equal to the remaining of the diagonal elements by implicitly assuming that a finite portion of the reservoir atoms is included in the central region to ensure that the diagonal elements of Φ are the same. If one end was free, i.e. if Φ 11 = ω 2 0 , then Φ eff 11 = ω 2 0 (η 2 /2 − iη 1 − η 2 /4) and G 11 = (ω 2 − Φ eff 11 ) −1 and the surface density of states  The radial components (f r ) are along the x-axis, while the in-plane transverse components (t ti ) and the out-of-plane transverse components (f to ) are parallel to the y -and z-axis, respectively.
That is, the density of states is finite at ω = 0 and decreases monotonically with increasing frequency.

A minimal force constant model for graphene based materials
DFT based methods yield reliable force constant matrices which predict phonon dispersions in quite good agreement with available experimental data. On the other hand, for complex systems of considerable sizes, one requires computationally feasible and reasonably accurate methods for determining the force constant matrices. The force constant model originally proposed by Saito et al for graphene based structures suits this purpose very well [52]. The model was implemented to various cases including carbon nanotubes, graphene nanoribbons, defected graphene structures, etc [26,[53][54][55][56][57][58][59][60][61][62][63][64]. Also, a reoptimized version of the force constant parametrization is available [65]. The method relies on parametrization of interatomic force constants up to fourth nearest neighbors (4NNFC approach). For each neighboring distance (n = 1, . . . , 4) three force constant parameters f to are defined, which stand for the radial (bond stretching), in-plane and out-of-plane tangential (bond bending) directions. The diagonal form of the force constant matrix, when both atoms lie on the x-axis (see figure 6), reads When the direction from a given atom to its neighbor is rotated around the z-axis by θ p to account for atom Bp, the force constant tensor is rotated as where For carbon nanotubes, that is when the atoms do not lie on the xy-plane, the tensor should be rotated around the axis of the tube (y -axis). For A type atoms, the rotation angle ϕ i is the polar angle between A1 and Ai around the circumference and , n a being the number of atoms per unit cell. For B type atoms, the tensors shuld be rotated around z by π, The parameters, as optimized by Zimmermann et al are given in table 1. The phonon dispersions for graphene with these parameters is compared against DFT results in figure 7, which display excellent agreement for the low energy acoustic modes and a fairly good agreement for the optical modes. The quadratic out-of-plane acoustic ZA mode is of particular importance, and its dispersion is reproduced perfectly with 4NNFC.

2.10.Applications of the Green function method
The advantage of GF methods over others is that, having obtained the GF for the region of interest, one has access to all physical quantities decomposed to their frequency components. This is of particular importance for phonons, where for a given device length different transport regimes take place at different parts of the spectrum, all of which contribute to heat current. GF method has been used for a range of systems, some of which are atomic chains [22], molecular junctions [66][67][68][69] carbon nanotubes and graphene based materials [63,[70][71][72][73][74][75][76][77][78], silicon nanowires [79][80][81], novel 2D crystals [82,83]. Below we review some of these works, where the advantageous sides of the GF methodology has been useful. In section 2.10.1, a complementation of MD and GF methods for defected carbon nanotubes is reviewed. Section 2.10.2 considers a controlling scheme which incoporates antiresonances in the transmission spectrum and section 2.10.3 examplifies the use of GFs to analyse the combined effects in hybrid nano-structuring, section 2.10.4 investigates the effects of molecular functionalization of grain boundaries on phonon scattering. Sections 2.10.5 and 2.10.6 are reserved for anistropic transport and strian engineering of phonons.

Defected carbon nanotubes. Carbon nanotubes (CNTs)
have been the subject of keen interest since the earliest days of nanotechnology research. Not only because it is an ideal example for a quasi 1D material, but also due to its unprecedented physical properties like extreme mechanical strength [84], ballistic electron transmission over long distances [85], and record thermal conductivity values [86]. Despite the fact that thermal conductivities of CNTs are loosely dependent on their diameter and chirality, a large range of values have been reported, displaying a large variance [87][88][89][90][91][92]. The effects of different defect types on armchair and zigzag CNTs of various lengths were studied using nonequilibrium molecular dynamics and it was found that small changes in the defect concentration give rise to large changes in thermal conductivities, when the defect concentration is low ( figure 8). On the other hand, at higher concentrations, thermal conductivity becomes insensitive to changes in concentration [71]. Another interesting issue regards divergence of conductivity with length. Thermal conductivity values as obtained from NEMD simulations keep increasing with length and no convergent value was observed even for systems as long as 1 µm. (see figure 9) This rather odd property was attributed to very long acoustic mean-free-paths.
Both properties, namely the peculiar dependence on defect concentration and seemingly divergent behavior of conductivity with length were possible to explain through an analysis using the GF formalism within the 4NNFC approach. Relying on the fact that the thermal conductivities weakly depend on defect type, randomly distributed monovacancy defects were used. Defect concentrations (d = N defect /N atom ×100) ranging between 0.1% to 1.5% were considered for CNTs as long as 4 µm. Transmission spectrum for d = 1% is shown as a function of frequency for various lengths in figure 10 (left). The zeroth of the length axis stands for ballistic transport, where the transmission changes take place stepwise. Transmission decays rapidly with length for ω>600 cm −1 . This is also evident from the mean-free-paths as shown in the right panel. At lower frequencies (ω < 300 cm −1 ) transmission is considerably larger at longer distances. For low energy acoustic modes (ω < 150 cm −1 ) transmission stays almost intact at 4 µm of length. This is in agreement with the finding that elastic is much longer than a µm for the first transmission plateau of acoustic phonons.
The transmission and elastic data in figure 10 sheds light on the NEMD results. Conductivity drops rapidly with defect concentration because elastic is reduced by almost an order of magnitude with concentration. More importantly, the drop is  The discrepancy in the experimental data for the thermal conductivity of CNTs could stem from the fact that the deviations are large at low defect concentrations. For a fixed value of variance in defect concentration, if the mean is shifted from 0.05 to 0.5 deviation in conductivity is lowered from 92 W m −1 K −1 to 4 W m −1 K −1 . Therefore, increased d is expected to enable a standardized value for the conductivity, which is still higher than 100 W m −1 K −1 . The GF calculations for CNTs with vacancies have shown that the low frequency acoustic modes have mean-free-paths much longer than sizes which can fit to any MD simulation. The phonon-phonon scatterings do not affect this behavior because of ω −2 dependence of the anharmonic mean-freepath [93].

2.10.2.
Graphene with molecular antiresonances. The large surface to volume ratio of 2D materials enable very effective engineering schemes for physical properties by surface functinalizations [94][95][96]. One possible route is to adsorbe molecules covalently. Each adsorbed molecule acts as a scatterer for phonons by creating tranmission antiresonances [97,98], therefore is a source of thermal resistance. A fundamental difference of phonon transport from electron transport is that, it takes place at the entire spectrum unlike electronic transport, which involves a small energy window around the Fermi level.  This gives rise to a striking difference in resistance summation rule, as it was shown in [78]. First, phonon transmission spectra of graphene with single molecular scatterers of different species were computed, where the force constants were obtained using the density functional tight binding (DFTB) approach [99]. CH 3 , C 2 H, C 3 H 3 , C 4 H, C 5 H 3 , benzene, biphenyl and para-terphenyl were used as adsorbates. It was shown that antiresonances in the transmission occur, which are due to vibrational modes that are localized on the molecules. These localized modes are also observed in the DOS as isolated peaks (see figure 11(a)-inset.) Besides the differences in their detailed structure, the phonon transmission spectra follow the same trend in the whole spectrum, independent of the molecular species ( figure 11(a)). Consequently, thermal conductance with a single adsorbant is also similar for all species ( figure 11(b)). Transmission due to multiple scatterers was calculated at the cascade scattering approximation [100,101] using where ζ o is the transmission of the pristine system, ζ i is transmission due to adsorbant species i, and n i is the number of those molecules. Keeping the total number of molecules constant and equal to 100, the change in n i is made in multiples of 10. While changing the number of species N sp from 1 to 8, transmission and conductance values were computed for every possible combination of n i . It was found that the average conductance values decrease as the number of species is increased at all temperatures (see figure 12). This result proves that thermal resistances are not additive, in general. It was further concluded that it is possible to obtain desired thermal transport properties by using appropriate combinations of molecular species. The key ingredient of this scheme is to obtain antiresonances in the transmission spectrum. Having distinct antiresonaces suppresses conduction more effectively. Antiresonances at lower frequencies are more effective to suppress conductivity, whereas strong bonding results in wider antiresonances.

2.10.3.A hybrid nano-structuring scheme using geometri-
cal structuring and isotopical disorder. Thermoelectric efficiency is a challenging topic because it requires diverse physical properties. Specifically, one needs not only a good electrical conductance and high Seebeck coefficient, but also a low phonon thermal conduction. Meeting those requirements in a single system could be achieved by using hybrid nanostructuring strategies. One of the proposals is using geometrical structuring and random distribution of clustered heavy isotopes, which suppress thermal conduction and increase thermoelectric figure of merit when combined [74].
Two types of graphene nanoribbons (GNRs) were considered, one having a straight geometry (s-GNR), whereas the other is chevron-type (c-GNR), both of which were synthesized via bottom-up methods (figure 13) [102]. Using the GF method with the DFTB approach, it was found that the pristine conductance values of s-GNR and c-CNR are quite different (5.0 and 2.3 nW K −1 nm −2 at 300 K, respectively), which is only due to difference in geometrical structure. Computing phonon transport using GF technique for ensembles of isotopically disordered structures, the authors reported that conductance of s-GNR at room temperature was reduced from 0.566 nW K −1 nm −2 to 0.460 nW K −1 nm −2 when the isotope density was the same (50%) but the distribution of isotopes was changed from atomic to precursor (clustered). The reduction is even more dramatic for c-GNR, it drops from 0.131 nW K −1 nm −2 to 0.060 nW K −1 nm −2 .
The physical mechanism of the effect is explained as follows. The distribution type effects the mean-free-paths of low and high frequency phonons differently. The precursor distribution scatters low-frequency phonons more effectively and the high-frequency phonons less effectively than the atomic distribution, a manifestation of Rayleigh scattering. The chevron geometry, on the other hand, divides the phonon dispersions into mini-bands with many singularities in the DOS. According to Fermi golden rule, scattering is proportional to the DOS. That is, the high frequency part of the spectrum is  suppressed due to flattened phonon bands. When combined, chevron geometry with precursor distribution of isotopes has an one order of magnitude less coductance than its straight counterpart with atomic distribution (see figure 14).

2.10.4.Functionalized graphene grain boundaries.
Tailoring the properties of a material primarily requires an understanding of how the material can be modified in a controlled manner. One possible strategy to influence the properties of graphene grain boundaries (GBs) can be their functionalization with molecules or ad-atoms [104,105], which can act as electronically active dopants controlling their electronic and phonon transport properties and, hence, its thermoelectric behavior. Using GF techniques combined with DFTB theory [106], Medrano-Sandonas et al [103] explored the influence of ad-atoms (Hydrogen and Oxygen) and chemisorbed molecules (Hydroxyl-OH, Methyl-CH 3 , and nitrophenyl-NO 2 C 6 H 4 (NPD)) on the thermal transport properties of two possible graphene GBs, symmetric (model I) and asymmetric (model II). Model I consists of two graphene sheets with the same orientation (8 • ) but in opposite directions, clockwise and counter clockwise, whereas model II has been built with a graphene sheet oriented in the armchair direction and another in the zigzag (10 • in the clockwise) direction (see figure 15). Both GBs are composed by pentagons and heptagons, and model II also presents a nonagon [107,108]. As shown in figure 15, the phonon thermal conductance, κ ph , for both graphene GBs is lower than the corresponding one to pristine graphene for the whole range of temperatures explored in this study (up to 800 K) because of local structural defects which induce additional phonon scattering. Moreover, when comparing κ ph for the two GBs, we see that model II has the smallest phonon thermal conductance. This is related to the more disordered structure of the GB region in model II. It has also been found that H atoms prefer to be located on top of carbon sites (C-H bond length is 1.14 Å ), while O atoms are attached to bridge (B) sites (the C-O bond length is 1.49 Å ) for both GBs, in good agreement with results reported by Mu et al [107]. The phonon transmission of these grain boundaries is reduced over the whole spectrum of frequencies by the inclusion of ad-atoms (not shown). H and O adatoms strongly affect out-of-plane modes which correspond to ω < 750 cm −1 . The analysis of the projected vibrational density of states reveals that the lighter H atoms have a stronger  influence at high frequencies (ω ∈ [1000-1500] cm −1 ) and O atoms at lower frequencies (ω < cm −1 ). The major contribution to the phonon transmission comes from a spectral region between 250 cm −1 and 750 cm −1 . Since this is also the range where oxygen functionalization shows up in phonon transport, we see that attaching oxygen atoms has a stronger influence on the thermal conductance of both GBs (see figure 15). Furthermore, OH, CH 3 , and NPD molecules are covalently bonded to carbon atoms on the grain boundaries on T sites. For both GB models, the largest influence of the functionalization is seen in the low frequency domain (ω 400 cm −1 ) as well as at higher frequencies ω > 1100 cm −1 (not shown). Besides of this reduction on the phonon transmission, thermal conductance of GB model II turns out to be less sensitive to changes in the type of molecule for T > 200 K (see figure 15(b)). The variations at low temperatures come from the strong suppression of out-of-planes modes produced by the functionalization group. For GB model I, CH 3 and NPD molecules affect in a very similar magnitude its phonon transmission which is reflected in the thermal conductance values (see figure 15(a)). Whereas OH groups show a higher decrement of the thermal conductance due to the fact that high frequency in-plane modes were suppressed after bonding. Deviations from the conductance of non-functionalized GB clearly increase with temperature as molecular higher frequency modes become thermally activated.

2.10.5.Thermal anisotropy in planar puckered struc-
tures. Motivated by the intrinsic anisotropic physical properties found in the successfully synthesized phosphorene monolayers [109,110], Medrano Sandonas et al. analyzed how the structural anisotropy affects phonon transport properties of homo-and heteroatomic 2D materials with puckered structures: arsenene and SnS monolayers (see figure 16(a)) and how do they compare with those in phosphorene [111]. Based on the DFTB phonon dispersion analysis (not shown), all studied systems were mechanically stable and did not show imaginary modes. Acoustic branches display the typical dispersion of 2D materials: the longitudinal (LA) and transverse acoustic branches have linear dispersion as the wave vector approaches the Γ point, while out-of-plane ZA branches exhibit quadratic dispersion due to the rapid decay of transversal forces. It should be noted that phonon branches for homo-atomic puckered materials are almost identical with the exception of their maximum frequency value, which is a consequence of the mass difference between As (∼75 uma) and P (∼31 uma). It has been also found that the LA branch in phosphorene shows a group velocity of 8.35 km s −1 and 4.74 km s −1 along the Γ-X (ZZ) and Γ-Y (AC) directions, respectively, which are very close to previous DFT results [110,[112][113]. The values for arsenene, ZZ-5.01 km s −1 and AC-2.71 km s −1 , are also in agreement with those reported by Zeraati et al [114] SnS monolayer gives group velocities of ZZ-6.48 km s −1 and AC-2.14 km s −1 . GF based transport calculations were carried out using periodic boundary conditions in the perpendicular direction to transport and with the same number of unit cells. The reservoirs consist of the same material as the scattering region. With this, intrinsic transport features of the different systems can be revealed. As it is expected from the phonon dispersion analysis, phosphorene and arsenene have relatively similar phonon transmission functions (not shown). The only differences are the phonon bandgaps (∼67 cm −1 and ∼45 cm −1 , respectively) and the maximum frequency mode. Likewise, phonon transmission in the ZZ direction is larger than in the AC direction for almost the whole frequency range. In particular, SnS monolayer turns out to display the strongest thermal anisotropy. As a consequence of this, the preferred direction for phonon thermal conductance of these puckered 2D materials was found to be along the ZZ direction (see figure 16(b)), [110,[113][114][115][116][117]. Moreover, at low temperatures the anisotropic effect on κ ph (measured by P AZ = κ ph−ZZ /κ ph−AC ) rapidly enhances and, after T ∼ 200 K, it keeps the same magnitude since the whole phonon spectrum has been already covered by the frequency integration. As a result of strong anisotropy in its phonon transmission function and group velocities, SnS monolayer displays the highest thermal anisotropy at 300 K, P AZ ∼ 2.5, followed by phosphorene ( P AZ ∼ 1.6), which is slightly more thermal anisotropic than arsenene ( P AZ ∼ 1.35). The P AZ value for phosphorene is slightly different from that obtained by employing DFT calculations [113,118], while for arsenene there is more discrepancy from that reported by Zeraati et al [114], which is mainly due to the methodology used. In the latter work self-consistent calculations were used to solve Boltzmann's transport equation.
2.10.6.Strain engineered 2D materials. Novel two-dimensional (2D) materials show unusual physical properties, which combined with strain engineering open up the possibility of new potential device applications in nanoelectronics. In particular, transport properties have been found to be very sensitive to applied strains. Thus, in [119] the GF method was applied to address the influence of strain engineering of the transport setup on the electron and phonon transport properties of 2D materials, focusing on hexagonal boron-nitride (hBN), phosphorene, and MoS 2 monolayers. Based on the partitioning scheme shown in figure 17(a), two possible theoretical setups were considered, which may mimic different experimental ones: (I) the device region is uniaxially strained (along ZZ direction) while the contact regions are not and (II) both the device and contact regions have the same strain level (i.e. homogeneously strained system). The same box length along the y -direction for all strain levels in x-direction (ZZ) was considered, which produces an extra force in the periodic direction after increasing the strain. The results have shown that among these three 2D materials, hBN monolayers display the highest thermal conductance κ ph at 300 K (∼310 W mK −1 ), followed by phosphorene (∼90 W mK −1 ). These values are close to those reported in other experimental and theoretical works [112,115,118,120,121]. It was also found that phosphorene shows the weakest bonds and hBN the strongest ones. Moreover, based on the strain dependence of bond lengths in the hBN monolayer (not shown), it is expected that with the exception of the B-N bond perpendicular to the transport direction, the strength of the first and second neighbor bonds will decrease after increasing the strain. Hence, the components of the dynamical matrix associated to these bonds (the strongest ones) become smaller and, then, κ ph continuously decreases with increased applied strain, as it can be seen in figure 17(b) (left panel). This reduction of κ ph is higher for setup I because of the presence of an additional interface resistance between contact and device regions, which strongly blocks phonon transfer at higher in-plane modes, >800 cm −1 . It is worth noting that by considering setup II, κ ph at room temperature (300 K) shows a slight increment. Similar results were reported by Zhu and Ertekin [120] using non-equilibrium molecular dynamics and Boltzmann transport equation methods.
In the case of phosphorene and MoS 2 monolayer the influence of the strain is qualitatively different due to the presence of additional transport channels related to high frequency out-of-plane modes, which have their origin in the additional atomic layers. κ ph in phosphorene monotonously decreases (setup I) and increases (setup II) with the strain, see figure 17(b) (center panel). The increment in strain reduces the interlayer distance and, hence, the strength of the out-of-plane bonds increases, becoming higher than the inplane ones. As a result of this, the transmission probability increases at low frequencies (in-plane modes) and so does the thermal transport, as it was also observed by Ong et al [118]. For MoS 2 monolayer, despite the relatively small change of first neighbor bonds (∼10 −2 Å ), their strength was found to increase considerably with the strain (not shown). This sensitivity to the strain level is more pronounced for out-of-plane bonds which increase 8 times their initial strength at 13.8% of strain. Thus, the range of transmitting frequencies and the vibrational band gap increase and, after ∼7.3% of strain, new transmission peaks emerge at the edge of the acoustic branch. Consequently, κ ph decreases for low strain levels and then it increases reaching values close to the initial ones at zero strain independently of the temperature (see right panel in figure 17(b)). In the case of unstrained contact regions, κ ph only decreases with increasing strain level because of the continuous suppression of high frequency in-and out-ofplane modes. Note that the strain dependence of κ ph , considering setup II, for hexagonal boron-nitride and phosphorene is similar to that obtained by employing the standard model for unixial strain (see dashed lines in figure 17(b)). However, for MoS 2 monolayer an increase of the transmission probability at low frequencies was obtained and, hence, κ ph only increases. In fact, by imposing the extra force in the periodic direction (y -direction), transport channels corresponding to low frequencies were blocked, which strongly influences the conduction through out-of-plane modes.

2.10.7.Substrate and curvature effects on phonon transmis-
sion across graphene. Graphene is generally supported by a substrate, which affects its vibrational properties in more than one ways. Van der Waals interaction couples the vibrational modes of the substrate to those of graphene, which is a source for scatterings. Also, graphene adheres well to the substrate and significant curvature is induced due to irregularities of the substrate. Sometimes steps or holes are introduced intentionally to control the electronic properties. These give rise to formation of bends, ripples, wrinkles, and bubbles [122][123][124][125][126], which not only affect the electronic properties but vibrational properties are also altered significantly.
These effects were investigated by using GF techniques and a structureless substrate model. Graphene was modeled within the DFTB approach, whereas the coupling to the substrate was included with a Lennard-Jones potential as [75] where r stands for the distance between a carbon atom to the substrate, and σ determines the equilibrium distance between the layer and the substrate. The atomic positions are optimized using force contributions from the LJ and carbon-carbon interactions. For a flat substrate, the coupling affects the transmission of low frequency portion of the out-of-plane phonons only. This is because V LJ is weak compared to the inter-atomic interactions, and in-plane and out-of-plane vibrational modes are decoupled in flat graphene. There appears a gap for out-ofplane acousitc modes around zero energy. The size of the gap is proportional to LJ 1/2 (see figure 18(a)). Imposing a step at the substrate, graphene is bent. The local curvature depends on the step height h s and the interaction strength LJ . Namely, the radius of curvature decreases with increasing h s and stronger LJ . Smaller radius of curvature couples in-plane and out-of-plane modes stronger and a significant reduction of conductance takes place (figures 18(b) and (c)). Transmission across kinked graphene is also shown in figure 18(c), where the kink is induced by adsorption of lines of hydrogen atoms [12,126,127]. For h s = 10 Å , and LJ = 40 meV, κ is reduced by 44%, 11%, 8%, and 7% at 50 K, 300 K, 500 K, and 1000 K, respectively. The smaller h s , the less is the reduction of conductance. Increasing the interaction strength, increases the reduction. Reductions of 5%, 11%, and 47% are reported for room temperature conductivity for LJ = 20 meV, 40 meV, and 160 meV, respectively. One observes that, for a given LJ conduction reduction is more pronounced at lower temperatures for all h s values ( figure 19(a)). This indicates that low frequency phonons are scattered more effectively at substrate steps. Reductions as large as 45% are possible when LJ is 160 meV ( figure 19(b)). Assuming that resistances due to distant steps are additive, it should be possible to reduce thermal conductance by at least an order of magnitude using stepped substrates.

Generalized quasi-classical Langevin (QCL) approach
Next we briefly discuss the quasi-classical generalized Langevin equation (GLE) which can connect the non-equilibrium Green function and the MD approaches to thermal transport and give a rather intuitive view on the NEGF quantities. The quasi-classical GLE provide an alternative path to interactions in the phonon system which is not limited to anharmonicity/phonon-phonon scattering, but can also be extended to include for example electron-phonon interaction effects such as Joule heating.
Molecular dynamics based on either a force-field model or DFT calculations can provide an 'unbiased' approach to anharmonic couplings in thermal transport in the sense that one does not have to pick the most important anharmonic terms 'by hand'. With standard non-equilibrium molecular dynamics (NEMD) the atomic motion is described using the classical Newtons equations while the two (finite) 'lead' reservoirs are treated as heat-baths in contact with a thermostat [128,129]. This approach has for example been used to examine non-linear thermal effects such as thermal rectification in nano-junctions [130,131]. For similar and many other important problems in nanodevices, one may assume the anharmonicity is restricted to a central 'device' region defined by a general non-linear force, F c , while the leads are harmonic each with a linear coupling to the central region. For a semi-infinite harmonic crystal one may obtain the Green function for the classical equation of motion (EOM) and include the effect on the equation of motion restricted to the central region, much in the same spirit as the one-particle quantum transport with NEGF. The retarded Green function for the harmonic crystal (α ∈ L, R) is the same for classical and quantum harmonic oscillators, In both the classical case one may eliminate the harmonic degrees of freedom of the lead reservoirs, u α (used for both left and right), following Adelman and Doll and others [12,37,132]. This is accomplished by solving their equations of motion using the retarded Green function (matrix), g α , together with the initial conditions, u α (0) and u α (0), and 'driving force' given by the motion of the central region, u C , coupled into the lead by the spring, Φ αC , and propagated, (90) This may be inserted into the EOM for u C . The resulting equation is the GLE for the central coordinates u C (mass-scaled), (91) where we have used the retarded self-energy, (92) and introduced the propagated initial conditions of the reservoir as a driving force, ξ, in the central region, The central point is to treat the initial conditions of the reservoirs as stochastic variables determined from the thermal equilibrium assumption of the reservoirs with temperatures T α . In this way the GLE is a stochastic differential equation. The stochastic forces are for the classical GLE described by zero-mean Gaussian noise with a correlation function based on Boltzmann statistics, Figure 19. Reduction of thermal conductance due to step height (a) and interaction strength (b) are given as the ratio of conductance of deformed graphene to that of flat graphene. The figure is reprinted from [75], with the permission of AIP Publishing.
This also establishes the fluctuation-dissipation theorem relating the deterministic and stochastic forces due to the reservoir coupling, namely the friction and fluctuating forces. In the limit of a fast relaxing reservoir (compared to the time-scale of the motion in the central region) we get Γ(ω) ≈ ∂ ω Γ(0) ω, the usual white noise and time-local (frequency independent) friction. However, for a realistic description of the reservoirs one should use 'colored noise', that is, time-correlated noise representing the correct motion in the reservoirs. This has been used for classical simulations of inter-atomic interactions at surfaces and heat transport [46,133,134]. To obtain the heat flux to a reservoir one calculates the (time averaged) power between the central region and a reservoir as the work done by the force originating from central region, This can also be written as the work done by the forces, F α , from the α reservoir on the central system so one only consider variables in the central region, with, However at low temperatures compared to the phonon energies, which could be room temperature for e.g. carbon nanosystems or molecular junctions, the occupation of modes should be described quantum mechanically using the Bose-Einstein distribution. Due to the identical equation of motion for classical and quantum harmonic variables one may derive the GLE in the quantum regime following the same algebra as for the classical variables and obtain the Landauer/Caroli formula (see equation (50)) for harmonic, steady-state transport [37]. However, the non-commutative nature of the operators result in a non-classical noise, ξ , not compatible with molecular dynamics simulation. On the other hand, following Wang [12,135] one may adopt a quasi-classical [136] approach where the noise is treated by a symmetrized expression which allows for GLE molecular dynamics simulations with a colored noise including quantum fluctuations given by, This approach yields the ballistic Landauer regime for low temperatures and the classical, diffusive regime at high temper atures. However, the treatment of the anharmonic force in the quantum regime is approximative, but it seems to require an artificially strong anharmonicity for it to break down which is probably not often important for realistic systems and potentials [137][138][139].
In practice one can use the velocity Verlet propagation using a noise history generated brute-force from random numbers in frequency space which is Fourier transformed to time domain [135]. A damping may be applied to limit the length of the correlation in time [140,141] Alternatively the noise can Local temperatures in an atomic carbon chain between graphene electrodes (initially at room temperature) with a heat source due to Joule heating by an electronic current calculated using MD simulations with the GLE and the Brenner potential [146]. ((a) and (b)) Calculated using the harmonic approximation, while ((c) and (d)) including anharmonicity. In (b) and (c) we show the local heating as a function of the applied voltage giving rise to the current evaluated in the regions shown in the inset. Reproduced from [147]. CC BY 4.0. be generated 'on the fly' [142], or the noise and friction can be represented using auxilary degrees of freedom [143,144]. The motion at zero temperature reproduce the zero-point motion of the atoms in the harmonic limit.
The quasi-classical GLE can thus handle both low temperature quantum freeze out of phonons as well as the anharmonic interaction approximately beyond perturbation theory. Both of these effects need to be included for carbon systems where the Debye temperature is above 2000 K. It has been employed for graphene [140], nanotubes, ribbons with edge roughness [145], and constrictions [46]. For example for a 5-dimer wide graphene constriction it was found that although the local temperature profile agreed well for classical (equation (94)) and quasi-classical MD (equation (98)) around room temperature, the heat current was overestimated by more than 100% in the classical MD.
One may include the coupling to electrons in the GLE as friction and fluctuating forces relevant for thermoelectrics. In this case the additional deterministic and fluctuating forces are evaluated for electron reservoirs out of equilibrium e.g. in the presence of an electronic current. In the presence of current the fluctuating force contains an additional component which can be associated to Joule heating, while additional deterministic 'electron-wind' forces are due to momentum transfer [148,149]. Using the GLE one may study the heat distribution and heat flow, not due to a temperature difference between two reservoirs, but instead generated by the local Joule heating where the electron-phonon scattering is active. This is illustrated in figure 20 where the GLE is used to calculate the local heating of a carbon atomic chain between graphene in the presence of an strong electrical current. The scattering of electrons will 'pump' certain high-energy phonons in the chain region which is lost by coupling to the phonons in the leads. The inclusion of anharmonicity is especially important as seen by comparing harmonic and anharmonic simulations in figures 20(a)-(d), respectively. It is seen here how the coupling to low-energy phonons enable more efficient cooling.
We refer the reader to a recent review on QCL method where further examples on thermal transport as well as applications to fields other than thermal transport are addressed [16].

Kubo-Greenwood method
The Kubo-Greenwood method is devised for investigating bulk transport properties of disordered materials. It is a real-space implementation of the Kubo formula for phonon propagation, by which phonon dynamics and thermal transport are interconnected. Dynamical information is extracted from the time evolution of an initial wavepacket, where the time evolution operator is expanded in terms of Chebyshev polynomials. A substantial increase in computational efficiency is achieved by employing the Lanczos technique, which avoids matrix inversions. The method has been originally developed for and tested on carbon based materials [59,60,62].
In the Kubo-Greenwood approach only the harmonic part of the Hamiltonian (see equation (1)) is taken into account. Following Allen and Feldman [150], conductivity σ along x-direction is defined as where Ω is system volume, ω stands for the frequency of the applied ac temperature gradient, and J x is the x-component of the flux operator, Here, X i is the equilibrium position of the atom, which the ith degree of freedom belongs, x i is the displacement operator, and ẋ j is the velocity operator. J x can be written in terms of annihilation and creation operators as J x = m,n J x mn a † m a n , where with X being the diagonal matrix of equilibrium positions and |n is the nth egenstate of Φ. It was also shown by Allen and Feldman [150] that the conductivity can be written as The above relation can also be expressed as where √ Φ = n ω n |n n| is used. Since the dynamical matrix is positive definite, it has precisely one positive-definite square root. We also make use of the fact that Φ|n = ω 2 n |n , and the matrices Φ and √ Φ have simultaneous eigenstates.
√ Φ], thermal conductance of a 1D system is obtained as (104) Recalling the definition of conductance in terms of Green functions (see equation (51)), one obtains the transmission function as which is of the same form with the electronic transmission function as derived from the Kubo-Greenwood formula [151,152] E being energy and H el the electronic Hamiltoian.
In the diffusion regime, transmission function is expressed as where D max (ω) is the maximum value of time-dependent diffusion coefficient D(ω, t). It is defined in terms of the meansquare displacement as where χ 2 (ω, t) = (X(t) − X(0)) 2 ω . Computation of D(ω, t) makes it possible to deduce the transport coefficients such as D max (ω), the mean-free-path and the transmission coefficient ζ(ω). Therefore computation of mean-square displacement lies at the heart of Kubo-Greenwood approach. It can be expressed as where U(t) is the time evolution operator. The numerator can be computed through the average of random phase states as N ψ|[X, U(t)] † δ(ω 2 − Φ)[X, U(t)]|ψ , N being the number of degrees of freedom. One should notice that the bra-ket is nothing but a projected DOS associated with the vector [X, U(t)]|ψ . The bra-kets are to be computed using the Lanczos method and the time development is calculated using Chebyshev expansion of U(t). These methods are explained below.

Lanczos tridiagonalization and the continued fraction methods
Computationally efficient algorithms are required in order to be able to simulate large systems. The Lanczos method [153] is a powerful method which is originally developed and widely used in electronic structure calculations and it is equivalently useful in phonon calculations. The first step of the method is tridiagonalization of the dynamical matrix, or the Hamiltonian in the case of electrons. Starting with a seed vector |ψ 1 , a basis set {|ψ n } is built recursively as a n = ψ n |Φ|ψ n (111) |ψ n+1 = Φ|ψ n − a n |ψ n − b n−1 |ψ n−1 , where n > 1 and b 0 = 0. The recursion coefficients a n and b n are the diagonal and the off-diagonal elements of the dynamical matrix in this basis, respectively.
One chooses the seed vector as a random phase state j representing the degree of freedom, and ϕ j is a random number within [0, 2π]. With a few number of such random phase states, the desired accuracy for the DOS is obtained. Using one can show that The continued fraction in equation (118) is terminated after the coefficients a n , b n converge at step N. One names the continued fraction as G 1 and defines G 2 with and similarly For the converged recursion coefficients a N and b N , the below equality should be satisfied (121) Consequently one finds As an example, the recursion coefficients and the DOS of 2D graphene are shown in figure 21.

Wavepacket propagation and Chebyshev expansion
In the Kubo-Greenwood approach, the quantity of central importance is the mean-square displacement, χ 2 (ω, t). Similar to DOS, χ 2 (ω, t) can also be computed using the Lanczos technique outlined above, provided that the time development is evaluated first. For the time evolution operator to be computed efficiently, Chebyshev polynomials of type I are employed. Since the domain of Chebyshev polynomials is [−1, 1], one needs to shift and scale the vibrational spectrum as Φ = (Φ − a)/2b, a being the mid-frequency of the spectrum and spectral range is 4b. Owing to the fact that dynamical matrix and Hamiltonian have the same spectrum, time evolution operator is defined in terms of Φ and can be expanded as with c n (t) = 2 π(1 + δ n,0 ) Hence the commutator takes the form [X, U(t)]|ψ = n c n (t)[X, Q n (Φ)]|ψ . Chebyshev polynomials fulfill the recurrence relation where Q 0 (Φ ) = 1 and Q 1 (Φ ) = Φ [154]. Consequently, the commutator takes the following form, Denoting one writes with |β 0 = 0 and |β 1 = [X, Φ ]|ψ by definition. For obtaining |β n , one needs |α n , which satisfies with |α 0 = |ψ and |α 1 = Φ |ψ . As a result, the computation of time evolution for a finite ∆t is equivalent to obtaining the following, The upper limit N poly is to be chosen depending on the evolution step and the band width. Once U(m∆t)|ψ and [X, U(m∆t)]|ψ are obtained, one proceeds as and That is, the time evolution is computed step by step using the expansion coefficient c n defined in equation (129). The expansion coefficients are computed numerically on a Chebyshev-Gauss grid of N points using the quadrature formula

Applications of the Kubo-Greenwood method
Carbon nanostructures and gaphene-based materials have extremely high mechanical strength and stiffness, which gives usually rise to high thermal conductivities [155][156][157]. In view of thermoelectric applications, lower thermal conductivity is sought, while preserving electronic transport, which remains a very challenging quest. Here we explore several applications of the methods developed in prior sections and show to which extent, it is possible to tune the phonon transport properties in the adiabatic approximation for low-dimensional mat erials. Although this remains a restricted quantitative view of the problem, since elastic phonon scattering processes also contribute to the overall thermal conductance, the obtained scaling properties of the phonon mean free path with the disorder nature and its strength provide guidelines to optimize the material structure for desired thermal/thermoelectrical response. In this part, we discuss the phonon mean free path frequency-dependent profiles in low-dimensional materials such as CNTs [60,158], graphene-based materials [159], and graphene-boron nitride hybrid 2D membranes (including polycrystalline structures) [59,160], obtained with Kubotype and Green function methods. It is known that isotopic or Anderson-type disorder can reduce thermal conductivity in quasi-one-dimensional structures of sp 2 hybridized carbon, and actually tuning the thermal conductivity of disordered CNT with isotope disorder has been proposed to reach the insulating regime of thermal glass [161]. One key characteristic of thermal transport is the phonon mean free path, of which frequency dependence brings essential information concerning the impact of disorder on vibrational properties. Similarly to carbon nanotubes, GNRs, which look like unfolded CNTs, are possible to synthesize either from top-down or bottom-up approaches, but ribbon edge irregularities should suppress the thermal conductance for small width GNRs [63]. Finally 2D heterostructures made from mixing graphene with disordered (isotope impurities) GNRs [162], or graphene-boron nitride membranes and polycrystalline mixture [163] offer more freedom for trying monitoring vibrational properties. We note that the force constant matrices were built using the 4NNFC method in all examples of the KG scheme in this review.
We first consider a CNT(7,0) with 10.7% of 14 C impurities, investigated by Savic and coworkers with the GF method [161]. Figure 22 (left panel) illustrates the disordered nanotube while figure 22 (right panel) shows the phonon wavepacket mean free path (MFP) versus mode frequency. Note that MFP is extracted from the saturation of the diffusion coefficient, D(ω, t), that is computed using the KG scheme (not shown here). Indeed such saturation of D(ω, t) to a maximum value characterizes diffusive transport. The results compare very well with those obtained with the GF approach, as well as with the analytical result (ω) = (12aN uc N ch (ω))/(π 2 f | ∆M M | 2 ρ 2 uc (ω)ω 2 ) where a is the length of the lattice vector in the translational direction, N uc is the number of atoms in each unit cell, ρ uc is the density of states per unit cell, and f is the percentage of isotopic impurities having mass difference ∆M [60].
We note that at the singularities of the phonon spectrum, the disorder induces a broadening of states which impact on the numerical MFP, and it remains slightly larger than the  one obtained with the analytical expression, but such difference remains minnor compared to the overall agreement of the results. As a general observation, one also remarks that (ω) generally increases for lowering frequency and for a frequency range between [250-2500] cm −1 , its variation is about three orders of magnitudes. The low-energy acoustic modes display much longer MFP, which is out of reach from our simulations and finite size effects. Optic modes with frequency above 1500 cm −1 have MFP in the order of few nanometers, corresponding to elastic scattering times of few femtoseconds. It is clear that optic modes are strongly affected by isotope impurity scattering but also that the thermal conductance will be dominated by low-energy modes, weakly sensitive to isotopic disorder.
Interestingly, a power law behavior of phonon lifetimes due to 13 C isotopic doping in graphene was found theoretically, with 1/τ ph-imp ∼ ω 3 (resp. ω 2 ) for in-plane (out-ofplane) acoustic modes, whereas phonon lifetimes of optical phonon modes were also found to be considerably smaller [169] (close to the Γ-point). At the Γ-point for the in-plane optical modes, the estimated value of the phonon lifetime due to 13 C impurities was τ ph-imp 0.83 ps for an isotopic atomic density of 0.5, while τ ph-imp 16.7 ps for naturally occurring carbon [164]. Besides the electron-phonon impact on phonon lifetime is τ ph-imp 0.6 ps so for high concentrations of isotope impurities, the scattering rates for the isotopic impurity is as large as that for the electron-phonon interaction.
We next contrast the case of disordered nanotubes with zigzag (ZGNR) and armchair (AGNR) graphene nanoribbons of different widths and with edge disorder. Figure 23 shows the comparison of (ω) in AGNR and ZGNR of approximately equal widths ( 17 nm) and both with similar amount of edge disorder. The ribbon widths are defined with the number of zigzag chains N z , whereas the edge disorder is dictated by the density of edge defects (removed carbon atoms at the edges) and is chosen to be 10% of edge C-atoms (figure 23). One observes large differences in the MFPs between ZGNR and AGNR. Indeed if MFP remains well below 200 nm for frequencies below 800 cm −1 , substantial fluctuations of MFP are seen for ZGNR, with peaked values above 500 nm at much lower frequencies. Additionally, differently to the case of ZGNR, AGNR exhibit MFPs as large as several micrometers around 850-1100 cm −1 when MFP for ZGNR is minimum for such frequency range. Overall the variations of (ω) are strikingly opposite for AGNR and ZGNR, but the interpretation of such feature remains elusive.
At a fixed disorder strength (ω) decays with decreasing the ribbon width (not shown here). This behavior can be rationalized with the fact that the scattering rate decreases with increasing width, a behavior generic for electron transport in both disordered CNTs and GNRs. One notes that, for low frequency modes, the MFP are several hundreds of nanometers, and due to large values of the transport channels the possibility to observe Anderson localization and exponentially damped thermal conductance is out of reach, as previously discussed [161]. The thermal conductances for AGNR and ZGNR are shown in figure 24. The difference between the pristine thermal conductances of AGNRs (dashed lines) and ZGNRs (solid lines) stem from the anisotropy observed in the phonon dispersion. Actually, the phonon MFPs of AGNRs are smaller than those of ZGNRs for the frequencies dominating the thermal conductance at low temperatures (see figure 23). These two factors explain why the AGNR thermal conductance is smaller than that of ZGNRs for a given edge disorder and ribbon geometry. Although the difference is found to be reduced as the ribbon width increases (not shown), coherent phonon propagation is clearly sensitive to the ribbon edge shape.
Large width ZGNR are explored in figure 25. We first report evolution of the wave packet dynamics for different frequencies for ZGNR with 10% edge disorder, through the timedependent diffusion coefficient (D(ω, t)) at selected phonon modes (inset). At short times, the linear increase in D(ω, t) indicates ballistic transport, which is followed by a saturation and further a decay of D(ω, t) owing to scattering and interference effects. The low decay observed for D(ω, t) however suggest very long localization lengths. The saturation of D(ω, t) to a maximum value dictates the value of the elastic MFP. In figure 25 (main frame), one also sees a decay of (ω) with decreasing the ribbon width for a fixed edge disorder. This behavior is explained by the reduction of the scattering rate with increasing width, an effect previously derived for electron transport in both disordered CNTs and GNRs [159].
We note that the conductivity of edge disordered GNRs ican be roughly estimated using σ = κ L/A, where A denotes the cross section area of the ribbons (taken as the interplane distance of graphite layers). At room tempetaure, σ = 1006 W m −1 K −1 and σ = 671 W m −1 K −1 with edge disorder 10% and 15%, respectively for L = 2 µm.
We next investigate the phonon propagation and thermal conductivity in hybrid boron nitride (BN) and graphene sheets (with a particular geometry depicted in figure 26(inset)). By inserting BN islands inside a graphene membrane, one expects to tune the heat transfer properties with less damage to electron transport. We here show the possible variation that can be obtained for such hypothetical morphologies, whereas the case of hybrid polycrystalline hBN-graphene (more realistic for large-scale 2D material growth) will be analyzed afterwards.
The calculated transport MFPs are plotted in figure 26 (bottom panels). The DOS is very low around 537 cm −1 for graphene, corresponding to the intersection of the dispersions of the ZA and ZO modes at the K point of the Brillouin zone. For pristine BN, these modes repel each other owing to the underlying broken-sublattice symmetry, which results in a lower DOS compared to pristine graphene [165]. Accordingly, the out-of-plane modes are significantly scattered at these frequencies with increasing the BN density. The additional source of disorder under consideration is an isotopic distribution on a graphene host with a 50% BN concentration with domain size d = 2 nm. Boron has two stable isotopes, 11 B and 10 B with natural abundances of 80.1% and 19.9%, respectively. Figure 26 (bottom-panel right) show that, in the low frequency regime, the MFPs for the pure 10 B case are longer than those for the pure 11 B case. Though the mass difference between 10 B and 12 C is larger than the difference between 11 B and 12 C, the negative mass difference is suppressed by the positive mass difference between 14 N and 12 C, leading to larger scattering by 11 B than by 10 B. Hence, low-frequency modes are strongly affected to the average mass of the BN clusters while at high frequencies, the deviation of atomic masses from the average are more effective. One also finds that the curves for the mixture case do not simply fall in between the curves of the isotopically pure cases. At low frequencies, they follow the 11 B curve, which is more abundant, but at higher frequencies, this tendency is lost.
Large-scale coplanar G-hBN heterostructures were successfully fabricated using chemical vapor deposition (CVD), with varying domains of sizes ranging from tens of nanometers to millimeters [166,167]. Fast CVD growth results in polycrystalline materials with varying grain sizes and structure morphologies, which demand to characterize how electronic and thermal properties are affected by grain boundaries, and if any scaling transport behavior actually takes place in such structures [156,168,169]. Non-equilibrium MD simulations were carried out to compute the thermal conductivity for 16 grain sizes between 1-1000 nm while changing the concentration of hBN [160]. For small average grain size, the minimum value of thermal conductivity occurs near 70% hBN agreeing with prior KG results [59]. This is explained by the fact that the thermal conductance for the G-hBN interface is lower than that of the hBN-hBN and G-G interfaces. For larger grain sizes, where the GBs no longer dominate the thermal transport, a monotonic scaling of κ with hBN grain density is found, consistent with the fact that the thermal conductivity of pristine hBN is lower than that of pristine graphene.

Summary and outlook
In summary, computational methods on quantum thermal transport have proven to be suitable for investigating all transport regimes, from ballistic to diffusion and localization. The reviewed methods have their own advantages, depending on the size and dimensionality of the system. For molecular junctions, the GF approach is extremely useful because the dominant scatterings are elastic, which are fully accounted for within the GF methodology. For quasi-one dimensional systems, the GF approaches can be implemented very efficiently using the recursion decimation algorithms outlined in section 2.7. These algorithms are also useful for 2D systems, when periodic boundary conditions are applicable. With the use of recursion algorithms, low dimensional systems with more than 10 6 atoms can be treated numerically and exactly with GF methods. Elastic mean-free-paths can be obtained quite accurately by employing GF and KG methods. The KG scheme is capable of handling even larger systems, and without any dimensional restrictions. It has been applied successfully to transport problems in disordered systems. The sources of disorder can be grain boundaries, point defects or even the glass-like structure. The main motivation to use this method is the fact that this KG method scales linearly, whereas most of the methods including GF scale with N 3 , N being the number of degrees of freedom in the system.
On the other hand, the reviewed methods cannot address inelastic processes on the same footing with the elastic scattering mechanisms. But since inelastic mean-free-paths are generally at the order of micrometers, GF and KG methods should be good approximations at the nano-scale. QCL approach, on the other hand, offers a powerful way to incorporate inelastic processes for relatively small systems. For extended systems, multi-phonon processes should be accounted for either by extending the reviewed methodology or externally, e.g. by applying Mathissen's rule.