This is the accepted version of the article: Sevinçli H., Roche S., Cuniberti G., Brandbyge M., Gutierrez R., Sandonas L.M.. Green function, quasi-classical Langevin and Kubo-Greenwood methods in quantum thermal transport. Journal of Physics Condensed Matter, (2019). 31. 273003: - . 10.1088/1361-648X/ab119a. Available at: https://dx.doi.org/10.1088/1361-648X/ab119a Green Function, Quasi-Classical Langevin and Kubo-Greenwood Methods in Quantum Thermal Transport H Sevinçli†‡ , S Roche§k , G Cuniberti¶+ , M Brandbyge∗ , R Gutierrez¶ , L Medrano Sandonas¶\ † Department of Materials Science and Engineering, Izmir Institute of Technology, 35430, Urla, Izmir, Turkey. § Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain k ICREA, Instituci o Catalana de Recerca i Estudis Avançats, 08070 Barcelona, Spain ¶ Institute for Materials Science and Max Bergmann Center of Biomaterials, TU Dresden, 01062 Dresden, Germany + Dresden Center for Computational Materials Science and Center for Advancing Electronics Dresden, TU Dresden, 01062 Dresden, Germany ∗ Department of Micro- and Nanotechnology (DTU Nanotech), Center for Nanostructured Graphene (CNG), Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark \ Max Planck Institute for the Physics of Complex Systems, 01187 Dresden, Germany E-mail: ‡ haldunsevincli@iyte.edu.tr Abstract. 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 materials design. Shrinkage of length scales has also changed the experimental and theoretical methods to study thermal transport. Not surprisingly, novel approaches has 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, also 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. GF, QCL and KG methods in quantum thermal transport 2 Contents 0.1 List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1 Introduction 4 2 Green function method 5 2.1 The Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Six phonon Green functions . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2.1 Contour ordering . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Partioning scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4 Expression for the heat current and Landauer-type formula . . . . . . 10 2.5 Generalization to multiple reservoirs. . . . . . . . . . . . . . . . . . . . 11 2.6 Self-consistent reservoir approach . . . . . . . . . . . . . . . . . . . . . 12 2.7 Decimation algorithms for efficient computation of Green functions . . 13 2.7.1 Recursive algorithm. . . . . . . . . . . . . . . . . . . . . . . . . 13 2.7.2 Renormalization-decimation algorithm. . . . . . . . . . . . . . 14 2.8 Surface Green functions (SGF). . . . . . . . . . . . . . . . . . . . . . . 14 2.8.1 Generic reservoir as a simple model for SGFs. . . . . . . . . . . 15 2.9 A minimal force constant model for graphene based materials . . . . . 16 2.10 Applications of the Green Function Method . . . . . . . . . . . . . . . 18 2.10.1 Defected carbon nanotubes. . . . . . . . . . . . . . . . . . . . . 18 2.10.2 Graphene with molecular antiresonances. . . . . . . . . . . . . 21 2.10.3 A hybrid nano-structuring scheme using geometrical structuring and isotopical disorder. . . . . . . . . . . . . . . . . . . . . . . 22 2.10.4 Functionalized graphene grain boundaries . . . . . . . . . . . . 24 2.10.5 Thermal anisotropy in planar puckered structures. . . . . . . . 25 2.10.6 Strain engineered two-dimensional materials. . . . . . . . . . . 27 2.10.7 Substrate and curvature effects on phonon transmission across graphene. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3 Generalized Quasi-Classical Langevin (QCL) approach 30 4 Kubo-Greenwood method 33 4.1 Lanczos tridiagonalization and the continued fraction methods . . . . 35 4.2 Wavepacket propagation and Chebyshev expansion . . . . . . . . . . . 36 4.3 Applications of the Kubo-Greenwood Method . . . . . . . . . . . . . . 38 5 Summary and Outlook 43 6 References 44 GF, QCL and KG methods in quantum thermal transport 3 0.1. List of Symbols Table 1. a Here is the list of symbols. Symbol Explanation xiα Displacement of ith atom along α-direction ẋiα Velocity of ith atom along α-direction mi Mass of ith atom. M Diagonal matrix of atomic masses. Vij,αβ = ∂ 2 E/∂xiα ∂xjβ 0  Vij,αβ Harmonic force constants. 3  Vijk,αβγ Third order force constants. Vijk,αβγ = ∂ E/∂xiα ∂xjβ ∂xkγ 0 4  Vijkl Fourth order force constants. Vijkl,αβγθ = ∂ E/∂xiα ∂xjβ ∂xkγ ∂xlθ 0 H µν (µ, ν) part of the Hamiltonian. µ, ν ∈ {L, C, R} √ uiα Mass normalized displacement. uiα = mi xiα √ u̇iα Mass normalized velocity. u̇iα = mi ẋiα Φij,αβ = ∂ 2 E/∂uiα ∂ujβ 0  Φij,αβ Mass normalized harmonic force constants. 3  Ψijk,αβγ Mass normalized third order force constants. Ψijk,αβγ = ∂ E/∂uiα ∂ujβ ∂ukγ 0 = ∂ 4 E/∂uiα ∂ujβ ∂ukγ ∂ulθ 0  Ψijkl,αβγθ Mass normalized fourth order force constants. Ψ Pijklαβγθ Bij,αβ (ij, αβ) element of diagonalizing matrix. ∗ jβ B Pij,αβ Bjk,βγ = δik δαγ uq Normal coordinate. uq = iα Bq,iα uiα Gr Retarded Green’s function. Ga Advanced Green’s function. Gt Time-ordered Green’s function. G< Lesser Green’s function. G> Greater Green’s function. 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. ΓL/R = i(ΣL/R − ΣL/R ) ζ(ω) Transmission spectrum. fB (ω, T ) Bose-Einstein distribution function. κ(T ) Conductance. J Thermal current. σ Thermal conductivity. `elastic Elastic mean-free-path ξ Driving force ξ µ (t) = Cµ µ µ R Φ g (t)u̇ (0) Sµ Gaussian noise S µ = dt eiωt hξ µ (t)ξ µ (0)i a Internal use, not for publication... GF, QCL and KG methods in quantum thermal transport 4 1. Introduction Thermal transport plays a central role in most of the state-of-the-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–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–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 real-space 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-of-plane 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 GF, QCL and KG methods in quantum thermal transport 5 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–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 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. 2. 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–24] 2.1. The Hamiltonian The total energy of the system Etot can be expanded in terms of the displacements of the constituent atoms from their equilibrium positions. Denoting mi and xiα as the mass of ith atom and its displacement along α-direction (α = x, y, z), one defines √ mass-normalized displacements as uiα = mi xiα and writes the Hamiltonian as 1X 2 1 X H= u̇iα + Φij,αβ uiα ujβ + · · · 2 i,α 2 ij,αβ 1 X 1 X + Ψijk,αβγ uiα ujβ ukγ + Ψijkl,αβγθ uiα ujβ ukγ ulθ , (1) 3! 4! ijk ijkl αβγ αβγθ where Φij,αβ = ∂ 2 Etot /∂uiα ∂ujβ are the elements of the dynamical matrix. Ψijk,αβγ =∂ 3 Etot /∂uiα ∂ujβ ∂ukγ and Ψijkl,αβγθ =∂ 4 Etot /∂uiα ∂ujβ ∂ukγ ∂ulθ 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, Hint , whereas the harmonic part (Hharm ) defines an eigenvalue problem. GF, QCL and KG methods in quantum thermal transport 6 One can define |ui as the normalized column vector of displacements uiα , and hu| as its transpose such that hu|ui = 1, and denote the eigenvalue problem defined by the dynamical matrix, Φ, as Φ|uq i = ωq2 |uq i (2) with uq being the eigenvector and ωq2 is the corresponding eigenvalue with index q. The equivalent form without mass-normalization is a generalized eigenvalue (2) problem, V (2) |xq i = ω 2 M |xq i, where Vij,αβ = ∂ 2 Etot /∂xiα ∂xjβ are the harmonic force constants and M is the diagonal matrix of atomic masses. The annihilation and creation operators are defined as as r   ωq i aq = uq + u̇q , (3) 2~ ωq r   + ωq i aq = uq − u̇q , (4) 2~ ωq which satisfy the commutation relation [aq , a+ q 0 ] = δqq 0 , and diagonalize the harmonic 1 Hamiltonian, Hharm = q ~ωq (a+ P q a q + 2 ). Also note that [uiα , u̇jβ ] = i~δij δαβ . 2.2. Six phonon Green functions In equilibrium, it is enough to to define one Green function, which is generally chosen as the retarded GF. In nonequilirium cases, it is common practice to define six Green functions, which are not linearly independent. Their definitions are, [17] i G> (u1 , u2 ) = − hu1 (t1 ), u2 (t2 )i, (5) ~ i G< (u1 , u2 ) = − hu2 (t2 ), u1 (t1 )i, (6) ~ Gt (u1 , u2 ) = θ(t1 − t2 )G> (u1 , u2 ) + θ(t2 − t1 )G< (u1 , u2 )i, (7) Gt̃ (u1 , u2 ) = θ(t2 − t1 )G> (u1 , u2 ) + θ(t1 − t2 )G< (u1 , u2 )i, (8) Gr (u1 , u2 ) = Gt (u1 , u2 ) − G< (u1 , u2 ), (9) Ga (u1 , u2 ) = Gt (u1 , u2 ) − G> (u1 , u2 ), (10) 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] Gr = θ(t1 − t2 )(G> − G< ), (11) Ga = −θ(t2 − t1 )(G> − G< ), (12) Gt − Gt̃ = (θ(t1 − t2 ) − θ(t2 − t1 ))(G> − G< ), (13) Gt + Gt̃ = G> + G< , (14) Gr − Ga = G> − G< , (15) Gr + Ga = Gt − Gt̃ , (16) GF, QCL and KG methods in quantum thermal transport 7 Figure 1. Contour integration path. 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 Z C(t1 , t2 ) = dτ A(t1 , τ ) B(τ, t2 ). (17) C < In order to calculate C , one chooses a contour C such that t1 is before t2 (see Figure 2a). Then the contour is deformed as C → C1 + C2 (see Figure 2b), so that t1 is on C1 and t2 is on C2 . Using these, one obtains the expression for C < as Z Z C < (t1 , t2 ) = dτ A(t1 , τ ) B < (τ, t2 ) + dτ A< (t1 , τ ) B(τ, t2 ). (18) C1 C2 < Here, one has B in the first term because all τ on C1 are before t2 , which is on C2 . Similarly, one has A< in the second term because t1 is always before τ , which is on C2 . Now, we do the the first integral.On C1 , τ can be before or after t1 , which is expressed as Z dτ A(t1 , τ ) B < (τ, t2 ) = C1 Zt1 −∞ Z > < dt A (t1 , t) B (t, t2 ) + dt A< (t1 , t) B < (t, t2 ). (19) −∞ t1 Likewise, the C2 integral takes the form Z dτ A< (t1 , τ ) B(τ, t2 ) = C2 (a) (b) C C C Figure 2. The contour C, on which t1 is before t2 , is deformed as C → C1 + C2 . GF, QCL and KG methods in quantum thermal transport 8 Zt2 −∞ Z < > dt A (t1 , t) B (t, t2 ) + dt A< (t1 , t) B < (t, t2 ). (20) −∞ t2 Using the identities Ar (t1 , t) = θ(t1 − t)(A> − A< ) and B a (t, t2 ) = −θ(t2 , t)(B > − B < ) (see Equations 11-12), one obtains Z∞ < dt Ar (t1 , t)B < (t, t2 ) + A< (t1 , t)B a (t, t2 ) .  C (t1 , t2 ) = (21) −∞ It is possible to use a matrix representation for the GFs and the self energies as [25]  t −G<  t −Σ<   G Σ G̃ = , Σ̃ = . (22) G> −Gt̃ Σ> −Σt̃ Using the matrix notation, à = B̃C̃, one can easily write −A< = −B t C < + B < C t̃ . making use of the definitions B t = B r + B < and C t̃ = C < − C a (see Equations 9- 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 g α,< (ω, Tα ) = fBα (ω, Tα )(g α,r (ω) − g α,a (ω)), (23) g α,> (ω, Tα ) = (1 + fBα (ω, Tα ))(g α,r (ω) − g α,a (ω)). (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. In this review, we limit our discussions with the the linear part of the Hamiltonian. Applying perturbation expansion, Dyson’s equation with contour ordering can be obtained as Z Z G(τ, τ 0 ) = g(τ, τ 0 ) + dτ1 dτ2 g(τ, τ1 )Σ(τ1 , τ2 )G(τ2 , τ 0 ), (25) and using the Langreth rules one obtains G≶ (ω) = Gr (ω)Σ≶ (ω)Ga (ω) (26) ≶ with Σ originating from couplings between subsystems. 2.3. 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 TL > TR , 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 GF, QCL and KG methods in quantum thermal transport 9 LEFT CENTER RIGHT Figure 3. Partitioning scheme. The infinite system is partitioned into three as the semi-infinite reservoirs (left and right regions) and the central scattering region. Note that the semi-infinite regions consist of identical layers, which makes them scattering free, perfect reservoirs. 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  LL H LC  H 0 H =  H CL H CC H CR  , (27) 0 H RC H RR where 1 X 2 1 X H µµ = u̇iα + Φij,αβ uiα ujβ , (28) 2 i∈µ,α 2 ij∈µ,αβ 1 X H µν = Φij,αβ uiα ujβ . (29) 2 i∈µ,j∈ν αβ The retarded Green function, also referred to as the resolvent, can be written as Gr (ω) = (ω + iδ)2 1 − Φ  −1 , (30) where δ is an infinitesimal positive number and 1 is the unit matrix. The elements of the resolvent read Grqq0 = δqq0 /(ω 2 − ωq2 + iδ) in the diagonal basis, δqq0 representing Kronecker-δ function. Using the partitioning scheme, the dynamical matrix and the resolvent can be written as 3 × 3 matrices, and Equation 30 reads (ω + iδ)2 1 − ΦLL −ΦLC   0  −ΦCL (ω + iδ)2 1 − ΦCC −ΦCR × 0 −Φ RC (ω + iδ)2 1 − ΦRR  LL,r GLC,r GLR,r  G  GCL,r GCC,r GCR,r  = 1. (31) GRL,r GRC,r GRR,r GF, QCL and KG methods in quantum thermal transport 10 The second column of the product matrix consists of the following equalities. (g L,r )−1 GLC,r − ΦLC GCC,r = 0, (32) −ΦCL GLC,r + [(ω + iδ)2 1 − ΦCC,r ]GCC,r − ΦCR GRC,r = 1, (33) R,r −1 RR,r RC RC,r (g ) G −Φ G = 0, (34) where g L/R,r = [(ω+iδ)2 1−ΦLL/RR ]−1 are the retarded Green functions for uncoupled reservoirs. Using Equations 32, 33, 34, Green function for the central region is obtained as −1 GCC,r = (ω + iδ)2 1 − ΦCC − ΣL,r − ΣR,r  , (35) 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 h i ΓL/R = Im ΣL/R,a − ΣL/R,r . (36) 2.4. 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, d J L→C (t) = − hH LL (t)i, dt i  = h H LL (t), H LC (t) + H CL (t) i,  (37) ~ where the Heisenberg equation of motion is used and it is assumed that the degrees of freedom of the left and right reservoirs are not coupled to the degrees of freedom of the central region that are involved in the non-linear part of the Hamiltonian, namely [H LL/RR , Hint ] = 0. This can always be assured with including a high enough number of principal layers of the reservoirs inside the central region. (see Figure 3) Using the commutation relation [ul1 , u̇l2 ] = i~δl1 ,l2 , one finds X J L→C (t) = ΦLC L C l1 c1 hu̇l1 (t)uc1 (t)i, (38) l1 c1 X d CL,< = i~ ΦLC l1 c1 lim G (t, t0 ). (39) 0 t →t dt0 c1 l1 l1 c1 Here, the indices l1 , c1 etc include both the atom index and the Cartesian components (cf Equation 6). Performing contour ordering within the partition scheme, one has GCL,< = GCC,r ΦCL g LL,< + GCC,< ΦCL g LL,a , (40) which leads to the expression Z d X J L→C (t) = i~ lim dt00 GCC,r 00 L,< 00 0 c1 c2 (t, t )Σc2 c1 (t , t ) + · · · (41) t0 →t dt0 c1 c2 00 00 0 GCC,< L,a  c1 c2 (t, t )Σc2 c1 (t , t ) (42) Applying Fourier transformation and symmetrizing the real valued function as J = (J + J ∗ )/2, one finds Z∞ L→C dω ~ω Tr ΣL,> (ω)GCC,< (ω) − ΣL,< (ω)GCC,> (ω) . (43)   J = 2π 0 GF, QCL and KG methods in quantum thermal transport 11 Using Equation 23 for the uncoupled GFs, one can write ΣL,< = ΦCL g L,< ΦLC , = fBL ΦCL (g L,r − g L,a )ΦLC , = fBL (ΣL,r − ΣL,a ), = ifBL ΓL . (44) Similarly, using Equation 24, one gets ΣL,> = i(1 + fBL )ΓL . (45) CC,< CC,> Employing Dyson’s Equation (cf Eqn. 26) for G and G as CC,< CC,r fBL ΓL + fBR ΓR GCC,a ,  G = iG (46) CC,> CC,r fBL )ΓL + (1 + fBR )ΓR GCC,a ,  G = iG (1 + (47) one obtains ΣL,> GCC,< − ΣL,< GCC,> = (fBL − fBR ) ΓL GCC,r ΓR GCC,a . (48) This equality reproduces the Landauer-type formula given Rego and Kirczenow [2]. Z ∞ dω J= ~ω [fB (ω, TL ) − fB (ω, TR )] ζ(ω), (49) 0 2π where TL > TR are reservoir temperatures and ζ(ω) is the transmission spectrum. The expression for the transmission function was obtained by Yamamoto and Watanabe [24] as † ζ(ω) = Tr[ΓL (ω)GCC,r (ω)ΓR (ω)GCC,r (ω)]. (50) 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 TL −TR  T = (TL +TR )/2, thermal conductance is κ = J/(TL −TR ), hence at given temperature k2 T ∞ Z κ(T ) = B dx p(x) ζ(x), (51) h 0 with x=~ω/kB T , and p(x)=−x2 ∂fB /∂x=x2 ex (ex −1)−2 . The function R ∞ p(x) = x2 ex /(ex − 1)2 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 κ=M κ0 , where κ0 = π 2 kB 2 T /3h (52) is the temperature dependent quantum of thermal conductance and M is the number of acoustic modes. [1–3] 2.5. Generalization to multiple reservoirs. The Landauer-type expression can be generalized to multi-terminal 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 † h i ζ µν (ω) = Tr Γµ GCC,r Γν GCC,r , (53) GF, QCL and KG methods in quantum thermal transport 12 where Γµ is the mode broadening due to coupling to reservoir µ, Γµ = Im [Σµ,a − Σµ,r ] (cf. Equation 36). GCC P collects self energy contributions from all reservoirs, GCC,r =[(ω + iδ)2 1 − ΦCC − µ Σµ,r ]−1 . Thermal current from µ to ν is then Z dω J µν = ~ω [fBµ (ω, Tµ ) − fBν (ω, Tν )] ζ µν (ω). (54) 2π The multi-terminal systems are of particular interest for applications such as thermal rectification. [26–31] 2.6. Self-consistent reservoir approach A method which is based on the multiple reservoirs is the self-consistent reservoir (SCR) approach. The method is originally developed for a classical harmonic chain consisting of N atoms. [32–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 T1 and TN , which correspond to the temperatures of the left and right reservoirs, T1 =TL and TN =TR . The rest of the temperatures from T2 to TN −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–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. Figure 4. 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 flow from them. GF, QCL and KG methods in quantum thermal transport 13 2.7. 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   Φ1,1 Φ1,2 0 Φ =  Φ1,2 Φ2,2 Φ2,3  . (55) 0 Φ3,2 Φ3,3 The retarded GF is defined as G(ω) = ((ω + iδ)2 1 − Φ)−1 (56)   G1,1 G1,2 G1,3 =  G2,1 G2,2 G2,3  . (57) G3,1 G3,2 G3,3 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 G2,2 (ω) = ((ω + iδ)2 1 − Φ2,2 )−1 (0) (58) from which the self energy terms arise as     Σ1,1 Σ1,3 Φ1,2 (0)  Σ= = G2,2 Φ2,1 Φ2,3 , (59) Σ3,1 Σ3,3 Φ3,2 and the 2×2 effective dynamical matrix is obtained as   eff Φ1,1 + Σ1,1 Σ1,3 Φ = , (60) Σ3,1 Φ3,3 + Σ3,3 where Φeff is a function of ω. 2.7.1. Recursive algorithm. It is straightforward to generalize the above scheme for an extended tridiagonal matrix,   Φ1,1 Φ1,2 0 ··· 0 0  Φ2,1 Φ2,2 Φ2,3 · · · 0 0     0 Φ3,2 Φ3,3 · · · 0 0  Φ= . (61)   . . .. . .  .. .. .. .. ..   .    0 0 0 · · · ΦN −1,N −1 ΦN −1,N  0 0 0 ··· ΦN,N −1 ΦN,N 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 (n) G(n) = ((ω + iδ)2 − Φn+2,n+2 )−1 , (62) (n+1) (n) (n) (n) Φ1,1 = Φ1,1 + Φ1,n+2 G(n) Φn+2,1 , (63) (n+1) (n) (0) Φ1,n+2 = Φ1,n+2 G(n) Φn+1,n+2 , (64) (n+1) (0) (n) Φn+2,1 = Φn+2,n+1 G(n) Φn+2,1 , (65) (n+1) (0) (0) (0) Φn+2,n+2 = Φn+2,n+2 + Φn+2,n+1 G(n) Φn+1,n+2 , (66) with n=0, . . . , N − 3. Repeating the decimation N −2 times, one arrives at the 2×2 effective GF for the central region consisting of G1,1 , G1,N , GN,1 and GN,N . GF, QCL and KG methods in quantum thermal transport 14 2.7.2. 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   .. .. ..  . . .   0 Φ0,−1 Φ0,0 Φ0,1 0  Φ= , (67)   0 Φ1,0 Φ1,1 Φ1,2 0   .. .. .. . . . where Φi,i =Φ0,0 , Φi,i+1 =Φ0,1 , and Φi+1,i =Φ0,−1 . 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   .. .. ..  . . .   0 Φ(1) (1) (1) 0,−2 Φ0,0 Φ0,2 0  (1) Φ = , (68)   (1) (1) (1)   0 Φ 2,0 Φ 2,2 Φ 2,4 0   .. .. .. . . . with (1) Φ0,0 = Φ0,0 + Φ0,−1 G(0) Φ−1,0 + Φ0,1 G(0) Φ1,0 , (69) (1) Φ0,2 = Φ0,1 G(0) Φ1,2 , (70) (1) (0) Φ2,0 = Φ2,1 G Φ1,0 . (71) (n) Repeating this procedure by using G(n) = ((ω + iδ)2 − Φ0,0 )−1 , the submatrices are computed as (n+1) (n) (n) (n) (n) (n) Φ0,0 = Φ0,0 + Φ2n ,0 G(n) Φ0,2n + Φ0,2n G(n) Φ2n ,0 , (72) (n+1) (n) (n) Φ0,2n+1 = Φ0,2n G(n) Φ0,2n , (73) (n+1) (n) (n) Φ0,−2n+1 = Φ2n ,0 G(n) Φ2n ,0 , (74) which converge very fast compared to the standard decimation method and incorporation of a thousands of layers is possible within a few steps. 2.8. Surface Green functions (SGF). For semi-infinite reservoirs, dynamical matrices have the following form. ..   .. .. .. ..  . . . . .  LL LL LL ΦLL =  . . . 0 Φ10 Φ00 Φ01 0  , (75)   ... 0 0 ΦLL10 ΦLL00 ΦLL 01  LL LL ... 0 0 0 Φ10 Φ00  RR RR  Φ00 Φ01 0 0 0 ...  ΦRR RR ΦRR  10 Φ00 01 0 0 ...  ΦRR =  0 RR RR RR , (76)  Φ10 Φ00 Φ01 0 ...  ..   .. .. .. .. . . . . . GF, QCL and KG methods in quantum thermal transport 15 and the surface terms are obtained as LL,(n+1) LL,(n) LL,(n) LL,(n) Φ0,0 = Φ0,0 + Φ2n ,0 GLL,(n) Φ0,2n , (77) RR,(n+1) RR,(n) RR,(n) RR,(n) Φ0,0 = Φ0,0 + Φ0,2n GRR,(n) Φ2n ,0 . (78) The convergence of the iteration is checked by comparing the difference between dynamical matrices at two consecutive steps. 2.8.1. 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 DOS is approximately constant around the Fermi energy. (see eg. Ref. [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 s=mω02 . The dynamical matrix is ..   .. .. ..  . . . .   ..   . 2ω 2 −ω 2 0  Φ= 0 0 . (79)  . .   . −ω 2 2 2 0 2ω0 −ω0  ... 0 −ω02 2ω02 It was shown by Müller et al. [51] that the elements of the GF are given as 1 eiθ(m+n) − eiθ|m−n| Gmn (ω) = − , (80) ω02 2i sin θ where θ = arccos(1 − ω 2 /2ω02 ), and {m, n} ∈ [1, ∞ ) with m=n=1 corresponding to the on-site term for the end atom. Therefore the SGF is obtained as r ! 1 η2 η2 G11 (ω) = − 2 1 − + iη 1 − , (81) ω0 2 4 GF, QCL and KG methods in quantum thermal transport 16 Figure 5. Real and imaginary parts of the surface Green function as obtained from the generic reservoir model. 1 Re[G 11 ] Im[G 11 ] 0.5 G 11 ( ) 0 2 0 -0.5 -1 0 0.5 1 1.5 2 2.5 3 / 0 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 ω02 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 elementsp of Φ are the same. If one end was free, i.e. if Φ11 = ω02 , −1 then Φ11 = ω0 (ηp/2 − iη 1 − η /4) and G11 = (ω − Φeff eff 2 2 2 2 11 ) and the surface density 2 of states ρ11 = 4 − η /πω0 . That is, the density of states is finite at ω = 0 and decreases monotonically with increasing frequency. 2.9. 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. [53–65] Also, a reoptimized version of the force constant parametrization is available. [66] The method relies on parametrization of interatomic force constants up to fourth nearest neighbors (4NNFC approach). For each neighboring distance (n = 1, . . . , 4) (n) (n) (n) three force constant parameters fr , fti and fto 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  (n)  fr 0 0 (n) F (A,B1) =  0 fti 0 . (82)   (n) 0 0 fto When the direction from a given atom to its neighbor is rotated around the z-axis GF, QCL and KG methods in quantum thermal transport 17 fto fti A B1 fr z y x Figure 6. Atom A is at the origin and B1 lies on the x-axis. The radial components (fr ) are along the x-axis, while the in-plane transverse components (tti ) and the out-of-plane transverse components (fto ) are parallel to the y- and z-axis, respectively. by θp to account for atom Bp, the force constant tensor is rotated as F (A,Bp) = Rz−1 (θp ) · F (A,B1) · Rz (θp ), (83) where   cos θp sin θp 0 Rz (θp ) =  − sin θp cos θp 0 . (84) 0 0 1 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 F (Ai,Bp) = Ry−1 (ϕi ) · F (A1,Bp−i+1) · Ry (ϕi ). (85) If (p − i + 1) ≤ 0 then it is replaced with (na /2 + p − i + 1), na being the number of atoms per unit cell. For B type atoms, the tensors shuld be rotated around z by π, F (Bi,Ap) = Ry−1 (ϕi ) · Rz−1 (π) · F (A1,Bp−i+1) · Rz (π) · Ry (ϕi ). (86) Figure 7. Phonon dispersion of two-dimensional graphene as obtained from density functional theory (dotted red) and 4NNFC parametrization (solid blue). The figure has been adapted from Ref. [66]. GF, QCL and KG methods in quantum thermal transport 18 Table 2. Fourth nearest neighbor force constant parameters for graphene based materials. The parameters are given in units of 104 dyn/cm = 10 N/m. (n) (n) (n) fr fti fto 1 41.8 15.2 10.2 2 7.6 -4.35 -1.08 3 -0.15 3.39 1.0 4 -0.69 -0.19 -0.55 The parameters, as optimized by Zimmermann et al. are given in Table 2. 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 [67–70] carbon nanotubes and graphene based materials [64,71–79], silicon nanowires [80–82], novel two-dimensional crystals [83, 84]. 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. 2.10.1. Defected 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 one-dimensional material, but also due to its unprecedented physical properties like extreme mechanical strength [85], ballistic electron transmission over long distances [86], and record thermal conductivity values. [87] 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. [88–93] 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. [72] Another interesting issue regards divergence of conductivity with length. Thermal conductivity GF, QCL and KG methods in quantum thermal transport 19 Figure 8. Thermal conductivity for (a) zigzag (10,0 ) and (b) armchair (10,10) CNTs for various defect types and concentrations. NEMD simulations were performed for different lengths, 200, 400 and 600 nm. The inset in (a) shows thermal conductivity for different distributions of defects, and the inset in (b) is a comparison of conductivity values for armchair and zigzag CNTs at 200, 300 and 400 K. The figure is reproduced from Ref. [72]. 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 Figure 9. Lattice thermal conductivities of defect free zigzag (10,0) and armchair (10,10) nanotubes at different temperatures as functions of length for T = 200, 300 and 400 K as obtained from NEMD simulations. The figure is reproduced from Ref. [72]. GF, QCL and KG methods in quantum thermal transport 20 −3 10 d = 0.01% −4 d = 0.02% 10 d = 0.04% d = 0.06% Mean Free Path (m) d = 0.08% −5 10 d = 0.10% 25 d = 0.15% −6 20 10 Transmission 15 −7 10 10 0 5 1 −8 0 10 0 200 2 400 600 Length (µm) −9 800 3 10 1000 0 200 400 600 800 1000 1200 1400 1600 1200 (cm-1) 1400 1600 4 ω(cm ) -1 Figure 10. Phonon transmission spectrum for CNT(10,0) with randomly distributed single vacancies (d = 1%) is given in the left panel as functions of length and frequency. Phonons display a contrasting behavior at high and low frequency regions. When the CNT is longer than 1 µm, high frequency phonons are suppressed almost completely. However, low frequency phonons remain quasi- ballistic for several µm. In the right panel, elastic mean free path, `elastic , is plotted as a function of frequency for different defect concentrations, d. `elastic scales as ω −2 , and includes details of phonon density of states. Approximately three orders of magnitude difference between the `elastic of high and low frequency modes is observed. Van Hove singularities in the DOS mark the frequencies at which the scattering rates are maximized and sharp dips in the `elastic are obtained. The figure is reproduced from Ref. [72]. (d=Ndefect /Natom ×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 more dramatic at lower concentrations. (see the right panel in Figure 10) This is in agreement with the NEMD results in Figure 8. The ballistic transmission values of low energy modes over long distances verifies the results in Figure 9 that larger conductivity values should be expected for longer tubes. 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 Wm−1 K−1 to 4 Wm−1 K−1 . Therefore, increased d is expected to enable a standardized value for the conductivity, which is still higher than 100 Wm−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 GF, QCL and KG methods in quantum thermal transport 21 (a) 12 (b) 3 ζ 10 2.5 DOS κ (nW K -1 nm -1 ) 8 0 ω 250 2 ζ (nm-1) 6 1.5 pristine C4 H 4 1 C5 H3 CH 3 Benzene C 2H 2 0.5 Biphenyl para- C3H3 Terphenyl 0 0 0 400 800 1200 1600 2000 0 250 500 750 1000 ω (cm-1 ) T (K) Figure 11. Width normalized transmission spectra for pristine and functionalized graphene (a). Details of DOS and transmission are shown in the inset. The figure is reproduced from Ref. [79]. this behavior because of ω −2 dependence of the anharmonic mean-free-path. [94] 2.10.2. Graphene with molecular antiresonances. The large surface to volume ratio of two-dimensional materials enable very effective engineering schemes for physical properties by surface functinalizations. [95–97] One possible route is to adsorbe molecules covalently. Each adsorbed molecule acts as a scatterer for phonons by creating tranmission antiresonances [98,99], 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 Ref. [79]. First, phonon transmission spectra of graphene with single molecular scatterers of different species were computed, where the force constants were obtained using Figure 12. Width normaized thermal conductance values for ensembles containing 100 adsorbants distributed over different numbers of species Nsp . The figure is reproduced from Ref. [79]. GF, QCL and KG methods in quantum thermal transport 22 the density functional tight binding (DFTB) approach. [100] CH3 , C2 H, C3 H3 , C4 H, C5 H3 , 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 11a-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 11a). Consequently, thermal conductance with a single adsorbant is also similar for all species (Figure 11b). Transmission due to multiple scatterers was calculated at the cascade scattering approximation [101, 102] using 1 1 X 1 1  = + ni − , (87) ζ ζo i ζi ζo where ζo is the transmission of the pristine system, ζi is transmission due to adsorbant species i, and ni is the number of those molecules. Keeping the total number of molecules constant and equal to 100, the change in ni is made in multiples of 10. While changing the number of species Nsp from 1 to 8, transmission and conductance values were computed for every possible combination of ni . 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 geometrical 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 nano- structuring 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. [75] Two types of graphene nanoribbons 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). [103] 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.131nW 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 GF, QCL and KG methods in quantum thermal transport 23 a precursor straight graphene nanoribbon (s-GNR) Br Br chevron type graphene nanoribbon (c-GNR) precursor Br Br b atomic distribution of isotopes precursor distribution of isotopes 7 6 Transmission 5 4 3 2 1 0 0 600 1200 1800 -1 Frequency (cm ) Figure 13. (a) Geometrical structures,(b) isotope distribution types of the studied graphene nanoribbon structures. (c) Phonon transmission spectra of pristine ribbons. The figure is reproduced from Ref. [75]. a 10 1 b 101 L = 5 µm s-GNR pristine s-GNR 10% atomic distribution s-GNR 50% atomic distribution κph /A (nW K nm ) κph /A (nW K nm ) -2 -2 0 0 s-GNR 10% precursor distribution 10 10 s-GNR 50% precursor distribution -1 -1 −1 −1 10 10 c-GNR pristine c-GNR 10% atomic distribution c-GNR 50% atomic distribution T = 300 K c-GNR 10% precursor distribution c-GNR 50% precursor distribution −2 −1 0 1 0 100 200 300 400 500 10 10 10 10 Temperature (K) Length (µm) Figure 14. Thermal conductance per cross section as functions of (a) temperature, and (b) length. The colors correspond to GNR types, straight (blue) and chevron (red). Isotope distribution types and percentages are plotted with different line types. The figure is reproduced from Ref. [75]. 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 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 GF, QCL and KG methods in quantum thermal transport 24 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, [105, 106] which can act as electronically active dopants con- trolling their electronic and phonon transport properties and, hence, its thermoelectric behavior. Using GF techniques combined with density functional tight-binding the- ory (DFTB), [107] Medrano-Sandonas et al. [104] explored the influence of ad-atoms (Hydrogen and Oxygen) and chemisorbed molecules (Hydroxyl-OH, Methyl-CH3 , and nitrophenyl-NO2 C6 H4 (NPD)) on the thermal transport properties of two possible graphene GBs, symmetric (model I) and asymmetric (model II). The model I con- sists of two graphene sheets with the same orientation (8◦ ) but in opposite directions, clockwise and counterclockwise, 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 Fig. 15). Both GBs are composed by pentagons and heptagons, and model II also presents a nonagon. [108, 109] As shown in Fig. 15, the phonon thermal conductance, κph , for both graphene GBs is lower than the corresponding one to pristine graphene for the whole range of tem- peratures 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.. [108] The phonon transmission of these grain boundaries is reduced over the whole spec- trum of frequencies by the inclusion of ad-atoms (not shown). H and O ad-atoms (a) (b) Figure 15. Temperature dependence of the phonon thermal conductance for graphene grain boundary (a) model I and (b) model II upon functionalization with ad-atoms (O and H atoms) and single molecules (OH, CH3 , and NPD). The results are divided respect the values corresponding to pristine graphene. The data has been taken from Ref. [104]. GF, QCL and KG methods in quantum thermal transport 25 (a) (b) armchair As, P Sn S zigzag Figure 16. (a) Atomistic view of the two-dimensional puckered structures. We highlighted the zigzag (ZZ) and armchair (AC) transport direction. (b) PAZ = κph−ZZ /κph−AC values as a function of the temperature for phosphorene, arsenene, and SnS monolayer. The data has been taken from Ref. [110]. 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 Fig. 15). Furthermore, OH, CH3 , 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 Fig. 16(b)). The vari- ations at low temperatures come from the strong suppression of out-of-planes modes produced by the functionalization group. For GB model I, CH3 and NPD molecules affect in a very similar magnitude its phonon transmission which is reflected in the thermal conductance values (see Fig. 15(a)). Whereas OH groups show a higher decre- ment 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 structures. Motivated by the intrinsic anisotropic physical properties found in the successfully synthesized phosphorene monolayers, [?,?] Ref. [110] analyzed how the structural anisotropy affects phonon transport properties of homo-and heteroatomic two-dimensional materials with puckered structures: arsenene and SnS monolayers (see Fig. 16(a)) and how do they compare with those in phosphorene. 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 (TA) acoustic branches have linear dispersion GF, QCL and KG methods in quantum thermal transport 26 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 longitudinal acoustic (LA) branch in phosphorene shows a group velocity of 8.35 km/s and 4.74 km/s along the Γ-X (ZZ) and Γ-Y (AC) directions, respectively, which are very close to previous DFT results. [?, ?, 116] The values for arsenene, ZZ-5.01 km/s and AC-2.71 km/s, are also in agreement with those reported by Zeraati et al. [?] SnS monolayer gives group velocities of ZZ-6.48 km/s and AC-2.14 km/s. 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 Fig. 16(b)), in good agreement to previous works. [?, ?, ?, ?, ?, ?, 115] Moreover, at low temperatures the anisotropic effect on κph (measured by PAZ = κ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, PAZ ∼ 2.5, followed by phosphorene (PAZ ∼ 1.6), which is slightly more thermal anisotropic than arsenene (PAZ ∼ 1.35). The PAZ value for phosphorene is slightly different from that obtained by employing DFT calculations, [?,114] while for (a) (b) Left Device Right contact region contact transport direction (X) B P Mo N S Setup I s=0 s = s0 s=0 Setup II s = s0 s = s0 s = s0 Figure 17. (a) Schematic representation of the partition scheme for transport calculation using Green’s function technique. The transport setups under consideration are also shown. (b) Variation of the thermal conductance as a function of the applied strain for hexagonal boron-nitride (left panel), phosphorene (center panel), and MoS2 monolayer (right panel). For comparison, the results for the standard uniaxial strain (dashed lines) are also shown. Extracted from Ref. [111]. GF, QCL and KG methods in quantum thermal transport 27 arsenene there is more discrepancy from that reported by Zeraati et al. [?], 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 two-dimensional 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 par- ticular, transport properties have been found to be very sensitive to applied strains. Thus, in Ref. [111] the GF method was applied to address the influence of strain en- gineering of the transport setup on the electron and phonon transport properties of two-dimensional materials, focusing on hexagonal boron-nitride (hBN), phosphorene, and MoS2 monolayer. Based on the partitioning scheme shown in Fig. 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), followed by phosphorene (∼ 90 W/mK). These values are close to those reported in other experimental and theoretical works. [112–116] 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 de- creases with increased applied strain, as it can be seen in Fig. 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 an slightly increment. Similar results were reported by Zhu and Ertekin [112] using non-equilibrium molecular dynamics and Boltzmann transport equation methods. In the case of phosphorene and MoS2 monolayer the influence of the strain is quali- tatively 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 Fig. 17(b) (center panel). The increment in strain reduces the inter- layer distance and, hence, the strength of the out-of-plane bonds increases, becoming higher than the in-plane ones. As a result of this, the transmission probability in- creases at low frequencies (in-plane modes) and so does the thermal transport, as it was also observed by Z.-Y. Ong et al.. [114] For MoS2 monolayer, despite the rela- tively 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 GF, QCL and KG methods in quantum thermal transport 28 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 Fig. 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-of-plane 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 Fig. 17(b)). However, for MoS2 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 transmission 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 (a) 40 Transmission (nm−2) 30 20 10 0 0 50 100 150 200 250 Energy (meV) (b) 50 Transmission (nm−2) 40 30 20 10 0 0 50 100 150 200 250 Energy (meV) (c) 40 Transmission (nm−2) 30 20 10 0 0 50 100 150 200 250 Energy (meV) Figure 18. Transmission per cross section area for graphene over a substrate step. Transmission with varying the step height, hs , for fixed coupling strength LJ = 40 meV (b), and for various coupling strength with a fixed step height hs = 10 Å. The figure is reproduced from Ref. [76]. GF, QCL and KG methods in quantum thermal transport 29 (a) (b) 1 1 0.9 0.9 0.8 T = 300 K 0.8 T = 500 K 0.7 0.7 T = 800 K 0.6 0.6 0.5 0.5 0 20 40 0 50 100 150 (meV) 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 reproduced from Ref. [76]. 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, [117–121] 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 [76]    σ 12  σ 6 VLJ = 4LJ − , (88) r r 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 VLJ 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-of-plane acousitc modes around zero energy. The size of the gap is proportional to LJ 1/2 (see Figure 18a). Imposing a step at the substrate, graphene is bent. The local curvature depends on the step height hs and the interaction strength LJ . Namely, the radius of curvature decreases with increasing hs and stronger LJ . Smaller radius of curvature couples in-plane and out-of-plane modes stronger and a a significant reduction of conductance takes place (Figure 18b-c). Transmission across kinked graphene is also shown in Figure 18c, where the kink is induced by adsorption of lines of hydrogen atoms. [121, 122] For hs = 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 hs , 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 hs values (Figure 19a). 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 19b). Assuming that resistances due to GF, QCL and KG methods in quantum thermal transport 30 distant steps are additive, it should be possible to reduce thermal conductance by at least an order of magnitude using stepped substrates. 3. 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 molecular dynamics (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 [123,124]. This approach has for example been used to examine non-linear thermal effects such as thermal rectification in nano-junctions [125, 126]. 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, Fc , 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, g α (ω) = [(ω + iδ)2 1 − Φαα ]−1 . (89) 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 [127–129]. 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, uC , coupled into the lead by the spring, ΦαC , and propagated, Z t uα (t) = g α (t)u̇α (0) + ġ α (t)uα (0) + g α (t − t0 )ΦαC uC (t0 )dt0 . (90) 0 C This may be inserted into the EOM for u . The resulting equation is the generalized Langevin equation(GLE) for the central coordinates uC (mass-scaled), Z t üC (t) = Fc (x) − Σ(t − t0 )uC (t0 )dt0 + ξ(t) , (91) 0 where we have used the retarded self-energy, X Σ(t) = ΦCα g α (t) ΦαC = ΣL (t) + ΣR (t) , (92) α=L,R GF, QCL and KG methods in quantum thermal transport 31 and introduced the propagated initial conditions of the reservoir as a driving force, ξ, in the central region, X ξ(t) = ΦCα g α (t)u̇α (0) + ΦCα ġ α (t)uα (0) = ξ L (t) + ξ R (t) . (93) α=L,R 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, Z  α  α iωt α α Γ (ω) S (ω) = e hξ (t)ξ (0)i dt = kB Tα , (94) ω where Γα = i[Σα (ω)−Σα ∗ (ω)]. This also establish 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 [130–132]. To obtain the heat flux to a reservoir one calculate the (time averaged) power between the central region and a reservoir as the work done by the force originating from central region, J α = Ė α = hu˙α · ΦαC uC it (95) The 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, J α = −hF α · u˙α it (96) with, Z t F α (t) = − Σα (t − t0 )uC (t0 )dt0 + ξ α (t) . (97) 0 However at low temperature 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 (cf. Equation 50) for harmonic, steady-state transport [127]. 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 [128,133] one may adopt a quasi-classical [134] 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, Z 1 S α (ω) = eiωt hξˆα (t)ξˆα (0) + ξˆα (0)ξˆα (t)i dt 2   α 1 = ~ nB (~ω/kB T ) + Γα (ω) . (98) 2 GF, QCL and KG methods in quantum thermal transport 32 This approach yield the ballistic Landauer regime for low temperatures and the classical, diffusive regime at high temperature. 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 [135–137]. 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 [133]. A damping may be applied to limit the length of the correlation in time [138, 139] Alternatively the noise can be generated “on the fly” [140], or the noise and friction can be represented using auxilary degrees of freedom [141, 142]. 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 [138], nanotubes, ribbons with edge roughness [143], and constrictions [131]. For example for a 5-dimer wide graphene constriction it was found that although the local temperature profile agreed well for classical (Eq.94) and quasi-classical MD (Eq. 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 a) b) c) d) Figure 20. Local temperatures in a 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 [144]. (a,b) Calculated using the harmonic approximation, while (c,d) including anharmonicity. In b,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. From [145]. GF, QCL and KG methods in quantum thermal transport 33 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 [146, 147]. 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 Fig. 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 Fig. 20(a,b) and Fig. 20(c,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] 4. 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. [60, 61, 63] In the Kubo-Greenwood approach only the harmonic part of the Hamiltonian (cf. Equation 1) is taken into account. Following Allen and Feldman [148], conductivity σ along x-direction is defined as Zβ Z∞ Ω σ= lim dλ dt ei(ω̃+iδ)t hJ x (−i~λ)J x (t)i, (99) T ω̃,δ→0 0 0 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, 1 X Jx = (Xi − Xj )Φij xi ẋj . (100) 2Ω i,j Here, Xi is the equilibrium position of the atom, which the ith degree of freedom belongs, xi 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 Jmn x a†m an , P where r r  x i~ ωm ωn Jmn = − + hm|[X, Φ]|ni (101) 4Ω ωn ωm with X being the diagonal matrix of equilibrium positions and |ni is the nth egenstate of Φ. It was also shown by Allen and Feldman [148] that the conductivity can be GF, QCL and KG methods in quantum thermal transport 34 written as πΩ X ∂fB x x σ= J J δ(ωm − ωn ). (102) ~T m,n ∂ωm mn nm The above relation can also be expressed as Z∞ √ √ o π ~ ∂fB n σ=− dω Tr [X, Φ]δ(ω − Φ)[X, Φ]δ(ω − Φ) , (103) 4Ω ω ∂T 0 √ P where Φ = n ωn |nihn| 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 Φ|ni = ωn2 |ni, √ and the matrices Φ and Φ have simultaneous eigenstates. Defining Vx = −i[X, Φ], thermal conductance of a one-dimensional system is obtained as Z∞ √ √ o π ∂fB n κ= 2 dω ~ω Tr Vx δ(ω − Φ)Vx δ(ω − Φ) . (104) L ∂T 0 Recalling the definition of conductance in terms of Green functions (see Equation 51), one obtains the transmission function as 2π 2 n √ √ o ζ(ω) = 2 Tr Vx δ(ω − Φ)Vx δ(ω − Φ) , (105) L which is of the same form with the electronic transmission function as derived from the Kubo-Greenwood formula [149, 150] 2π 2 ~2 ζel (E) = Tr {Vx δ(E − Hel )Vx δ(E − Hel )} , (106) L2 E being energy and Hel the electronic Hamiltoian. In the diffusion regime, transmission function is expressed as 2ωπ ζ(ω) = 2 Tr δ(ω 2 − Φ) Dmax (ω),  (107) L where Dmax (ω) is the maximum value of time-dependent diffusion coefficient D(ω, t). It is defined in terms of the mean-square displacement as χ2 (ω, t) D(ω, t) = , (108) t where χ2 (ω, t) = h(X(t) − X(0))2 iω . Computation of D(ω, t) makes it possible to deduce the transport coefficients such as Dmax (ω), 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  2 Tr (X(t) − X(0))2 δ(ω 2 − Φ) χ (ω, t) = , (109) Tr {δ(ω 2 − Φ)} Tr [X, U (t)]† δ(ω 2 − Φ)[X, U (t)]  = , (110) Tr {δ(ω 2 − Φ)} where U (t) is the time evolution operator. The numerator can be computed through the average of random phase states as N hψ|[X, U (t)]† δ(ω 2 − Φ)[X, U (t)]|ψi, 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)]|ψi. 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. GF, QCL and KG methods in quantum thermal transport 35 4.1. Lanczos tridiagonalization and the continued fraction methods Computationally efficient algorithms are required in order to be able to simulate large systems. The Lanczos method [151] 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 i, a basis set {|ψn i} is built recursively as an = hψn |Φ|ψn i (111) 0 |ψn+1 i = Φ|ψn i − an |ψn i − bn−1 |ψn−1 i, (112) q bn 0 = hψn+1 0 |ψn+1 i, (113) 1 0 |ψn+1 i = |ψ i, (114) bn n+1 where n > 1 and b0 = 0. The recursion coefficients an and bn are the diagonal and the off-diagonal elements of the dynamical matrix in this basis, respectively.   a1 b1 0 0 · · ·  b1 a2 b2 0     0 b2 a3 b3 Φ= (115)    0 0 b3 a4  ..   .. . . One chooses the seed vector as a random phase state N 1 X iϕj |ψ1 i = √ e |uj i, (116) N j=1 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 1 1 hψ1 |δ(ω 2 − Φ)|ψ1 i = − lim Imhψ1 | |ψ1 i, (117) π δ→0 (ω + iδ)2 − Φ PDOS (arbitrary unit) an bn 0 20 40 60 80 100 Recursion number n 0 200 400 600 800 1000 1200 1400 1600 −1 ω (cm ) Figure 21. Phonon DOS of two-dimensional graphene as obtained using the Lanczos continued fraction method. Converge of the recursion coefficients are shown in the inset. Figure is reproduced from Li et al., Phys. Rev. B 83, 155416 (2011). GF, QCL and KG methods in quantum thermal transport 36 one can show that 1 1 hψ1 | |ψ1 i = b21 (118). (ω + iδ)2 − Φ (ω + iδ)2 − a1 − b2 2 (ω+iδ)2 −a2 − b2 3 (ω+iδ)2 −a3 − .. . The continued fraction in Equation 118 is terminated after the coefficients an , bn converge at step N . One names the continued fraction as G1 and defines G2 with 1 G1 = , (119) (ω + iδ) − a1 − b21 G2 2 and similarly 1 Gn = . (120) (ω + iδ)2 − an − b2n Gn+1 For the converged recursion coefficients aN and bN , the below equality should be satisfied 1 1 GN +1 = = (121) . (ω + iδ)2 − aN − b2N GN +2 (ω + iδ)2 − aN − b2N GN +1 Consequently one finds p (ω + iδ)2 − aN − i (2bN )2 − ((ω + iδ)2 − aN )2 GN +1 = . (122) 2b2N As an example, the recursion coefficients and the DOS of two-dimensional graphene are shown in Figure 21. 4.2. 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 Φ0 = (Φ − 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 √ ∞ X U (t) = ei Φt = cn (t)Qn (Φ0 ) (123) n=0 with Z1 √ 2 1 0 cn (t) = dx0 √ e−it a+2bx Qn (x0 ) (124) π(1 + δn,0 ) 1−x 02 −1 P Hence the commutator takes the form [X, U (t)]|ψi= n cn (t)[X, Qn (Φ)]|ψi. Cheby- shev polynomials fulfill the recurrence relation Qn+1 (Φ0 ) = 2Φ0 Qn (Φ0 ) − Qn−1 (Φ0 ), (125) GF, QCL and KG methods in quantum thermal transport 37 where Q0 (Φ0 ) = 1 and Q1 (Φ0 ) = Φ0 . [152] Consequently, the commutator takes the following form, [X, Qn+1 (Φ0 )] = 2[X, Φ0 Qn (Φ0 )] − [X, Qn−1 (Φ0 )]. (126) Denoting |αn i = Qn (Φ0 )|ψi, (127) |βn i = [X, Qn (Φ0 )]|ψi, (128) one writes |βn+1 i = 2Φ0 |βn i + 2[X, Φ0 ]|αn i − |βn−1 i, (129) 0 with |β0 i = 0 and |β1 i = [X, Φ ]|ψi by definition. For obtaining |βn i, one needs |αn i, which satisfies |αn+1 i = Qn+1 (Φ0 )|ψi, (130) 0 = 2Φ |αn i − |αn−1 i, (131) 0 with |α0 i = |ψi and |α1 i = Φ |ψi. As a result, the computation of time evolution for a finite ∆t is equivalent to obtaining the following, Npoly X U (∆t)|ψi = cn (∆t)|αn i, (132) n=0 Npoly X [X, U (∆t)]|ψi = cn (∆t)|βn i. (133) n=0 The upper limit Npoly is to be chosen depending on the evolution step and the band width. Once U (m∆t)|ψi and [X, U (m∆t)]|ψi are obtained, one proceeds as U ((m + 1)∆t)|ψi = [X, U (∆t)]U (m∆t)|ψi, (134) and [X, U ((m + 1)∆t)] |ψi = [X, U (∆t)] U (m∆t) |ψi + U (∆t) [X, U (m∆t)] |ψi. (135) That is, the time evolution is computed step by step using the expansion coefficient cn defined in Equation 129. The expansion coefficients are computed numerically on a Chebyshev-Gauss grid of N points π(k + 21 )   0 xk = cos , k = 0, . . . , N − 1 (136) N using the quadrature formula Z1 N −1 f (x0 ) π X dx0 √ = f (x0k ). (137) 1 − x02 N k=0 −1 4.3. 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 [153–155]. In view of thermoelectric applications, lower thermal conductivity is sought, while preserving electronic transport, which remains a very challenging quest. Here we GF, QCL and KG methods in quantum thermal transport 38 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 materials. 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 carbon nanotubes (CNTs) [61, 156], graphene- based materials [157], and graphene-boron nitride hybrid two-dimensional membranes (including polycrystalline structures) [60, 158], obtained with Kubo-type and Green function methods. It is known that isotopic or Anderson-type disorder can reduce thermal conductivity in quasi-one-dimensional structures of sp2 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 [159]. 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, graphene nanoribbons (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 [60]. Finally two-dimensional heterostructures made from mixing graphene with disordered (isotope impurities) GNRs [160], or graphene-boron nitride membranes and polycrystalline mixture [161] 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 [159]. Figure 22 (left panel) illustrates the disordered nanotube while Fig. 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 Figure 22. Right panel: frequency-dependent MFP in CNT with isotopic disorder (10%) with the GF method, the KG method and the analytical formula (inset). Left panel: sketch of a disordered carbon nanotube with isotopic disorder. Figure is reproduced from Li et al., Phys. Rev. B 83, 155416 (2011). GF, QCL and KG methods in quantum thermal transport 39 Figure 23. Left panel: Ball-and-stick representations of GNRs. Right panel: elastic MFP for ZGNR with Nz =80 and for the Na =138 (widths in order of 15 nm) with disorder density of 10%. Figure is reproduced from Li et al., Phys. Rev. B 83, 155416 (2011). value characterizes diffusive transport. The results compare very well with those obtained with the GF approach [159], as well as with the analytical result `(ω) = (12aNuc Nch (ω))/(π 2 f | ∆M 2 2 2 M | ρuc (ω)ω ) where a is the length of the lattice vector in the translational direction, Nuc 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 [61]. 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 1500cm−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- of-plane) acoustic modes, whereas phonon lifetimes of optical phonon modes were also found to be considerably smaller [162] (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 [162]. 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. Fig. 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 Nz , whereas the edge disorder is dictated GF, QCL and KG methods in quantum thermal transport 40 by the density of edge defects (removed carbon atoms at the edges) and is chosen to be 10% of edge C-atoms (Fig. 23). One observes large differences in the MFPs between ZGNR and AGNR. Indeed if MFP remains well below 200 nm for frequencies below 800cm−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-1100cm−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 [159]. The thermal conductances for AGNR and ZGNR are shown in Fig. 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 Fig. 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 Fig. 25. We first report evolution of the wave packet dynamics for different frequencies for ZGNR with 10% edge disorder, through the time-dependent 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 40 pristine 35 L=100 nm L=500 nm 30 L=1 µm 25 κ (nW/K) L=2 µm 20 15 10 5 0 0 100 200 300 400 500 Temperature (K) Figure 24. Thermal conductance for the ZGNR (Nz =80, solid lines) and the AGNR (Na = 138, dashed lines) with edge disorder of 10% and various ribbon lengths. Figure is reproduced from Ref. [61]. GF, QCL and KG methods in quantum thermal transport 41 800 Diffusion Coefficient (nm /ps) ZGNR(20)10% ZGNR(40)10% 2 1200 700 ω = 1130 cm -1 ZGNR(80)10% ω = 250 cm -1 Mean free path (nm) ZGNR(80)15% 800 ω = 50 cm -1 600 400 500 0 0 500 1000 1500 2000 t (ps) 400 300 200 100 0 70 200 400 600 800 1000 1200 1400 1600 -1 (cm ) Figure 25. Left panel: short portion of an edge-disordered ZGNR (width 17 nm, length L=50 nm) with 10% edge disorder . Right panel: MFP for ZGNR of widths 4.26 nm, 8.52 nm, and 17 nm (the length of ZGNR is 985 nm), respectively with disorder density of 10%, and also for the 17 nm wide ZGNR with 15% disorder. Inset: time-dependent diffusion coefficients for three chosen frequencies. Figure is reproduced from Ref. [63]. long localization lengths. The saturation of D(ω, t) to a maximum value dictates the value of the elastic MFP. In Fig. 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 [157]. 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 Wm−1 K−1 and σ = 671 Wm−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 Fig.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 two-dimensional material growth) will be analyzed afterwards. The calculated transport MFPs are plotted in Fig.26 (bottom panels). The density of states (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 [163]. Accordingly, the out-of-plane modes are significantly scatter 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. Fig.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 GF, QCL and KG methods in quantum thermal transport 42 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 [164, 165]. 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 [154, 166, 167]. 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. [158] For small average grain size, the minimum value of thermal conductivity occurs near 70% hBN agreeing with prior KG results [60]. 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. 4 10 4 80% 10 11 3 60% B 10 10 3 10 50% B Mean Free Path (nm) Mean Free Path (nm) 11 10 10 2 40% 2 B: B = 0.801:0.199 10 20% 1 1 10 (a) 10 0 (a) 0 10 4 10 10 4 10 3 10 3 10 2 2 10 10 1 1 10 (b) 10 (b) 0 0 10 10 200 400 600 800 1000 1200 1400 1600 200 400 600 800 1000 1200 1400 1600 w (cm-1 ) w (cm -1) Figure 26. Top panel: Sketch of a graphene-BN heterostructures ( BN clusters with 2-nm diameters are inserted inside graphene). Bottom-left panel: Mean free paths for different BN islands concentrations for (a) out-of-plane and (b) in-plane modes. Bottom-right panel: same for MFPs for different B-isotopic compositions at a 50% density. Figure is reproduced from Ref. [60]. GF, QCL and KG methods in quantum thermal transport 43 5. 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 is Section 2.7. These algorithms are useful also for two-dimensional systems, when periodic boundary conditions are applicable. With the use of recursion algorithms, low dimensional systems with more than 106 atoms can be treated numerically 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 methods scale linearly, whereas most of the methods including GFs 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 hence 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. Acknowledgment HS acknowledges support from TÜBİTAK (115F445, 117F480) and Bilim Akademisi -the Science Academy, Turkey. The authors would like to thank Profs. Lü and Wang and Dr. Gunst for useful discussions on the Langevin approach. GF, QCL and KG methods in quantum thermal transport 44 6. References [1] Schwab K, Henriksen E A, Worlock J M and Roukes M L 2000 Nature 404 974–977 [[4]], [[11]] [2] Rego L G C and Kirczenow G 1998 Physical Review Letters 81 232–235 [[4]], [[11]] [3] Molenkamp L W, Gravier T, van Houten H, Buijk O J A, Mabesoone M A A and Foxon C T 1992 Physical Review Letters 68 3765–3768 [[4]], [[11]] [4] Rego L G C and Kirczenow G 1999 Physical Review B 59 13080–13086 [[4]] [5] Li N, Ren J, Wang L, Zhang G, Hänggi P and Li B 2012 Reviews of Modern Physics 84 1045–1066 [[4]] [6] Cahill D G, Ford W K, Goodson K E, Mahan G D, Majumdar A, Maris H J, Merlin R and Phillpot S R 2003 Journal of Applied Physics 93 793–818 [[4]] [7] Cahill D G, Braun P V, Chen G, Clarke D R, Fan S, Goodson K E, Keblinski P, King W P, Mahan G D, Majumdar A, Maris H J, Phillpot S R, Pop E and Shi L 2014 Applied Physics Reviews 1 011305 [[4]] [8] Volz S, Ordonez-Miranda J, Shchepetov A, Prunnila M, Ahopelto J, Pezeril T, Vaudel G, Gusev V, Ruello P, Weig E M, Schubert M, Hettich M, Grossman M, Dekorsy T, Alzina F, Graczykowski B, Chavez-Angel E, Sebastian Reparaz J, Wagner M R, Sotomayor-Torres C M, Xiong S, Neogi S and Donadio D 2016 The European Physical Journal B 89 15 [[4]] [9] Zhou Y, Fan Z, Qin G, Yang J Y, Ouyang T and Hu M 2018 ACS Omega 3 3278–3284 [[4]] [10] Tewary V K and Yang B 2009 Physical Review B 79 125416 [[4]] [11] Broido D A, Malorny M, Birner G, Mingo N and Stewart D A 2007 Applied Physics Letters 91 231922 [[4]] [12] Wang J S, Wang J and Lü J T 2008 The European Physical Journal B 62 381–404 [[5]] [13] Wang J S, Agarwalla B K, Li H and Thingna J 2014 Frontiers of Physics 9 673–697 [[5]] [14] Mingo N 2009 Green’s function methods for phonon transport through nano-contacts Thermal Nanosystems and Nanomaterials (Topics in Applied Physics no 118) ed Volz S (Springer Berlin Heidelberg) pp 63–94 ISBN 978-3-642-04257-7 978-3-642-04258-4 [[5]] [15] Chen X, Liu Y and Duan W 2018 Small Methods 2 1700343 [[5]] [16] Lü J T, Hu B Z, Hedegard P and Brandbyge M Progress in Surface Science URL https: //doi.org/10.1016/j.progsurf.2018.07.002 [[5]], [[33]] [17] Mahan G D 2000 Many-Particle Physics (Springer) [[5]], [[6]] [18] Schwinger J 1961 Journal of Mathematical Physics 2 407–432 [[5]] [19] Keldysh L V 1965 Sov. Phys. JETP 20 1018–1026 [[5]] [20] Kadanoff L and Baym G 1962 Quantum statistical mechanics: Green’s function methods in equilibrium and nonequilibrium problems Frontiers in physics (W.A. Benjamin) [[5]] [21] Meir Y and Wingreen N S 1992 Physical Review Letters 68 2512–2515 [[5]] [22] Ozpineci A and Ciraci S 2001 Physical Review B 63 [[5]], [[11]], [[18]] [23] Mingo N and Yang L 2003 Physical Review B 68 [[5]], [[11]] [24] Yamamoto T and Watanabe K 2006 Physical Review Letters 96 [[5]], [[11]] [25] Craig R A 1968 Journal of Mathematical Physics 9 605–611 [[8]] [26] Xie Z X, Li K M, Tang L M, Pan C N and Chen K Q 2012 Applied Physics Letters 100 183110 [[12]] [27] Peng X F, Xiong C, Wang X J, Chen L Q, Luo Y F and Li J B 2013 Computational Materials Science 77 440–444 [[12]] [28] Zhang L, Wang J S and Li B 2010 Physical Review B 81 [[12]] [29] Ouyang T, Chen Y, Xie Y, Wei X L, Yang K, Yang P and Zhong J 2010 Physical Review B 82 [[12]] [30] Ding X and Ming Y 2014 Chinese Physics Letters 31 046601 [[12]] [31] Ming Y, Wang Z X, Ding Z J and Li H M 2010 New Journal of Physics 12 103041 [[12]] [32] Rich M and Visscher W M 1975 Physical Review B 11 2164–2170 [[12]] [33] Bolsterli M, Rich M and Visscher W M 1970 Physical Review A 1 1086–1088 [[12]] [34] Bonetto F, Lebowitz J L and Lukkarinen J 2004 Journal of Statistical Physics 116 783–813 [[12]] [35] Büttiker M 1986 Physical Review B 33 3020–3026 [[12]] [36] Büttiker M 1986 Physical Review Letters 57 1761–1764 [[12]] [37] Dhar A and Roy D 2006 Journal of Statistical Physics 125 801–820 [[12]] [38] Roy D 2008 Physical Review E 77 [[12]] [39] Dhar A 2008 Advances in Physics 57 457–537 [[12]] [40] Visscher W M and Rich M 1975 Physical Review A 12 675–680 [[12]] [41] Neto A F, Lemos H C F and Pereira E 2007 Physical Review E 76 031116 [[12]] [42] Pereira E 2010 Physics Letters A 374 1933–1937 [[12]] GF, QCL and KG methods in quantum thermal transport 45 [43] Pereira E, Lemos H C F and Ávila R R 2011 Physical Review E 84 061135 [[12]] [44] Wu L A and Segal D 2011 Physical Review E 83 051114 [[12]] [45] Bandyopadhyay M and Segal D 2011 Physical Review E 84 011151 [[12]] [46] Sääskilahti K, Oksanen J and Tulkki J 2013 Physical Review E 88 012128 [[12]] [47] Sancho M P L, Sancho J M L, Sancho J M L and Rubio J 1985 Journal of Physics F: Metal Physics 15 851–858 [[14]] [48] Guinea F, Tejedor C, Flores F and Louis E 1983 Physical Review B 28 4397–4402 [[14]] [49] Verzijl C J O, Seldenthuis J S and Thijssen J M 2013 The Journal of Chemical Physics 138 094102 [[15]] [50] Kurt G and Sevinçli H 2018 in preparation [[15]] [51] Müller W, Schiller R and Nolting W 2000 The European Physical Journal B-Condensed Matter and Complex Systems 16 705–718 [[15]] [52] Saito R, Dresselhaus G and Dresselhaus M S 1998 Physical Properties of Carbon Nanotubes (Imperial College Press) [[16]] [53] Yan Y, Liang Q F, Zhao H, Wu C Q and Li B 2012 Physics Letters A 376 2425–2429 [[16]] [54] Yan Y, Liang Q F, Zhao H and Wu C Q 2012 Physics Letters A 376 1154–1158 [[16]] [55] Chang P H and Nikolić B K 2012 Physical Review B 86 [[16]] [56] Yan Y and Zhao H 2012 Journal of Applied Physics 111 113531 [[16]] [57] Pan C N, Xie Z X, Tang L M and Chen K Q 2012 Applied Physics Letters 101 103115 [[16]] [58] Xie Z X, Tang L M, Pan C N, Li K M, Chen K Q and Duan W 2012 Applied Physics Letters 100 073105 [[16]] [59] Xie Z X, Li K M, Tang L M, Pan C N and Chen K Q 2012 Applied Physics Letters 100 183110 [[16]] [60] Sevinçli H, Li W, Mingo N, Cuniberti G and Roche S 2011 Physical Review B 84 205444 [[16]], [[33]], [[38]], [[42]], [[43]] [61] Li W, Sevinçli H, Roche S and Cuniberti G 2011 Physical Review B 83 155416 [[16]], [[33]], [[38]], [[39]], [[40]] [62] Aksamija Z and Knezevic I 2011 Applied Physics Letters 98 141919 [[16]] [63] Li W, Sevinçli H, Cuniberti G and Roche S 2010 Physical Review B 82 041410 [[16]], [[33]], [[41]] [64] Sevinçli H and Cuniberti G 2010 Physical Review B 81 113401 [[16]], [[18]] [65] Chen Y, Jayasekera T, Calzolari A, Kim K W and Nardelli M B 2010 Journal of Physics: Condensed Matter 22 372202 [[16]] [66] Zimmermann J, Pavone P and Cuniberti G 2008 Physical Review B 78 [[16]], [[17]] [67] Nikolić B K, Saha K K, Markussen T and Thygesen K S 2012 Journal of Computational Electronics 11 78–92 [[18]] [68] Saha K K, Markussen T, Thygesen K S and Nikolić B K 2011 Physical Review B 84 041412 [[18]] [69] Markussen T 2013 The Journal of Chemical Physics 139 244101 [[18]] [70] Nozaki D, Sevinçli H, Li W, Gutiérrez R and Cuniberti G 2010 Physical Review B 81 [[18]] [71] Xu Y, Wang J S, Duan W, Gu B L and Li B 2008 Physical Review B 78 224303 [[18]] [72] Sevik C, Sevinçli H, Cuniberti G and Çağın T 2011 Nano Letters 11 4971–4977 [[18]], [[19]], [[20]] [73] Gunst T, Markussen T, Jauho A P and Brandbyge M 2011 Physical Review B 84 155449 [[18]] [74] Sevinçli H, Lehmann T, Ryndyk D A and Cuniberti G 2013 Journal of Computational Electronics 12 685–691 [[18]] [75] Sevinçli H, Sevik C, Çağın T and Cuniberti G 2013 Scientific Reports 3 1228 [[18]], [[22]], [[23]] [76] Sevinçli H and Brandbyge M 2014 Applied Physics Letters 105 153108 [[18]], [[28]], [[29]], [[30]] [77] Sevinçli H and Sevik C 2014 Applied Physics Letters 105 223108 [[18]] [78] Sandonas L M, Sevinçli H, Gutierrez R and Cuniberti G 2018 Advanced Science 5 1700365 [[18]] [79] Sevim K and Sevinçli H 2018 Carbon 140 603–609 [[18]], [[21]] [80] Markussen T, Jauho A P and Brandbyge M 2008 Nano Letters 8 3771–3775 [[18]] [81] Markussen T, Jauho A P and Brandbyge M 2009 Physical Review Letters 103 055502 [[18]] [82] Markussen T, Jauho A P and Brandbyge M 2009 Physical Review B 79 035415 [[18]] [83] Sevik C and Sevinçli H 2016 Nanotechnology 27 355705 [[18]] [84] Sevinçli H 2017 Nano Letters 17 2589–2595 [[18]] [85] Yakobson B I, Brabec C J and Bernholc J 1996 Physical Review Letters 76 2511–2514 [[18]] [86] White C T and Todorov T N 1998 Nature 393 240–242 [[18]] [87] Che J, Çagin T and III W A G 2000 Nanotechnology 11 65 [[18]] GF, QCL and KG methods in quantum thermal transport 46 [88] Li Q, Liu C, Wang X and Fan S 2009 Nanotechnology 20 145702 [[18]] [89] Pettes M T and Shi L 2009 Advanced Functional Materials 19 3918–3925 [[18]] [90] Hone J, Whitney M, Piskoti C and Zettl A 1999 Physical Review B 59 R2514–R2516 [[18]] [91] Pop E, Mann D, Wang Q, Goodson K and Dai H 2006 Nano Letters 6 96–100 [[18]] [92] Pop E, Mann D, Cao J, Wang Q, Goodson K and Dai H 2005 Physical Review Letters 95 [[18]] [93] Yu C, Shi L, Yao Z, Li D and Majumdar A 2005 Nano Letters 5 1842–1846 [[18]] [94] Klemens P and Pedraza D 1994 Carbon 32 735–741 [[21]] [95] Balandin A A 2011 Nature Materials 10 569–581 [[21]] [96] Rocha C G, Rümmeli M H, Ibrahim I, Sevincli H, Börrnert F, Kunstmann J, Bachmatiuk A, Pötschke M, Li W and Makharza S A M 2011 Tailoring the physical properties of graphene (CRC Press) [[21]] [97] Rümmeli M H, Rocha C G, Ortmann F, Ibrahim I, Sevincli H, Börrnert F, Kunstmann J, Bachmatiuk A, Pötschke M, Shiraishi M, Meyyappan M, Büchner B, Roche S and Cuniberti G 2011 Advanced Materials 23 4471–4490 [[21]] [98] Sevinçli H, Mukhopadhyay S, Senger R T and Ciraci S 2007 Physical Review B 76 [[21]] [99] Nozaki D, Sevinçli H, Avdoshenko S M, Gutierrez R and Cuniberti G 2013 Physical Chemistry Chemical Physics 15 13951 [[21]] [100] Aradi B, Hourahine B and Frauenheim T 2007 The Journal of Physical Chemistry A 111 5678–5684 [[22]] [101] Stewart D A, Savić I and Mingo N 2009 Nano Letters 9 81–84 [[22]] [102] Wang J, Li L and Wang J S 2011 Applied Physics Letters 99 091905 [[22]] [103] Cai J, Ruffieux P, Jaafar R, Bieri M, Braun T, Blankenburg S, Muoth M, Seitsonen A P, Saleh M, Feng X, Müllen K and Fasel R 2010 Nature 466 [[22]] [104] Medrano Sandonas L, Gutierrez R, Pecchia A, Dianat A and Cuniberti G 2015 J. Self-Assem. Mol. Electron. 3(1) 1–20 [[24]] [105] Huang P, Zhu H, Jing L, Zhao Y and Gao X 2011 ACS Nano 5 7945–7949 [[24]] [106] Zhang H, Lee G, Gong C, Colombo L and Cho K 2014 The Journal of Physical Chemistry C 118 2338–2343 [[24]] [107] Frauenheim T, Seifert G, Elsterner M, Hajnal Z, Jungnickel G, Porezag D, Suhai S and Scholz R 2000 Phys. Status Solidi B 217 41–62 ISSN 1521-3951 [[24]] [108] Meyer J C, Kisielowski C, Erni R, Rossell M D, Crommie M F and Zettl A 2008 Nano Letters 8 3582–3586 [[24]] [109] Huang P Y, Ruiz-Vargas C S, van der Zande A M, Whitney W S, Levendorf M P, Kevek J W, Garg S, Alden J S, Hustedt C J, Zhu Y, Park J, McEuen P L and Muller D A 2011 Nature 469 389–392 [[24]] [110] Medrano Sandonas L, Teich D, Gutierrez R, Lorenz T, Pecchia A, Seifert G and Cuniberti G 2016 J. Phys. Chem. C 120 18841–18849 [[25]] [111] Medrano Sandonas L, Gutierrez R, Pecchia A, Seifert G and Cuniberti G 2017 Phys. Chem. Chem. Phys. 19(2) 1487–1495 [[26]], [[27]] [112] Zhu T and Ertekin E 2015 Phys. Rev. B 91(20) 205429 [[27]] [113] Li W, Carrete J and Mingo N 2013 Appl. Phys. Lett. 103 253103 [[27]] [114] Ong Z Y, Cai Y, Zhang G and Zhang Y W 2014 J. Phys. Chem. C 118 25272–25277 [[27]], [[28]] [115] Jain A and McGaughey A J H 2015 Sci. Rep. 5 8501 [[27]] [116] Liu T H and Chang C C 2015 Nanoscale 7(24) 10648–10654 [[26]], [[27]] [117] Bunch J S and Dunn M L 2012 Solid State Communications 152 1359–1364 [[28]] [118] Lee J K, Yamazaki S, Yun H, Park J, Kennedy G P, Kim G T, Pietzsch O, Wiesendanger R, Lee S, Hong S, Dettlaff-Weglikowska U and Roth S 2013 Nano Letters 13 3494–3500 [[28]] [119] Hicks J, Tejeda A, Taleb-Ibrahimi A, Nevius M S, Wang F, Shepperd K, Palmer J, Bertran F, Le Fèvre P, Kunc J, de Heer W A, Berger C and Conrad E H 2013 Nature Physics 9 49–54 [[28]] [120] Low T, Perebeinos V, Tersoff J and Avouris P 2012 Physical Review Letters 108 096601 [[28]] [121] Rasmussen J T, Gunst T, Bøggild P, Jauho A P and Brandbyge M 2013 Beilstein Journal of Nanotechnology 4 103–110 [[28]], [[29]] [122] Özbal G, Falkenberg J T, Brandbyge M, Senger R T and Sevinçli H 2018 Phys. Rev. Materials 2(7) 073406 [[29]] [123] Cummings P and Evans D 1992 INDUSTRIAL & ENGINEERING CHEMISTRY RESEARCH 31 1237–1252 ISSN 0888-5885 [[30]] [124] Chantrenne P, Joulain K and Lacroix D 2009 Nanostructures THERMAL NANOSYSTEMS AND NANOMATERIALS (Topics in Applied Physics vol 118) ed Volz, S (Springer Berlin / Heidelberg) pp 17–62 ISBN 978-3-642-04257-7 [[30]] GF, QCL and KG methods in quantum thermal transport 47 [125] Wu G and Li B 2007 Phys. Rev. B 76 085424 ISSN 1098-0121 [[30]] [126] Cartoixà X, Colombo L and Rurali R 2015 Nano Lett. 15 8255–8259 ISSN 1530-6984 [[30]] [127] Dhar A and Roy D 2006 J. Stat. Phys. 125 801–820 ISSN 0022-4715 [[30]], [[32]] [128] Wang J S, Wang J S and Lü J T 2008 Eur. Phys. J. B - Condens. Matter Complex Syst. 62 381–404 ISSN 1434-6028 [[30]], [[32]] [129] Kantorovich L 2008 Phys. Rev. B 78 94304 ISSN 1098-0121 [[30]] [130] Karpov E G, Park H S and Liu W K 2007 Int. J. Numer. Methods Eng. 70 351–378 ISSN 00295981 [[31]] [131] Sääskilahti K, Oksanen J and Tulkki J 2013 Phys. Rev. E 88 012128 [[31]], [[32]] [132] Ness H, Genina A, Stella L, Lorenz C D and Kantorovich L 2016 PHYSICAL REVIEW B 93 ISSN 2469-9950 [[31]] [133] Wang J S 2007 Phys. Rev. Lett. 99 160601 ISSN 0031-9007 [[32]] [134] Schmid A 1982 J. Stat. Phys. 609–626 [[32]] [135] Dammak H, Chalopin Y, Laroche M, Hayoun M and Greffet J J 2009 Phys. Rev. Lett. 103 190601 ISSN 0031-9007 [[32]] [136] Barrozo A H and De Koning M 2011 Phys. Rev. Lett. 107 198901 ISSN 00319007 [[32]] [137] Dammak H, Hayoun M, Chalopin Y and Greffet J J 2011 Phys. Rev. Lett. 107 198902 ISSN 0031-9007 [[32]] [138] Wang J S, Ni X and Jiang J W 2009 Phys. Rev. B 80 224302 [[32]] [139] Gunst T, Lu J T, Markussen T, Jauho A P and Brandbyge M 2012 2012 15th Int. Work. Comput. Electron. 1–5 [[32]] [140] Barrat J L and Rodney D 2011 J. Stat. Phys. 144 679–689 ISSN 0022-4715, 1572-9613 (Preprint 1103.6108) [[32]] [141] Ceriotti M, Bussi G and Parrinello M 2010 J. Chem. Theory Comput. 6 1170–1180 ISSN 1549-9618 [[32]] [142] Stella L, Lorenz C D and Kantorovich L 2014 PHYSICAL REVIEW B 89 ISSN 1098-0121 [[32]] [143] Savin A V, Kosevich Y A and Cantarero A 2012 Phys. Rev. B 86 064305 ISSN 1098-0121 [[32]] [144] Brenner D W, Shenderova O a, Harrison J a, Stuart S J, Ni B and Sinnott S B 2002 J. Phys. Condens. Matter 14 783–802 ISSN 0953-8984 (Preprint S0953-8984(02)31186-X) [[33]] [145] Lü J T, Gunst T, Hedegård P and Brandbyge M 2011 Beilstein J. Nanotechnol. 2 814–23 ISSN 2190-4286 [[33]] [146] Lü J T, Brandbyge M, Hedeg\aard P, Todorov T N, Dundas D, Hedegård P, Todorov T N and Dundas D 2012 Phys. Rev. B 85 245444 ISSN 1098-0121 [[32]] [147] Lü J T, Christensen R B, Wang J s, Hedegård P and Brandbyge M 2015 Phys. Rev. Lett. 114 096801 ISSN 0031-9007 [[32]] [148] Allen P B and Feldman J L 1993 Physical Review B 48 12581–12588 [[33]], [[34]] [149] Roche S and Mayou D 1997 Physical Review Letters 79 2548 [[34]] [150] Torres L, Roche S and Charlier J 2014 Introduction to Graphene-Based Nanomaterials: From Electronic Structure to Quantum Transport Introduction to Graphene-based Nanomaterials: From Electronic Structure to Quantum Transport (Cambridge University Press) ISBN 978- 1-107-03083-1 [[34]] [151] Lanczos C 1950 J.Res.Natl.Bur.Stand.B 45 255–282 [[35]] [152] Arfken G, Weber H and Harris F 2013 Mathematical Methods for Physicists: A Comprehensive Guide (Elsevier) [[37]] [153] Zuev Y M, Chang W and Kim P 2009 Physical Review Letters 102 096807 [[38]] [154] Isacsson A, Cummings A W, Colombo L, Colombo L, Kinaret J M and Roche S 2017 2D Materials 4 012002 [[38]], [[43]] [155] Balandin A A, Ghosh S, Bao W, Calizo I, Teweldebrhan D, Miao F and Lau C N 2008 Nano Letters 8 [[38]] [156] Charlier J C, Blase X and Roche S 2007 Reviews of Modern Physics 79 677–732 [[38]] [157] Cresti A, Nemec N, Biel B, Niebler G, Triozon F, Cuniberti G and Roche S 2008 Nano Research 1 361–394 [[38]], [[41]] [158] Barrios-Vargas J E, Mortazavi B, Cummings A W, Martinez-Gordillo R, Pruneda M, Colombo L, Rabczuk T and Roche S 2017 Nano Letters 17 1660–1664 [[38]], [[43]] [159] Savić I, Mingo N and Stewart D A 2008 Physical Review Letters 101 165502 [[38]], [[39]], [[41]] [160] Mingo N, Esfarjani K, Broido D A and Stewart D A 2010 Physical Review B 81 [[38]] [161] Ci L, Song L, Jin C, Jariwala D, Wu D, Li Y, Srivastava A, Wang Z F, Storr K, Balicas L, Liu F and Ajayan P M 2010 Nature Materials 9 430–435 [[38]] [162] Rodriguez-Nieva J F, Saito R, Costa S D and Dresselhaus M S 2012 Physical Review B 85 245406 [[39]] GF, QCL and KG methods in quantum thermal transport 48 [163] Wirtz L, Rubio A, de la Concha R A and Loiseau A 2003 Physical Review B 68 045425 [[42]] [164] Levendorf M P, Kim C J, Brown L, Huang P Y, Havener R W, Muller D A and Park J 2012 Nature 488 627–632 [[43]] [165] Liu L, Park J, Siegel D A, McCarty K F, Clark K W, Deng W, Basile L, Idrobo J C, Li A P and Gu G 2014 Science 343 163–167 [[43]] [166] Van Tuan D, Kotakoski J, Louvet T, Ortmann F, Meyer J C and Roche S 2013 Nano Letters 13 1730–1735 [[43]] [167] Cummings A W, Duong D L, Nguyen V L, Van Tuan D, Kotakoski J, Barrios Vargas J E, Lee Y H and Roche S 2014 Advanced Materials 26 5079–5094 [[43]]