Yangyu Guo
email: yyguo@iis.u-tokyo.ac.jp

Marc Bescond
email: bescond@iis.u-tokyo.ac.jp

Zhongwei Zhang

Mathieu Luisier

Masahiro Nomura

Sebastian Volz
email: volz@iis.u-tokyo.ac.jp

Quantum mechanical modeling of anharmonic phonon-phonon scattering in nanostructures

The coherent quantum effect has become increasingly important in the heat dissipation bottleneck of semiconductor nanoelectronics with the characteristic size recently shrinking down to a few nanometers scale. However, the quantum mechanical model remains elusive for anharmonic phonon-phonon scattering in extremely small nanostructures with broken translational symmetry. It is a long-term challenging task to correctly simulate quantum heat transport including anharmonic scattering at a scale relevant to practical applications. In this article, we present a clarified theoretical formulation of anharmonic phonon nonequilibrium Green's function (NEGF) formalism for both one-and three-dimensional nanostructures, through a diagrammatic perturbation expansion and an introduction of Fourier's representation to both harmonic and anharmonic terms. A parallelized computational framework with first-principle force constants input is developed for large-scale quantum heat transport simulation. Some crucial approximations in numerical implementation are investigated to ensure the balance between numerical accuracy and efficiency. A quantitative validation is demonstrated for the anharmonic phonon NEGF formalism and computational framework by modeling cross-plane heat transport through a silicon thin film. The phonon-phonon scattering is shown to appreciably enhance the thermal resistance or conductance in extremely small homogeneous or heterogeneous thin films. The present methodology provides a robust platform for the device quantum thermal modeling, as well as a study on the transition from coherent to incoherent heat transport in nanophononic crystals. This work thus paves the way to understand and to manipulate heat conduction via the wave nature of phonons.

I. INTRODUCTION

The manipulation and control of heat transport in a dielectric crystal mediated by phonons has been a significant issue in modern technology and industrial applications. The aim of low phonon thermal conductivity is pursued in thermoelectric [START_REF] Minnich | Bulk nanostructured thermoelectric materials: Current research and future prospects[END_REF] and thermal barrier coating materials [START_REF] Mukhopadhyay | Two-channel model for ultralow thermal conductivity of crystalline Tl 3 VSe 4[END_REF], whereas high phonon thermal conductivity is desired for heat dissipation problems in nanoelectronics [START_REF] Moore | Emerging challenges and materials for thermal management of electronics[END_REF]. In the past several decades, fruitful efforts have been made in controlling the heat transport by tailoring the mean free path (MFP) of phonons through nanostructures [START_REF] Cahill | Nanoscale thermal transport[END_REF][START_REF] Cahill | Nanoscale thermal transport. II. 2003-2012[END_REF][START_REF] Chen | Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons[END_REF][START_REF] Zhang | Nano/Microscale Heat Transfer[END_REF]. These efforts are mainly based on the particle nature of phonons [START_REF] Ziman | Electrons and Phonons: The Theory of Transport Phenomena in Solids[END_REF][START_REF] Guo | Phonon hydrodynamics and its applications in nanoscale heat transport[END_REF], since the characteristic size of the system is usually much larger than the dominant width (or so-called coherence length) of phonon wave packets. With the development of nanofabrication and manufacturing technology in recent years, heat conduction tuning via the wave nature of phonons based on nanophononic crystals (PnCs) also becomes possible [START_REF] Maldovan | Sound and heat revolutions in phononics[END_REF][START_REF] Maldovan | Phonon wave interference and thermal bandgap materials[END_REF][START_REF] Volz | Nanophononics: State of the art and perspectives[END_REF][START_REF] Nomura | Thermal phonon engineering by tailored nanostructures[END_REF]. Distinctive behaviors of phonon heat transport have been experimentally demonstrated in this wave regime: the period-controlled lattice thermal conductivity minimum of superlattices (SLs) [START_REF] Ravichandran | Crossover from incoherent to coherent phonon scattering in epitaxial oxide superlattices[END_REF], the order of nano-PnC slowing down the heat propagation [START_REF] Maire | Heat conduction tuning by wave nature of phonons[END_REF], nanodot induced Anderson localization of thermal phonons in SLs [START_REF] Luckyanova | Phonon localization in heat conduction[END_REF], the aperiodicity of SLs reducing thermal conductivity via multiple localization [START_REF] Hu | Machine-Learning-Optimized Aperiodic Superlattice Minimizes Coherent Phonon Heat Conduction[END_REF], and so on [START_REF] Yu | Reduction of thermal conductivity in phononic nanomesh structures[END_REF][START_REF] Luckyanova | Coherent phonon heat conduction in superlattices[END_REF][START_REF] Zen | Engineering thermal conductance using a two-dimensional phononic crystal[END_REF].

The modeling of phonon heat transport in this wavelike coherent regime to particlelike incoherent regime remains, however, a challenging task. The widely adopted semiclassical phonon Boltzmann equation has become no longer valid [START_REF] Chen | Thermal conductivity and ballistic-phonon transport in the cross-plane direction of superlattices[END_REF][START_REF] Ravichandran | Coherent and incoherent thermal transport in nanomeshes[END_REF] as it describes only the evolution of phonon-particle population whereas it ignores the crucial phase information (interference effect) of phonon modes in this situation. There has been some theoretical effort to describe such coherent heat transport by a lattice dynamical model [START_REF] Maire | Heat conduction tuning by wave nature of phonons[END_REF][START_REF] Simkin | Minimum Thermal Conductivity of Superlattices[END_REF][START_REF] Yang | Partially coherent phonon heat conduction in superlattices[END_REF], with the anharmonic phonon scattering included by an imaginary wave vector component dependent on MFP in a phenomenological way. The molecular dynamics (MD) simulation is another popular approach to study the transition from coherent to incoherent heat transport in nano-PnC such as SLs [START_REF] Chen | Minimum superlattice thermal conductivity from molecular dynamics[END_REF][START_REF] Wang | Decomposition of coherent and incoherent phonon conduction in superlattices and random multilayers[END_REF][START_REF] Latour | Microscopic description of thermal-phonon coherence: From coherent transport to diffuse interface scattering in superlattices[END_REF][START_REF] Zhu | Phonon transport on two-dimensional graphene/boron nitride superlattices[END_REF]. As MD is a classical method, some essential quantum behavior of phonons would be lost at relatively low temperatures or in extremely small nanostructures. On the other hand, the empirical atomic interaction potential employed in MD might not be sufficiently accurate or even not available sometimes. Therefore, a full quantum mechanical model without any empirical input parameters is highly desired, which is the main focus of this work.

The nonequilibrium Green's function (NEGF) method [START_REF] Wang | Quantum thermal transport in nanostructures[END_REF][START_REF] Wang | Nonequilibrium Green's function method for quantum thermal transport[END_REF][START_REF] Chen | Thermal engineering in lowdimensional quantum devices: A tutorial review of nonequilibrium Green's function methods[END_REF][START_REF] Sadasivam | The atomistic Green's function method for interfacial phonon transport[END_REF] represents such a full quantum mechanical approach which could principally account for the coherent interaction and the incoherent scattering in the same footprint. In comparison to the Boltzmann transport theory, the NEGF method describes the evolution of both population and coherences of phonon modes [START_REF] Cohen-Tannoudji | Quantum Mechanics[END_REF] such that it will be able to capture the wavelike to particlelike behaviors of phonon heat transport. The NEGF formalism was originally proposed around the middle of the last century [START_REF] Martin | Theory of many-particle systems. I[END_REF][START_REF] Kadanoff | Quantum Statistical Mechanics[END_REF][START_REF] Keldysh | Diagram technique for nonequilibrium processes[END_REF] during the development of quantum field theory [START_REF] Rammer | Quantum Field Theory of Non-equilibrium States[END_REF] mainly for electrons, and has been relatively well established and widely applied in quantum transport modeling of nanoelectronic devices [START_REF] Datta | Nanoscale device modeling: The Green's function method[END_REF][START_REF] Anantram | Modeling of nanoscale devices[END_REF][START_REF] Haug | Quantum Kinetics in Transport and Optics of Semiconductors[END_REF][START_REF] Luisier | Atomistic full-band simulations of silicon nanowire transistors: Effects of electron-phonon scattering[END_REF]. Later the NEGF method was introduced to model phonon transport in the early years of this century attributed to several pioneering works [START_REF] Mingo | Phonon transport in nanowires coated with an amorphous material: An atomistic Green's function approach[END_REF][START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Wang | Nonequilibrium Green's function approach to mesoscopic thermal transport[END_REF][START_REF] Yamamoto | Nonequilibrium Green's Function Approach to Phonon Transport in Defective Carbon Nanotubes[END_REF]. Due to the large computational cost, the phonon NEGF has been often applied for ballistic heat transport through relatively simple structures like lowdimensional nanostructures (ultrathin nanowire [START_REF] Mingo | Phonon transport in nanowires coated with an amorphous material: An atomistic Green's function approach[END_REF][START_REF] Markussen | Electron and phonon transport in silicon nanowires: Atomistic approach to thermoelectric properties[END_REF][START_REF] Zhang | Simulation of phonon transport across a non-polar nanowire junction using an atomistic Green's function method[END_REF], nanotubes [START_REF] Yamamoto | Nonequilibrium Green's Function Approach to Phonon Transport in Defective Carbon Nanotubes[END_REF][START_REF] Savić | Phonon Transport in Isotope-Disordered Carbon and Boron-Nitride Nanotubes: Is Localization Observable?[END_REF], molecular junctions [START_REF] Galperin | Heat conduction in molecular transport junctions[END_REF], 2D structures [START_REF] Jiang | Minimum thermal conductance in graphene and boron nitride superlattice[END_REF][START_REF] Ouyang | Thermal transport of isotopic-superlattice graphene nanoribbons with zigzag edge[END_REF][START_REF] Zhang | A theoretical prediction of super high-performance thermoelectric materials based on MoS 2 /WS 2 hybrid nanoribbons[END_REF], etc.), interfaces [START_REF] Zhang | Simulation of interfacial phonon transport in Si-Ge heterostructures using an atomistic Green's function method[END_REF][START_REF] Tian | Enhancing phonon transmission across a Si/Ge interface by atomic roughness: First-principles study with the Green's function method[END_REF][START_REF] Li | Effect of lattice mismatch on phonon transmission and interface thermal conductance across dissimilar material interfaces[END_REF][START_REF] Ong | Efficient approach for modeling phonon transmission probability in nanoscale interfacial thermal transport[END_REF][START_REF] Latour | Ab initio study of mode-resolved phonon transmission at Si/Ge interfaces using atomistic Green's functions[END_REF][START_REF] Polanco | Phonon thermal conductance across GaN-AlN interfaces from first principles[END_REF], and SLs [START_REF] Hu | Machine-Learning-Optimized Aperiodic Superlattice Minimizes Coherent Phonon Heat Conduction[END_REF][START_REF] Tian | Green's function studies of phonon transport across Si/Ge superlattices[END_REF][START_REF] Mendoza | Anderson localization of thermal phonons leads to a thermal conductivity maximum[END_REF][START_REF] Cheaito | Interplay between total thickness and period thickness in the phonon thermal conductivity of superlattices from the nanoscale to the microscale: Coherent versus incoherent phonon transport[END_REF]. Very few works [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Wang | Nonequilibrium Green's function approach to mesoscopic thermal transport[END_REF][START_REF] Wang | Nonequilibrium Green's function method for thermal transport in junctions[END_REF][START_REF] Li | Anharmonic phonon transport in atomic wire coupled by thermal contacts with surface bond reconstruction[END_REF][START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF][START_REF] Lee | Anharmonic phonon-phonon scattering modeling of three-dimensional atomistic transport: An efficient quantum treatment[END_REF][START_REF] Dai | Rigorous formalism of anharmonic atomistic Green's function for three-dimensional interfaces[END_REF] directly take into account the anharmonic phonon-phonon scattering, yet often consider few-atom systems like the atomic chain or junction [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Wang | Nonequilibrium Green's function approach to mesoscopic thermal transport[END_REF][START_REF] Wang | Nonequilibrium Green's function method for thermal transport in junctions[END_REF][START_REF] Li | Anharmonic phonon transport in atomic wire coupled by thermal contacts with surface bond reconstruction[END_REF], and single-unit-cell interface [START_REF] Dai | Rigorous formalism of anharmonic atomistic Green's function for three-dimensional interfaces[END_REF]. An indirect treatment of incoherent phonon scattering was also proposed by the Büttiker probe approach [START_REF] Miao | Buttiker probes for dissipative phonon quantum transport in semiconductor nanostructures[END_REF][START_REF] Sadasivam | Thermal transport across metal silicide-silicon interfaces: First-principles calculations and Green's function transport simulations[END_REF][START_REF] Chu | Thermal boundary resistance predictions with non-equilibrium Green's function and molecular dynamics simulations[END_REF] inspired by its counterpart in electron NEGF, yet it relies on fitting the anharmonic scattering rates with empirical expressions. The challenge in NEGF modeling of anharmonic heat transport comes from not only the difficulty in numerical implementation, but also the less established theoretical formalism. As will be declared later in Sec. II, the crucial self-energy expressions for phonon-phonon scattering remain still diverse in several prototypical works in the literature [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Wang | Nonequilibrium Green's function approach to mesoscopic thermal transport[END_REF][START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF]. As a result, more work still needs to be done to form a solid quantitative validation of the anharmonic phonon NEGF formalism. Another important issue in phonon NEGF modeling is the input of the force constant (FC) matrix. The empirical atomic interaction potential is usually adopted to extract the harmonic FC matrix [START_REF] Mingo | Phonon transport in nanowires coated with an amorphous material: An atomistic Green's function approach[END_REF][START_REF] Yamamoto | Nonequilibrium Green's Function Approach to Phonon Transport in Defective Carbon Nanotubes[END_REF][START_REF] Markussen | Electron and phonon transport in silicon nanowires: Atomistic approach to thermoelectric properties[END_REF][START_REF] Zhang | Simulation of phonon transport across a non-polar nanowire junction using an atomistic Green's function method[END_REF][START_REF] Galperin | Heat conduction in molecular transport junctions[END_REF][START_REF] Jiang | Minimum thermal conductance in graphene and boron nitride superlattice[END_REF][START_REF] Ouyang | Thermal transport of isotopic-superlattice graphene nanoribbons with zigzag edge[END_REF][START_REF] Zhang | A theoretical prediction of super high-performance thermoelectric materials based on MoS 2 /WS 2 hybrid nanoribbons[END_REF][START_REF] Zhang | Simulation of interfacial phonon transport in Si-Ge heterostructures using an atomistic Green's function method[END_REF][START_REF] Li | Effect of lattice mismatch on phonon transmission and interface thermal conductance across dissimilar material interfaces[END_REF][START_REF] Ong | Efficient approach for modeling phonon transmission probability in nanoscale interfacial thermal transport[END_REF] and anharmonic FC matrix [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Wang | Nonequilibrium Green's function approach to mesoscopic thermal transport[END_REF][START_REF] Wang | Nonequilibrium Green's function method for thermal transport in junctions[END_REF][START_REF] Li | Anharmonic phonon transport in atomic wire coupled by thermal contacts with surface bond reconstruction[END_REF][START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF][START_REF] Lee | Anharmonic phonon-phonon scattering modeling of three-dimensional atomistic transport: An efficient quantum treatment[END_REF]. In recent years, attributed to the advance in density-functional theory (DFT) [START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF] and computational power, the first-principle input has been also introduced into phonon NEGF codes, yet mostly harmonic FC matrix in ballistic heat transport [START_REF] Savić | Phonon Transport in Isotope-Disordered Carbon and Boron-Nitride Nanotubes: Is Localization Observable?[END_REF][START_REF] Tian | Enhancing phonon transmission across a Si/Ge interface by atomic roughness: First-principles study with the Green's function method[END_REF][START_REF] Latour | Ab initio study of mode-resolved phonon transmission at Si/Ge interfaces using atomistic Green's functions[END_REF][START_REF] Polanco | Phonon thermal conductance across GaN-AlN interfaces from first principles[END_REF][START_REF] Tian | Green's function studies of phonon transport across Si/Ge superlattices[END_REF][START_REF] Mendoza | Anderson localization of thermal phonons leads to a thermal conductivity maximum[END_REF][START_REF] Cheaito | Interplay between total thickness and period thickness in the phonon thermal conductivity of superlattices from the nanoscale to the microscale: Coherent versus incoherent phonon transport[END_REF], except for a recent work [START_REF] Dai | Rigorous formalism of anharmonic atomistic Green's function for three-dimensional interfaces[END_REF] also considering the anharmonic one.

The aim of the present work includes several aspects. Firstly, we aim to present a clarified formulation of the anharmonic phonon NEGF method through a diagrammatic perturbation expansion of Green's function and a thorough comparison to existing results. As a further step, we will extend the formalism to anharmonic heat transport through 3D nanostructures with transverse periodicity by introducing a Fourier's representation to both the harmonic and anharmonic terms. Secondly, we aim to develop a firstprinciple-based numerical framework for large-scale quantum heat transport simulation by introducing advanced computation techniques as well as DFT input of harmonic and anharmonic FCs. Finally, a quantitative validation of the anharmonic phonon NEGF formalism and numerical framework will be demonstrated by several benchmark studies. The crucial approximations beyond the treatment in a previous work [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF] by one of the coauthors will be elucidated to ensure feasible yet still accurate large-scale simulation. The remainder of this article will be arranged as follows: The theoretical formulations and numerical implementation of the phonon NEGF method will be described in Sec. II; the results and related discussions will be shown in Sec. III; concluding remarks are made in Sec. IV.

II. MATHEMATICAL AND NUMERICAL MODELS

In this section, we will present a clarified anharmonic phonon NEGF formalism for each relevant physical system. Section II A will firstly focus on 1D nanostructures without periodicity. Then a Fourier's representation is introduced to the anharmonic formalism in Sec. II B for 3D nanostructures with transverse periodicity, which is a paradigm for quantum heat transport across thin films, interfaces, SLs, and multilayer nanostructures. In Sec. II C, the numerical implementation of the phonon NEGF formalism will be elucidated in detail, including the computational strategies and techniques for large-scale simulation. In Sec. II D, the input of harmonic and anharmonic FC matrices from first-principle (DFT) calculation will be introduced.

A. Anharmonic phonon NEGF formalism for 1D nanostructures

In this subsection, we firstly consider quantum heat transport through a 1D nanostructure without any translational periodicity, as shown in Fig. 1(a). The governing equation for the retarded Green's function of the device in a matrix notation is [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF]]

G R (ω) = [ω 2 I --R (ω)] -1 , ( 1
)
where I denotes the unity matrix, denotes the second-order FC matrix [cf. Eq. (A4) in Appendix A], and the superscript "-1" means the inverse of a matrix. Note that throughout this work the bold character represents a matrix, whereas its component is denoted by a regular character. ω represents the frequency of the phonon modes through the system. The retarded self-energy matrix includes the contribution from both of the two contacts and the anharmonic phonon-phonon scattering within the device:

R (ω) = R 1 (ω) + R 2 (ω) + R s (ω). (2)
The retarded contact self-energy matrix R 1(2) (ω) in Eq. ( 2) is related to the surface Green's function of the contacts, which are calculated by the decimation technique [START_REF] Sadasivam | The atomistic Green's function method for interfacial phonon transport[END_REF][START_REF] Guinea | Effective twodimensional Hamiltonian at surfaces[END_REF]. For brevity, the subscript "1(2)" is introduced to refer to contact 1 or contact 2 hereafter. The calculation of the retarded scattering self-energy R s (ω) will be introduced later. For ballistic heat transport problems, usually we only need to resolve the retarded Green's function of the device in Eq. (1) to get the transmission through the system using the Caroli formula [START_REF] Wang | Quantum thermal transport in nanostructures[END_REF][START_REF] Mingo | Phonon transport in nanowires coated with an amorphous material: An atomistic Green's function approach[END_REF][START_REF] Yamamoto | Nonequilibrium Green's Function Approach to Phonon Transport in Defective Carbon Nanotubes[END_REF]. However, for anharmonic heat transport, we also have to compute the greater/lesser Green's function as [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF]

G >, < (ω) = G R (ω) >, < (ω)G A (ω), (3)
where the advanced Green's function of the device is the Hermitian conjugate of the retarded Green's function: [G R (ω)] † . The greater/lesser self-energy matrix also includes the contribution from both of the two contacts and the anharmonic phonon-phonon scattering in the device:

G A (ω) =
>,< (ω) = >,< 1 (ω) + >,< 2 (ω) + >,< s (ω). ( 4
)
The greater/lesser contact self-energy matrix in Eq. ( 4) is related to the retarded contact self-energy matrix as

> 1(2) (ω) = [1 + f eq (ω, T 1(2) )] R 1(2) (ω) -R 1(2) (ω) † , ( 5
) < 1(2) (ω) = f eq (ω, T 1(2) ) R 1(2) (ω) -R 1(2) (ω) † , (6)
where f eq (ω, T 1(2) ) denotes the Bose-Einstein equilibrium phonon distribution at the contact temperature T 1 [START_REF] Mukhopadhyay | Two-channel model for ultralow thermal conductivity of crystalline Tl 3 VSe 4[END_REF] . Equations ( 5) and ( 6) are crucial relations to introduce the isothermal boundary conditions to the steady-state heat transport in anharmonic phonon NEGF simulation.

A crucial part of the anharmonic phonon NEGF formalism is the calculation of the greater/lesser scattering self-energy matrix >,< s (ω) in Eq. ( 4). The retarded scattering self-energy matrix in Eq. ( 2) is then related to the greater/lesser scattering self-energy matrix as [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF] R

s (ω) = 1 2 [ > s (ω) -< s (ω)] + iP ∞ -∞ dω 2π > s (ω ) -< s (ω ) ω -ω . ( 7
)
The second term on the right-hand side of Eq. ( 7) is the Cauchy principal integral, and is often neglected for computational simplicity [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF]. This term represents the frequency shift of phonon modes due to the anharmonic phonon-phonon scattering and is usually very small for most applications. In terms of the greater/lesser scattering self-energy, we find similar expressions as available in the literature yet with different coefficients [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Wang | Nonequilibrium Green's function approach to mesoscopic thermal transport[END_REF][START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF], as summarized in Table I. For a clear comparison, we keep the notation of third-order FC in the original literature whereas we rewrite the notation of other variables (atomic displacement, Green's function, and self-energy) to be consistent with that in this work. This discrepancy might be due to the very complicated derivation process of the anharmonic scattering self-energy through the diagrammatic perturbation expansion of phonon Green's function, the details of which are often not provided [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Wang | Nonequilibrium Green's function approach to mesoscopic thermal transport[END_REF][START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF]. Since there is no concretely validated and widely accepted expression of anharmonic scattering self-energy for heat transport in nanostructures, our strategy to resolve this issue includes a twofold procedure: (1) Firstly we conduct a diagrammatic perturbation derivation of the self-energy expression by ourselves; (2) then we demonstrate a quantitative validation of the theoretical formalism.

Through the diagrammatic perturbation expansion, we obtain the expression of the greater/lesser self-energy matrix for third-order anharmonic phonon scattering as

>, <i j s, ll (ω) = 1 2 i h l 1 l 2 l 3 l 4 j 1 j 2 j 3 j 4 ∞ -∞ dω 2π i j 1 j 2 ll 1 l 2 j j 3 j 4 l l 3 l 4 × G >, < j 1 j 4 l 1 l 4 (ω )G >, < j 2 j 3 l 2 l 3 (ω -ω ), (8)
where i j 1 j 2 ll 1 l 2 is the component of the third-order FC matrix [cf. Eq. (A5) in Appendix A], whereas the subscripts l (l 1 , l 2 , • • • ) denote the atomic index, and the superscripts i, j denote the Cartestian coordinate index (x, y, z). The detailed derivation TABLE I. Summary of expressions for greater/lesser anharmonic phonon scattering self-energy.

Third-order term in Hamiltonian

Scattering self-energy References

i jk V (3) i jk u i u j u k i < s,in (ω) = h jklm ∞ -∞ V (3) i jk G < jl (ω )G < km (ω -ω )V (3) lmn dω Equation (27) in Ref. [43] 1 3 i jk T i jk u i u j u k < s, jk (ω) = 2i h lmrs T jlm T rsk ∞ -∞ G < lr (ω )G < ms (ω -ω ) dω 2π
First term of Eq. ( 13) in Ref. [START_REF] Wang | Nonequilibrium Green's function approach to mesoscopic thermal transport[END_REF]

1 3! i jk lmn dV (3)i jk lmn u i l u j m u k n >, <i j s, nm (ω) = 2i h k 1 k 2 k 3 k 4 l 1 l 2 l 3 l 4 ∞ -∞ dω 2π dV (3)ik 1 k 2 nl 1 l 2 dV (3)k 3 k 4 j l 3 l 4 m ×G >, <k 1 k 3 l 1 l 3 (ω + ω )G <, >k 4 k 2 l 4 l 2 (ω ) Equation (7) in Ref. [64] 1 3! nml i jk i jk nml u i n u j m u k l >, <i j s, ll (ω) = 1 2 i h l 1 l 2 l 3 l 4 j 1 j 2 j 3 j 4 ∞ -∞ dω 2π i j 1 j 2 ll 1 l 2 j j 3 j 4 l l 3 l 4 × G >, < j 1 j 4 l 1 l 4 (ω )G >, < j 2 j 3 l 2 l 3 (ω -ω )
Present process of Eq. ( 8) is given in Appendix A. In comparison to the expressions of anharmonic scattering self-energy in Refs. [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF], there exists a factor difference of π /9 and 4, respectively, with respect to the present result of Eq. [START_REF] Ziman | Electrons and Phonons: The Theory of Transport Phenomena in Solids[END_REF].

Considering that the third-order FC (T i jk in Table I) is half of the present one ( i jk nml ), the part of the self-energy corresponding to the three-phonon anharmonic scattering in Ref. [START_REF] Wang | Nonequilibrium Green's function approach to mesoscopic thermal transport[END_REF] is actually consistent with the present expression. A further theoretical corroboration of Eq. ( 8) will be shown in Sec. II B when we compare it to the expression of anharmonic phonon scattering self-energy in bulk materials. Up to now, Eqs. (1), (3), [START_REF] Zhang | Nano/Microscale Heat Transfer[END_REF], and (8) constitute a closed set of coupled matrix equations, which have to be solved through an iterative process called the self-consistent Born approximation (SCBA) [START_REF] Luisier | Atomistic full-band simulations of silicon nanowire transistors: Effects of electron-phonon scattering[END_REF][START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF]. Note Eqs. ( 1) and ( 3) are the steady-state limit of the equations of motion for real-time Green's functions, which have to be solved for transient problems such as ultrafast heat transport [START_REF] Guo | Phonon hydrodynamics for nanoscale heat transport at ordinary temperatures[END_REF][START_REF] Hua | Space-time dependent thermal conductivity in nonlocal thermal transport[END_REF]. It is computationally very challenging to model transient quantum heat transport by phonon NEGF formalism due to the memory effect in temporal integration of Green's function and self-energy. There has been much progress in electron NEGF for modeling transient phenomena, such as the generalized Kadanoff-Baym ansatz framework [START_REF] Hirsbrunner | A review of modeling interacting transient phenomena with non-equilibrium Green functions[END_REF]. It is a good strategy to borrow the techniques from the transient electron NEGF, which are well summarized in a very recent review [START_REF] Hirsbrunner | A review of modeling interacting transient phenomena with non-equilibrium Green functions[END_REF].

B. Anharmonic phonon NEGF formalism for 3D nanostructures

In this subsection, we consider quantum heat transport through 3D nanostructures with transverse periodicity as shown in Fig. 1(b). Attributed to this periodicity in the transverse direction, the phonon NEGF formalism in Sec. II A can be rewritten into a Fourier's representation by introducing a transverse wave vector q ⊥ and the following Fourier transform for the Green's function:

G i j ll (ω) = 1 N q ⊥ exp (iq ⊥ • R ⊥ )G i j l x l x (ω; q ⊥ ), (9)
G i j l x l x (ω; q ⊥ ) = R ⊥ exp (-iq ⊥ • R ⊥ )G i j ll (ω), ( 10
)
where N is the number of transverse wave vectors, and R ⊥ = (l yl y )a 2 + (l zl z )a 3 with a 2 and a 3 the lattice vectors along the y and z directions, respectively. Here l ≡ (l x , l y , l z ) denotes the index of the lattice unit cell in the device along the x, y, and z directions. Note that in the Green's function (and related variables such as self-energy matrix and FCs matrix), the subscript l includes the index of atoms within the corresponding lattice unit cell, i.e., G i j ll (ω) ≡ G i j lκ, l κ (ω) with κ the atomic index in a lattice unit cell. For the sake of clarity, we only keep the index of the lattice unit cell.

With the help of Eqs. ( 9) and ( 10), the formalism in Sec. II A can be rewritten into the Fourier's representation, and the main governing equations Eqs. ( 1), ( 3), [START_REF] Zhang | Nano/Microscale Heat Transfer[END_REF], and ( 8) become

G R (ω; q ⊥ ) = [ω 2 I -˜ (q ⊥ ) -R (ω; q ⊥ )] -1 , (11)
G >,< (ω; q ⊥ ) = G R (ω; q ⊥ ) >,< (ω; q ⊥ )G A (ω; q ⊥ ), (12) R s
(ω; q ⊥ ) = 1 2 [ > s (ω; q ⊥ ) -< s (ω; q ⊥ )] + iP ∞ -∞ dω 2π > s (ω ; q ⊥ ) -< s (ω ; q ⊥ ) ω -ω , ( 13
)
<, >i j s, l x l x (ω; q ⊥ ) =

1 2 i h l 1x l 2x l 3x l 4x j 1 j 2 j 3 j 4 1 N q ⊥ ˜ i j 1 j 2 l x l 1x l 2x (q ⊥ , q ⊥ -q ⊥ ) ˜ j j 3 j 4 l x l 3x l 4x (q ⊥ -q ⊥ , -q ⊥ ) × ∞ -∞ dω 2π G <, > j 1 j 4 l 1x l 4x (ω ; q ⊥ )G <, > j 2 j 3 l 2x l 3x (ω -ω ; q ⊥ -q ⊥ ). (14)
In Eq. ( 11), the Fourier's representation of the harmonic FC matrix is defined as

˜ i j l x l x (q ⊥ ) = R ⊥ i j ll exp (-iq ⊥ • R ⊥ ), (15)
which is consistent with the definition in the previous ballistic NEGF formalism [START_REF] Zhang | Simulation of interfacial phonon transport in Si-Ge heterostructures using an atomistic Green's function method[END_REF][START_REF] Polanco | Phonon thermal conductance across GaN-AlN interfaces from first principles[END_REF]. In Eq. ( 14), the Fourier's representation of the third-order anharmonic FC matrix is defined as

˜ i j 1 j 2 l x l 1x l 2x (q ⊥ , q ⊥ ) = R ⊥ R ⊥ i j 1 j 2 ll 1 l 2 exp (-iq ⊥ • R ⊥ ) × exp (-iq ⊥ • R ⊥ ), (16)
with R ⊥ = (l yl 1y )a 2 + (l zl 1z )a 3 , R ⊥ = (l yl 2y )a 2 + (l zl 2z )a 3 here. The third-order FC matrix has a dependence on two transverse wave vectors since the three-body interaction depends on two relative displacements ( R ⊥ and R ⊥ ) in the cross-sectional direction. Due to the translational invariance along the transverse direction, the right-hand summation term in Eq. ( 16) will no longer be dependent on the (l y , l z ) pair. The basic idea for the derivation of Eq. ( 14) is to put Eq. ( 9) and the inverse transform of Eq. ( 16) into Eq. ( 8). The detailed derivation process is a bit cumbersome and is shown in Appendix B.

In a recent contribution to the anharmonic phonon NEGF formalism for 3D interfaces [START_REF] Dai | Rigorous formalism of anharmonic atomistic Green's function for three-dimensional interfaces[END_REF], a different third-order tensor Fourier's decomposition using the P matrix is developed for the anharmonic FC matrix. The authors obtain the anharmonic phonon scattering self-energy matrix below [START_REF] Dai | Rigorous formalism of anharmonic atomistic Green's function for three-dimensional interfaces[END_REF]:

<, > s, ur (ω; q ⊥ ) = i h vwpq q ⊥ Ṽuvw (q ⊥ , q ⊥ ) Ṽpqr (q ⊥ , q ⊥ ) ∞ -∞ dω × G <, > qv (ω ; q ⊥ )G <, > w p (ω -ω ; q ⊥ ), (17)
where we have kept most of the notations therein. In spite of a similar mathematical form, our result of Eq. ( 14) is different from Eq. ( 17) in terms of the following three aspects:

(1) a factor difference since their development is based on the anharmonic phonon scattering self-energy in Ref. [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF], as is already shown in Table I; (2) both the energy conservation [ω = ω + (ω -ω )] and the momentum conservation [q ⊥ = q ⊥ + (q ⊥q ⊥ )] are automatically satisfied in our expression (as explicitly shown in the Feynman diagram in Fig. 2) whereas only the energy conservation is ensured in Eq. ( 17); (3) a normalization over the number of transverse wave vectors is also included before the sum over q ⊥ in our formulation.

The expression of anharmonic phonon-phonon scattering self-energy in Eq. ( 14) can be further corroborated by its counterpart in bulk material. In comparison to the present 3D nanostructures with transverse periodicity, the bulk material also has periodicity in the transport direction. Through a diagrammatic perturbation expansion of the phonon Green's function similar to that in Appendix A, one could obtain the self-energy for three-phonon anharmonic scattering in bulk material as follows [START_REF] Xu | Nonequilibrium Green's function method for phonon-phonon interactions and ballistic-diffusive thermal transport[END_REF]:

< s (ω; q) = 2i h q 1 q 2 ∞ -∞ dω 1 2π ∞ -∞ dω 2 2π 2πδ(ω -ω 1 -ω 2 )
× |F (-q, q 1 , q 2 )| 2 G (0)< (ω 1 ; q 1 )G (0)< (ω 2 ; q 2 ), [START_REF] Yu | Reduction of thermal conductivity in phononic nanomesh structures[END_REF] where q is the 3D wave vector in bulk material, |F (-q, q 1 , q 2 )| 2 = F (-q, q 1 , q 2 ) × F * (-q, q 1 , q 2 ) with F * (-q, q 1 , q 2 ) = F (q, -q 1 , -q 2 ). Here F (q 1 , q 2 , q 3 ) is the third-order anharmonic dynamic matrix contributing to the anharmonic term in the Hamiltonian as [START_REF] Xu | Nonequilibrium Green's function method for phonon-phonon interactions and ballistic-diffusive thermal transport[END_REF]

V = 1 3 q 1 q 2 q 3 F (q 1 , q 2 , q 3 )A q 1 A q 2 A q 3 , ( 19
)
where A q = h/2ω q (a q + a † -q ) is the phonon normal coordinate operator, with a q , a † q the phonon destruction and creation operator, respectively. Due to the translational symmetry of the 3D lattice, F (q 1 , q 2 , q 3 ) includes a delta function as (q 1 + q 2 + q 3 ) which will be nonvanishing only when q 1 + q 2 + q 3 = 0 or a reciprocal lattice vector [START_REF] Xu | Nonequilibrium Green's function method for phonon-phonon interactions and ballistic-diffusive thermal transport[END_REF][START_REF] Maradudin | Theory of Lattice Dynamics in the Harmonic Approximation[END_REF]. Thus Eq. ( 18) can be rewritten into

< s (ω; q) = 2i h q 1 ∞ -∞ dω 1 2π F (-q, q 1 , q -q 1 )
× F (q, -q 1 , q 1q)G (0)< (ω 1 ; q 1 )

× G (0)< (ω -ω 1 ; q -q 1 ). ( 20
)
Considering the different notations in the third-order term of the Hamiltonian between Eq. ( 19) and that [Eq. (A3) in Appendix A] in the present work (a factor 1/3 versus 1/3!), Eq. ( 20) will be exactly consistent with Eq. ( 14) once the TABLE II. General relations between the nonequilibrium phonon Green's functions for 1D and 3D nanostructures.

1D nanostructures [START_REF] Wang | Nonequilibrium Green's function method for thermal transport in junctions[END_REF] 3D nanostructures (present result)

[G > (ω)] T = G < (-ω) [ G > (ω; q ⊥ )] T = G < (-ω; -q ⊥ ) [G < (ω)] † = -G < (ω) [ G < (ω; q ⊥ )] † = -G < (ω; q ⊥ ) [G R (ω)] * = G R (-ω) [ G R (ω; q ⊥ )] * = G R (-ω; -q ⊥ ) [G R (ω)] † = G A (ω) [ G R (ω; q ⊥ )] † = G A (ω; q ⊥ )
periodicity of the transport direction is released. Note that the normalization factor 1/N is included in |F (-q, q 1 , q 2 )| 2 [START_REF] Xu | Nonequilibrium Green's function method for phonon-phonon interactions and ballistic-diffusive thermal transport[END_REF].

Finally, some general symmetry relations between the nonequilibrium phonon Green's functions for 3D nanostructures with transvese peridocity are derived:

[G > (ω; q ⊥ )] T = G < (-ω; -q ⊥ ), (21)
[G < (ω; q ⊥ )] † = -G < (ω; q ⊥ ), (22)
[G R (ω; q ⊥ )] * = G R (-ω; -q ⊥ ), (23)
[G R (ω; q ⊥ )] † = G A (ω; q ⊥ ). ( 24
)
The superscript "T" in Eq. ( 21) denotes the transpose of a matrix. We provide a detailed proof of Eq. ( 21) in Appendix C. The other relations can be proved in a similar way and are not shown here. A comparison to similar general relations between the phonon Green's function for 1D nanostructures [START_REF] Wang | Nonequilibrium Green's function method for thermal transport in junctions[END_REF] is summarized in Table II.

C. Numerical implementation

The numerical implementation of the anharmonic phonon NEGF formalism by the SCBA iterative solution of Eqs. ( 11)-( 14) is challenging for large-scale simulations because of the computational time cost and memory cost. The large time cost is mainly due to the intensive calculation of multiple summations and integrations in the anharmonic scattering self-energy in Eq. ( 14). Furthermore, the memory cost and time cost related to the matrix storage and operation will be proportional to (N x N d ) 2 and (N x N d ) 3 , respectively, if the full matrix of the device is directly resolved [START_REF] Anantram | Modeling of nanoscale devices[END_REF]. Here the device consists of N d slabs as shown in Fig. 3, with the matrix size for each slab being N x . Two advanced computational techniques are introduced in this subsection to retrieve the situation: (1) the recursive algorithm in Sec. II C 1; (2) the parallelization scheme in Sec. II C 2. We will also briefly introduce the macroscopic variable calculation in Sec. II C 3. FIG. 3. Schematic of recursive algorithm for numerical implementation of phonon NEGF. "LC" and "RC" denote left contact (contact 1) and right contact (contact 2) respectively, and the device is divided into N d slabs.

Recursive algorithm

The recursive algorithm was originally proposed in electron NEGF codes for large-scale nanoelectronic device modeling [START_REF] Svizhenko | Two-dimensional quantum mechanical modeling of nanotransistors[END_REF], and has been also applied in anharmonic phonon NEGF modeling of heat transport in nanowires [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF]. A similar variant of the recursive algorithm was also developed in ballistic phonon NEGF [START_REF] Li | Effect of lattice mismatch on phonon transmission and interface thermal conductance across dissimilar material interfaces[END_REF]. The motivation of the recursive algorithm arises from the fact that only the diagonal blocks and first off-diagonal blocks of Green's functions are required to compute the relevant macroscopic observables (i.e., local density and current) of the system. The recursive algorithm only deals with matrices for each slab and its adjacent slabs such that the memory cost and time cost will be linearly proportional to the system size along the transport direction as (N x ) 2 N d and (N x ) 3 N d , respectively [START_REF] Anantram | Modeling of nanoscale devices[END_REF]. In this work, we adopt the recursive algorithm [START_REF] Anantram | Modeling of nanoscale devices[END_REF][START_REF] Svizhenko | Two-dimensional quantum mechanical modeling of nanotransistors[END_REF] to compute the N x × N x diagonal blocks and first off-diagonal blocks of the retarded Green's function and the greater/lesser Green's functions. We also make other approximations in the numerical implementation as follows: (I) only the interaction within neighboring slabs in Fig. 3 is considered for harmonic FC; (II) only the atomic nearest-neighbor interaction is considered for third-order anharmonic FC; (III) only the N x × N x diagonal blocks of anharmonic scattering self-energy are considered, and within each diagonal block only the dominant terms between each atom and its nearest-neighbor atoms are computed; (IV) the anharmonic phonon-phonon scattering is considered only in the device region whereas the contacts are harmonic. The approximation (II) has also been adopted in a previous work [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF] based on the empirical valence-force-field model, and shall be more or less reasonable for weakly anharmonic materials like silicon considered in this work, as will be shown later in Sec. II D. In principle, more neighbors should be considered in the third-order FC for strongly anharmonic materials such as oxides, yet it remains computationally very challenging for large-scale simulations. Concerning approximation (III), an even simpler approximation is assumed in Ref. [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF] where only the 3 × 3 diagonal blocks of the anharmonic scattering self-energy matrix are computed. We will show in Sec. III C that the approximation in Ref. [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF] will significantly overestimate the anharmonic scattering rate and thus underestimate the thermal conductance. In terms of approximation (IV), it has been shown that the harmonic or anharmonic contact has a negligibly small influence on the thermal transport properties [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF].

Parallelization scheme

With a multiple-CPU (central processing unit) computer facility, the iterative solution of Eqs. ( 11)-( 14) can be parallelized based on the message-passing-interface (MPI) standard [START_REF] Gropp | A highperformance, portable implementation of the MPI message passing interface standard[END_REF]. The MPI parallelization scheme has been widely adopted in electron NEGF for large-scale device simulations [START_REF] Luisier | Atomistic full-band simulations of silicon nanowire transistors: Effects of electron-phonon scattering[END_REF][START_REF] Luisier | Proceedings of the 2010 ACM/IEEE International Conference for High Performance Computing, Networking, Storage and Analysis[END_REF], and is also used in anharmonic phonon NEGF in a previous work [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF] by one of the coauthors. In comparison to the parallelization scheme in Ref. [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF] for nanowire (1D nanostructures) simulations that require data exchange between different frequency points as shown in Eq. ( 8), here we also have to exchange the data between different transverse wave vectors for simulation of 3D nanostructures as shown in Eq. ( 14). Thus, we build a dual-level parallelization scheme for the transverse wave vector and frequency points, both of which are uniformly discretized as shown in Fig. 4. In the first level, all the transverse wave vectors are parallelly treated, i.e., one wave vector will be allocated to each CPU. In the second level, the frequency points are divided into several intervals which are parallelly treated, i.e., each CPU will receive a segment with several frequency points. For ballistic heat transport simulation, the problem is embarrassingly parallelized, i.e., the retarded Green's function and transmission of each mode (ω; q ⊥ ) is independently calculated in each CPU without the need of data exchange during computation. For anharmonic heat transport simulation, the situation is more complicated as the Green's function of one mode (ω; q ⊥ ) is coupled with many other modes as inferred from Eq. ( 14). Since the phonon Green's functions for each mode are often distributed in different CPUs, data exchange is needed when computing the anharmonic scattering self-energy. We design an algorithm for the data exchange, the details of which are provided in Appendix D.

Macroscopic variable calculation

The local heat flow from the qth slab to the (q + 1)th slab shown in Fig. 3 is related to the first off-diagonal blocks of the lesser phonon Green's function [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF] and can be computed as

J q→q+1 = - n∈q m∈q+1 i j 1 N q ⊥ ∞ 0 hω ˜ i j nm (q ⊥ ) × G <, ji mn (ω; q ⊥ ) -G <,i j nm (ω; q ⊥ ) ˜ ji mn (q ⊥ ) dω 2π (25)
Equation.( 25) can be written as

J q→q+1 = 1 N q ⊥ ∞ 0 J (ω; q ⊥ ) dω 2π
, with the heat flow contributed by each mode expressed as

J (ω; q ⊥ ) = - n∈q m∈q+1 i j hω ˜ i j nm (q ⊥ )G <, ji mn (ω; q ⊥ ) -G <,i j nm (ω; q ⊥ ) ˜ ji mn (q ⊥ ) . ( 26
)
The local heat flow distribution is used as a criterion of convergence of the SCBA scheme. After each iteration, the local heat flow is computed based on Eq. [START_REF] Chen | Minimum superlattice thermal conductivity from molecular dynamics[END_REF]. The SCBA convergence is reached when the conservation of heat flow is fullfilled along the transport direction, i.e., the local heat flow is the same for all the slabs within a certain allowed numerical error (1% in the present work).

For purely coherent heat transport, the definition of local temperature is usually not relevant. However, when considering anharmonic phonon-phonon scattering, the local temperature can be computed based on the local energy conservation condition within one slab:

∞ 0 ρ(ε, R n ) 2 hω h2 hω d ( hω) 2π = ∞ 0 LDOS(ε, R n ) f eq (ω, T eff ) 2 hω h2 hω d ( hω) 2π , ( 27
)
where R n denotes the spatial position of the qth slab, and ε ≡ ω 2 is the eigenvalue of the harmonic FC matrix. f eq (ω, T eff ) is the Bose-Einstein equilibrium phonon distribution at an effective local temperature T eff . The local phonon number in Eq. ( 27) is related to the diagonal blocks of the lesser phonon Green's function [START_REF] Rhyner | Atomistic modeling of coupled electron-phonon transport in nanowire transistors[END_REF] and is computed as

ρ(ε, R n ) = Tr 1 N q ⊥ iG < nn (ω; q ⊥ ) , (28)
with "Tr" denoting the trace of a square matrix, and the subscript "nn" represents the sub-block of G < corresponding to the qth slab. The local density of states (LDOS) in Eq. ( 27) is defined as LDOS(ε, R n ) = Tr[A nn (ω)], with the spectral function matrix computed as

A nn (ω) = 1 N q ⊥ A nn (ω; q ⊥ ) = 1 N q ⊥ i[G > nn (ω; q ⊥ ) -G < nn (ω; q ⊥ )]. ( 29
)
In a previous work on coupled electron-phonon NEGF modeling [START_REF] Rhyner | Atomistic modeling of coupled electron-phonon transport in nanowire transistors[END_REF], a different definition of local temperature was adopted based on the local phonon number conservation condition. This approach amounts to removing one hω in the integration on both sides of Eq. ( 27). Here we keep the convention of using local energy density to characterize the local temperature in nonequilibrium transport as currently done in the heat transport community [START_REF] Chen | Thermal conductivity and ballistic-phonon transport in the cross-plane direction of superlattices[END_REF][START_REF] Miao | Buttiker probes for dissipative phonon quantum transport in semiconductor nanostructures[END_REF][START_REF] Sadasivam | Thermal transport across metal silicide-silicon interfaces: First-principles calculations and Green's function transport simulations[END_REF][START_REF] Chu | Thermal boundary resistance predictions with non-equilibrium Green's function and molecular dynamics simulations[END_REF].

D. First-principle input

In this work, the material properties of silicon are adopted for the nanostructures and are obtained by first-principle (DFT) calculation. The DFT calculation is implemented in the open-source package QUANTUM ESPRESSO (QE) [START_REF] Giannozzi | QUANTUM ESPRESSO: A modular and opensource software project for quantum simulations of materials[END_REF] with norm-conserving pseudopotential and the LDA (local density approximation) exchange-correlation functional for silicon. A kinetic energy cutoff of 60 Ry is used for the wave function and a self-consistent convergence threshold of 10 -12 is adopted after independence check. Firstly, a relaxation process on a primitive unit cell is run to obtain an optimized lattice constant of 5.4018 Å with an electronic wave vector grid of 8 × 8 × 8. For the harmonic FC, the finite displacement method is used as implemented in the open-source package PHONOPY [START_REF] Togo | First principles phonon calculations in materials science[END_REF] combined with the DFT calculation in QE. A supercell of 3 × 3 × 3 conventional unit cells (216 atoms) is considered, with all the interactions within the supercell included. An electronic wave vector of 2 × 2 × 2 is used in DFT calculation. For the third-order anharmonic FC, the finite displacement method is also used as implemented in the open-source package THIRDORDER [START_REF] Li | SHENGBTE: A solver of the Boltzmann transport equation for phonons[END_REF] combined with the DFT calculation in QE. A supercell of 2 × 2 × 2 conventional unit cells (64 atoms) is adopted. In terms of the atomic interaction range, we consider, respectively, the first, second, and third nearest-neighbor shells to compare the results. An electronic wave vector of 1 × 1 × 1 is taken in DFT calculation. After extracting both the harmonic and third-order anharmonic FCs, we use them to calculate the bulk thermal conductivity of silicon (with nature isotope abundancy) in the open-source package SHENGBTE [START_REF] Li | SHENGBTE: A solver of the Boltzmann transport equation for phonons[END_REF]. The results at room temperature (300 K) calculated (based on the primitive unit cell) with a phonon wave vector 24 × 24 × 24 are, respectively, 120.69, 136.46, and 147.13 (W/m K) when the atomic interaction within the first, second, and third nearest-neighbor shells are separately considered in calculating the third-order FC. The thermal conductivity result is very close to the experimental value (148W/m K) when the third nearest-neighbor shell is considered, which demonstrates the good quality of DFT harmonic and anharmonic FCs. Due to the challenge of large computational cost, we consider only the first nearest-neighbor shell for the third-order FC in NEGF, which underestimates the bulk thermal conductivity at 300 K by about 20%. Since the present work is mainly focused on the demonstration of the anharmonic phonon NEGF methodology, such a simplified treatment captures the dominant third-order interaction and is acceptable from the perspective of microscopic modeling. The supercells in the DFT calculation of harmonic and anharmonic FCs have a limited size. The DFT FC matrices obtained for the supercells are then used to reconstruct the larger FC matrices for nanostructures as input into the large-scale phonon NEGF simulation. In our simulation, one conventional unit cell of silicon is chosen as one slab of the device in Fig. 3. The reconstruction procedure of both harmonic and third-order FC matrices for the device nanostructures is provided in Appendix E.

III. RESULTS AND DISCUSSIONS

In this section, we mainly aim to demonstrate the validity and application of the phonon NEGF formalism and computational framework introduced in Sec. II. Firstly, a simple case of ballistic heat transport at a Si/Ge interface is considered in Sec. III A. The quantitative validation of the anharmonic phonon NEGF formalism remains challenging due to its limited computational capability for extremely small structures and/or the usual empirical anharmonic FC input in previous works. In Sec. III B, we tackle this challenge by modeling anharmonic heat transport across a silicon thin film with a thickness larger than 10 nm through our DFTbased large-scale NEGF simulation. We discuss some crucial approximations in the treatment of anharmonic phonon scattering self-energy in Sec. III C. Finally in Sec. III D, we show an application of the anharmonic phonon NEGF to heat transport through a heterogeneous Si/Ge thin film and discuss some perspectives as well.

A. Validation: Ballistic heat transport

We consider ballistic heat transport through a Si/Ge interface as shown in Fig. 5(a), with only one unit cell (uc) in the device region. The lattice constant and harmonic FC of Ge are assumed the same as those of Si, with only the atomic mass difference taken into account, to be consistent with already reported studies [START_REF] Tian | Enhancing phonon transmission across a Si/Ge interface by atomic roughness: First-principles study with the Green's function method[END_REF][START_REF] Latour | Ab initio study of mode-resolved phonon transmission at Si/Ge interfaces using atomistic Green's functions[END_REF]. This approximate treatment is rather reasonable since the lattice structure of Ge is indeed very similar to that of Si. After a numerical convergence test, a frequency mesh of N m = 201 and transverse wave vector mesh of 20 × 20 are adopted in the present NEGF simulation with the effect of anharmonic phonon-phonon scattering turned off. Note that the mesh of both frequency and transverse wave vector for the present Si/Ge interface transport shall be much denser than that in homogeneous Si film to be discussed in Sec. III B. This is due to the very different cutoff frequencies (almost twofold difference) in the phonon bandstructures of Si and Ge. The ballistic transmission is calculated based on the Caroli formula [START_REF] Wang | Quantum thermal transport in nanostructures[END_REF][START_REF] Tian | Enhancing phonon transmission across a Si/Ge interface by atomic roughness: First-principles study with the Green's function method[END_REF]. The spectral transmission and transmissivity through the Si/Ge interface are shown in Figs. 5(b) and 5(c), respectively, which demonstrates an overall good agreement with respect to previous studies [START_REF] Tian | Enhancing phonon transmission across a Si/Ge interface by atomic roughness: First-principles study with the Green's function method[END_REF][START_REF] Latour | Ab initio study of mode-resolved phonon transmission at Si/Ge interfaces using atomistic Green's functions[END_REF]. The minor difference may arise from the slightly different DFT harmonic FCs calculated by different groups.

B. Validation: Anharmonic heat transport

In this subsection, we present a validation of the theoretical model of anharmonic scattering self-energy in the phonon NEGF formalism. A classical case of heat transport across a silicon thin film with a thickness d around T 0 is considered, as shown in Fig. 6. The temperatures of the left contact (contact 1) and the right contact (contact 2) are set at T 0 + T /2 and T 0 -T /2, respectively (with a temperature difference T). We consider cases of T 0 = 50, 100, 200, and 300 K, with the temperature difference T = 4, 4, 10, and 10 K respectively. The heat transport under such a temperature difference still lies within a linear regime based on our numerical test under a tiny temperature difference of 0.1 K. A frequency mesh of 101 and a transverse wave vector mesh 6 × 6 are adopted after a numerical convergence test at 300 K as summarized in Table III. Denser meshes produce thermal conductance within about 2% variation compared to the present mesh. We consider a series of thicknesses for the silicon thin film up to 24 uc (13 nm) due to computational cost limitations. For all the cases, the SCBA iterative solution of the governing equations in Sec. II B converges within ten iterations. Usually more iterations are needed for larger thickness or higher temperature due to stronger anharmonic phonon-phonon scattering. The computational time cost is about 1-36 h for all the cases in this subsection with 3636 CPUs. The position-dependent spectral heat flow and heat flow across the thin film with a thickness of 20 uc (10.8 nm) at 300 K are shown in Fig. 7, which demonstrates the good validity of heat flow conservation along the transport direction. After the convergence of SCBA iteration, the thermal conductance of the thin film is computed by G = J/(A c T ), with J the average heat flow across the thin film and A c the cross-section area of the device region (A c = a 2 • a 3 here). The effective thermal conductivity of the thin film is then calculated from the thermal conductance as κ eff = Gd.

To provide benchmark data to the anharmonic phonon NEGF simulation result, we conduct a Monte Carlo (MC) solution of the phonon Boltzmann equation for the same cross-plane heat transport with consistent DFT input. In principle, NEGF and Boltzmann formalisms should provide very similar results in the transport regime where the particle picture of phonons is valid and the coherent effects are negligible. For cross-plane heat transport through a silicon thin film at room temperature, a recent study has shown that the phonon Boltzmann equation could work down to a thickness of about 10 nm [START_REF] Hu | Unification of nonequilibrium molecular dynamics and the mode-resolved phonon Boltzmann equation for thermal transport simulations[END_REF]. This critical thickness is sound since the dominant phonon coherence length of silicon at room temperature is about 1 nm [START_REF] Chen | Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons[END_REF]. Thus for a thickness larger than 10 nm, the Monte Carlo solution could be a good benchmark for the anharmonic phonon NEGF result. We adopt an efficient energy-based deviational phonon MC scheme [START_REF] Péraud | Efficient simulation of multidimensional phonon transport using energy-based variance-reduced Monte Carlo formulations[END_REF] for solving the Boltzmann equation under the single-mode relaxation time (SMRT) approximation, which is a very good model for silicon with weak normal phonon scattering [START_REF] Guo | Phonon hydrodynamics and its applications in nanoscale heat transport[END_REF][START_REF] Li | SHENGBTE: A solver of the Boltzmann transport equation for phonons[END_REF]. The phonon dispersion and relaxation time are computed from the same DFT harmonic and anharmonic FCs as those for the phonon NEGF, and then converted into average isotropic ones [START_REF] Hua | Analytical Green's function of the multidimensional frequency-dependent phonon Boltzmann equation[END_REF] as input into the MC scheme. To treat the isothermal boundary conditions, the deviational particles emitted from the two contacts follow the Bose-Einstein equilibrium distribution [START_REF] Péraud | Efficient simulation of multidimensional phonon transport using energy-based variance-reduced Monte Carlo formulations[END_REF], which is consistent with the treatment in phonon NEGF by Eqs. ( 5) and [START_REF] Chen | Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons[END_REF]. We use 10 6 -5 × 10 6 deviational particles in the MC simulation at different film thicknesses and temperatures to ensure small statistical fluctuations. As a validation of the present MC code, the thickness-dependent cross-plane thermal conductivity of a silicon thin film at 300 K is shown in Fig. 8(a) when the third nearest-neighboring shell is considered in the third-order DFT FC. The MC result shows very good agreement with (i) the previous MC study considering DFT input [START_REF] Vermeersch | Cross-plane heat conduction in thin films with ab-initio phonon dispersions and scattering rates[END_REF]; (ii) a recent experimental measurement [START_REF] Jiang | Role of low-energy phonons with mean-free-paths >0.8 μm in heat conduction in silicon[END_REF]. Since only the first nearest-neighboring shell is considered in the third-order DFT FC input in the anharmonic phonon NEGF simulation, we also consider the same FC input in MC. The cross-plane thermal conductivity of a silicon thin film with a thickness of 20 uc (10.8 nm) and 24 uc (13 nm) at 300 K is shown in Table IV. It is seen that the result predicted by the present NEGF simulation is very close to the corresponding MC result. The ∼10% underestimation of the NEGF formalism comes from the fact that we still neglect some terms in the diagonal blocks and all the terms in the first off-diagonal blocks of the phonon scattering self-energy as explained in Sec. II C 1. Although these terms are small, including more terms in the scattering self-energy is expected to increase the thermal conductivity a bit, as will be inferred later in Sec. III C. Taking into account such effect, this comparison can act as a quantitative validation of the present anharmonic phonon NEGF formalism and numerical framework. As a comparison to the result by the scattering self-energy in previous anharmonic NEGF formalism in Refs. [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF], we conduct a simulation in the present numerical framework by including the corresponding factor of π/9 and 4 (cf. Sec. II A), respectively, before the right-hand side of Eq. ( 14). The silicon thin film with a thickness of 20 uc (10.8 nm) is simulated for the former comparison, and the thicknesses of 3 and 5 uc for the latter one. The scattering self-energy in Ref. [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF] overestimates the thermal conductance (867.20 MW/m 2 K) FIG. 8. Thickness-dependent cross-plane thermal conductivity of silicon thin film at 300 K: (a) comparison of the result by the present Monte Carlo (MC) with DFT input considering third nearest-neighboring shell in the third-order FC to previous result by MC with DFT input [START_REF] Vermeersch | Cross-plane heat conduction in thin films with ab-initio phonon dispersions and scattering rates[END_REF] and experimental data [START_REF] Jiang | Role of low-energy phonons with mean-free-paths >0.8 μm in heat conduction in silicon[END_REF]; (b) comparison of the present ballistic and anharmonic NEGF with DFT input with the present MC with DFT input considering first nearest-neighboring shell in the third-order FC (consistent with that in anharmonic NEGF). The dashed line represents the ballistic solution of the Boltzmann transport equation (BTE) with the same DFT input: κ eff = G ballistic d where the ballistic thermal conductance is computed by G ballistic = s vgx (q,s)>0 v gx (q, s) hω(q, s)∂ f eq /∂T dq/(2π ) 3 . for the considered two cases. The amount of underestimation is expected to increase with increasing film thickness due to stronger anharmonic scattering. In summary, the previous scattering self-energy expressions [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF] will underestimate or overestimate appreciably the phonon-phonon interaction strength. It is also interesting to note that the results of the phonon Boltzmann equation and the present phonon NEGF formalism also agree well for a film thickness smaller than 10 nm, as shown in Fig. 8(b). This indicates that the phonon Boltzmann equation seems to still work for cross-plane heat transport through a silicon thin film of a few nanometers. In principle, the phonon coherence shall play an important role in extremely small nanostructures, yet it does not manifest explicitly in the present cross-plane case. The underlying physical mechanism remains to be investigated in the future. Quantifying the coherence of phonons and its contribution to heat transport remains an open question. Since the computed thermal conductivities have some difference when the first and the third nearest-neighboring shell are, respectively, considered, we show the thickness-dependent cross-plane thermal conductivity normalized by the corresponding bulk value in Fig. 9. At small thickness, the present anharmonic phonon NEGF result agrees well with the general trend of nondimensional thermal conductivity from different reference data.

We also provide the results of local temperature distribution across the silicon thin film with a thickness of 20 uc (10.8 nm) and 24 uc (13 nm) at 300 K in Fig. 10. The results by the present anharmonic phonon NEGF simulation show a good agreement with those of the MC simulation with the same DFT input. A large temperature jump near the two contacts is obtained due to the strong nonequilibrium effect between the thin film and contacts, which is well known in the heat transport community and has been already reported [START_REF] Chen | Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons, Molecules, Phonons, and Photons[END_REF][START_REF] Zhang | Nano/Microscale Heat Transfer[END_REF]. This represents a further demonstration of the accuracy of the present formalism and numerical framework. Moreover, the large deviation of local temperature distribution from the result in the ballistic limit indicates that the phonon-phonon scattering is appreciable even at thicknesses around 10 nm. This is also clearly visible in the difference between ballistic and anharmonic results in Fig. 8(b), where about 20% reduction of thermal conductivity is obtained due to anharmonic phonon-phonon scattering. Physically this reduction mainly comes from the scattering of optical phonons and LA (longitudinal acoustic) phonons, as is inferred from the spectral heat flow distribution in Fig. 11.

Finally, we demonstrate the temperature-dependent results for the thin film with a thickness of 20 uc (10.8 nm) in Fig. 12. With decreasing temperature, the phonon-phonon scattering rate will decrease due to smaller amplitude of lattice atomic vibration and thus weaker anharmonicity. As a result, the reduction of thermal conductivity due to anharmonic scattering will also decrease with temperature, as shown in Fig. 12(a). Below 100 K, the heat transport across the silicon thin film becomes almost ballistic. This is also seen in the dimensionless temperature distribution in Fig. 12(b). As the phonon scattering weakens at lower temperature, the thermalization within the thin film is more difficult. Thus larger temperature jump occurs near the contacts and the temperature distribution turns to be closer to the uniform profile in the ballistic limit. The present anharmonic phonon NEGF formalism generally captures well the behaviors of heat transport from the low temperature regime, wherein the quantum effect is important, to the room temperature regime. It still works well at higher temperature, where the anharmonic phonon scattering becomes even stronger. We do not discuss such situation for this case due to the very large computational cost.

C. Assessment of approximations in numerical implementation

In this subsection, we discuss the numerical approximations implemented for the calculation of the anharmonic phonon scattering self-energy, which is crucial to ensure a balance between accuracy and efficiency in large-scale quantum heat transport simulation. As already mentioned in Sect. II C 1, a simpler treatment of the scattering self-energy matrix, by considering only the 3 × 3 diagonal blocks, is adopted in a previous work by one of the coauthors [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF]. Here we also include this approximation into our computational framework and compare the obtained results with those of the present approximation, i.e., considering the dominant terms in the N x × N x diagonal blocks (cf. Sec. II C 1). The heat transport across a silicon thin film with two thicknesses of 3 and 5 uc at room temperature is simulated. For the thickness of 3 uc, the thermal conductance predicted based on the present approximation and the previous one is, respectively, 939.72 and 867.65 (MW/m 2 K). For the thickness of 5 uc, the thermal conductance based on the present and previous approximation is, respectively, 890.97 and 784.09 (MW/m 2 K). Comparing to the present approximation, the 3 × 3 diagonal approximation will underestimate ∼8% and ∼12% of the thermal conductance of silicon thin film with a thickness of 3 and 5 uc, respectively. The underestimation will increase with increasing film thickness due to the stronger phonon-phonon scattering rate. The underlying reason can be understood from the spectral heat flow results shown in Fig. 13. It is seen that considering only the 3 × 3 diagonal blocks of the anharmonic scattering self-energy matrix will significantly overestimate the scattering rate of low-frequency TA (transverse acoustic) phonons which mainly contribute to heat transport. In principle, the incorporation of more terms into the scattering self-energy matrix in the present approximation will help in retrieving the underestimation of thermal conductivity shown in Table IV and Fig. 8(b). However, since the dominant terms have already been considered, the present treatment constitutes a good balance between the computational accuracy and efficiency. Note that a further approximation was made to the product of the greater/lesser Green's function matrix in the phonon scattering self-energy of Eq. ( 8) [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF]: where the matrix product is simplified into the product of their 3 × 3 diagonal blocks. Such kind of local diagonal approximation is usually assumed in electron NEGF [START_REF] Luisier | Atomistic full-band simulations of silicon nanowire transistors: Effects of electron-phonon scattering[END_REF][START_REF] Rhyner | Atomistic modeling of coupled electron-phonon transport in nanowire transistors[END_REF] and will speed up the simulation and reduce the computational cost. Yet based on our numerical test, it is a too large simplification for phonon NEGF which will lead to a further overestimation of the phonon-phonon scattering rate and will make the SCBA iteration diverge. The physical interpretation may be that the self-energy for electron-phonon scattering is linearly proportional to the electron Green's function, whereas the self-energy for phonon-phonon scattering is a quadratic function of the phonon Green's function, which is a more complicated nonlinear problem. This could be a possible explanation for the large scaling of empirical valence-force-field anharmonic FCs to fit the experimental thermal conductivity of silicon therein [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF].

>, <i j s, ll (ω) = 1 2 i h l 1 l 2 j 1 j 2 j 3 j 4 ∞ -∞ dω 2π i j 1 j 2 ll 1 l 2 j j 3 j 4 ll 2 l 1 × G >, < j 1 j 4 l 1 l 1 (ω )G >, < j 2 j 3 l 2 l 2 (ω -ω ), (30)

D. Application and perspective

The present nonequilibrium Green's function formalism is specially applicable for nanostructures with broken translational symmetry, where the usual phonon picture and Boltzmann transport theory do not work any more. In this subsection, we show a NEGF simulation of heat transport through a multilayer Si/Ge heterogeneous thin film with a thickness of only 3 uc, shown in Fig. 14(a). The embedded Ge layer plays the role of a potential barrier to the heat conduction. Both the second-and third-order FCs of the Ge are assumed to be the same as those of Si with only the atomic mass difference taken into account [START_REF] Dai | Rigorous formalism of anharmonic atomistic Green's function for three-dimensional interfaces[END_REF]. A frequency mesh of 121 and a transverse wave vector mesh of 8 × 8 are adopted in the simulation. In contrast to the usual reduction effect, here the anharmonic phonon-phonon scattering enhances appreciably the thermal conductance compared to the ballistic limit, as shown in Fig. 14(b). Such an enhancement increases with increasing temperature, and reaches close to 20% at 600 K. The underlying physical mechanism is that the anharmonicity opens transport channels across the Ge barrier layer via the interaction between phonons with different frequencies. The spectral contribution of heat flow in Fig. 14(c) elucidates indeed stronger interaction within 2-6 THz and 8-12 THz, which correspond to the frequency range with good overlap of local density of states in the Si layer and Ge layer, as given in Fig. 14(d).

Another important case where the phonon Boltzmann equation fails and the present anharmonic phonon NEGF method becomes indispensable is the heat conduction through superlattices with small periods. The anharmonic phononphonon scattering plays a crucial part in destroying the coherence of phonons, especially when the total thickness is large and the system temperature is high [START_REF] Luckyanova | Coherent phonon heat conduction in superlattices[END_REF], which was instead modeled by ballistic phonon NEGF [START_REF] Tian | Green's function studies of phonon transport across Si/Ge superlattices[END_REF] due to the limitation of methodology development. The present numerical framework provides an avenue to modeling such partially coherent heat transport. As this work is mainly focused on the methodology development and demonstration, the discussion of the superlattices case is beyond the scope of the current article and is pending in our future work. Note that recently a generalized transport equation based on the Wigner distribution function [START_REF] Simoncelli | Unified theory of thermal transport in crystals and glasses[END_REF] or Green-Kubo formalism [START_REF] Isaeva | Modeling heat transport in crystals and glasses from a unified lattice-dynamical approach[END_REF] beyond the Boltzmann equation is developed to model heat conduction in crystals and glasses in a unified way. One difference between this generalized transport model and the anharmonic NEGF is that the former still relies on a unit cell of the system to compute its eigenstates and eigenenergies as input. The present NEGF formalism only requires the FC matrix of the system, which may not have any translational periodicity to define a unit cell such as the case in Fig. 14(a) and the case of aperiodic superlattices [START_REF] Hu | Machine-Learning-Optimized Aperiodic Superlattice Minimizes Coherent Phonon Heat Conduction[END_REF]. In addition, the generalized transport model [START_REF] Simoncelli | Unified theory of thermal transport in crystals and glasses[END_REF][START_REF] Isaeva | Modeling heat transport in crystals and glasses from a unified lattice-dynamical approach[END_REF] is currently designed for the calculation of bulk thermal conductivity of crystals and glasses. To directly model heat conduction in finite-size nanostructures, the boundary treatment remains a challenging task due to the difficulty in specifying the off-diagonal terms of the generalized distribution function matrix. In comparison, it is quite convenient in NEGF to treat both the boundary scattering and anharmonic scattering via the self-energy concept.

IV. CONCLUSIONS

In summary, we present a nonequilibrium Green's function (NEGF) computational framework with the first-principle input for large-scale quantum heat transport simulations accounting for the anharmonic phonon-phonon scattering. The theoretical formulation of anharmonic scattering self-energy is clarified through a careful diagrammatic perturbation analysis, with a Fourier's representation further introduced, which satisfies both the energy and momentum conservation for nanostructures with transverse periodicity. A quantitative validation of the anharmonic phonon NEGF formalism is demonstrated through a comparison to the results by the Boltzmann equation in the particle transport regime based on a classical cross-plane heat transport in a silicon thin film. The phonon-phonon scattering is shown to be non-negligible even for films with a thickness as small as 10 nm and to introduce a 20% reduction of thermal conductivity at room temperature. The widely used local diagonal approximation of the scattering self-energy in electron NEGF is shown to significantly overestimate the phonon-phonon interaction strength. In heat conduction through an extremely thin heterogeneous film, the anharmonic phonon scattering opens transport channels and enhances the thermal conductance up to around 20% at 600 K. The present computational framework provides a potential platform for first-principle prediction of thermal boundary conductance at the interface, and for investigation of the transition from coherent to incoherent heat transport in nanophononic crystals such as superlattices, which is pending in the near future. This study thus opens opportunities for the understanding and tuning of phonon heat transport in the coherent quantum wave regime.

For heat transport in 1D nanostructures without any translational periodicity, the total Hamiltonian operator of the system in the Heisenberg representation is expressed as [START_REF] Ziman | Electrons and Phonons: The Theory of Transport Phenomena in Solids[END_REF]

H (t ) = H 0 (t ) + V (t ), ( A1
)
where H 0 (t ) and V(t) are the exactly solvable harmonic part and the third-order anharmonic perturbation part, respectively,

H 0 (t ) = 1 2 n,i ui n (t ) ui n (t ) + 1 2 nm i j i j nm u i n (t )u j m (t ), (A2) V (t ) = 1 3! nml i jk i jk nml u i n (t )u j m (t )u k l (t ), ( A3
)
where n, m, l denote the atomic index, whereas i, j, k denote the Cartesian coordinate index (x, y, z). u i n (t ) is the atomic displacement operator rescaled with atomic mass (u i n = r i n √ M n , with r i n the real atomic displacement operator and M n the atomic mass). The dot on the atomic displacement operator means the time derivative (corresponding to atomic velocity). The normalized second-order and third-order force constants i j nm and i jk nml are defined as

i j nm = 1 √ M n M m ∂ 2 E ∂r i n ∂r j m 0 , ( A4
) i jk nml = 1 √ M n M m M l ∂ 3 E ∂r i n ∂r j m ∂r k l 0 , ( A5
)
where E denotes the atomic interaction potential, and the subscript "0" denotes the equilibrium position of the crystal lattice.

The starting point of the diagrammatic perturbation expansion is the definition of the contour-ordered phonon Green's function [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Wang | Nonequilibrium Green's function approach to mesoscopic thermal transport[END_REF]:

G i j nm (τ, τ ) = -i T C u i n (τ )u j m (τ ) , ( A6
)
where T C is the contour ordering operator, and the times τ, τ are on the contour. Note that in the present derivation we adopt the convention of h = 1 [START_REF] Mahan | Many-Particle Physics[END_REF], which shall be recovered in the final expression. Transforming from the Heisenberg representation to the interaction representation of Eq. (A6), we obtain the following expression for the contour-ordered phonon Green's function: ρ 0 = exp(-βH 0 )/Tr[exp(-βH 0 )], where β = 1/k B T . The caret on the operators denotes the interaction representation. The evolution operator in Eq. (A7) is defined as [START_REF] Rammer | Quantum Field Theory of Non-equilibrium States[END_REF][START_REF] Haug | Quantum Kinetics in Transport and Optics of Semiconductors[END_REF]

G i j nm (τ, τ ) = -i T C S C ûi n (τ ) û j m (τ ) 0 , ( A7
S C = exp -i C dτ 1 V (τ 1 ) = ∞ n=0 (-i) n n! C dτ 1 C dτ 2 • • • C dτ n V (τ 1 ) V (τ 2 ) • • • V (τ n ), (A8)
where the third-order perturbation part of the Hamiltonian in the interaction representation is

V (τ ) = 1 3! nml i jk i jk nml ûi n (τ ) û j m (τ ) ûk l (τ ). (A9)
The evolution operator is the basis for the diagrammatic perturbation expansion. Considering the infinite series in Eq. (A8) within second order and substituting into Eq. (A7), we obtain

G i j nm (τ, τ ) = G (0)i j nm (τ, τ ) + G (1)i j nm (τ, τ ) + G (2)i j nm (τ, τ ), ( A10
)
where the zeroth-, first-, and second-order terms are, respectively,

G (0)i j nm (τ, τ ) = -i T C ûi n (τ ) û j m (τ ) 0 , (A11) G (1)i j nm (τ, τ ) = (-i) 2 C dτ 1 T C ûi n (τ ) V (τ 1 ) û j m (τ ) 0 , (A12) G (2)i j nm (τ, τ ) = (-i) 3 2 C dτ 1 C dτ 2 T C ûi n (τ ) V (τ 1 ) V (τ 2 ) û j m (τ ) 0 . (A13)
The first-order term in Eq. (A12) will be vanishing since it involves an odd number of atomic displacement operators. The full expression of the second-order term will be

G (2)i j nm (τ, τ ) = (-i) 3 2 C dτ 1 C dτ 2 1 3! • 3! n 1 m 1 l 1 i 1 j 1 k 1 n 2 m 2 l 2 i 2 j 2 k 2 i 1 j 1 k 1 n 1 m 1 l 1 i 2 j 2 k 2 n 2 m 2 l 2 × T C ûi n (τ ) ûi 1 n 1 (τ 1 ) û j 1 m 1 (τ 1 ) ûk 1 l 1 (τ 1 ) ûi 2 n 2 (τ 2 ) û j 2 m 2 (τ 2 ) ûk 2 l 2 (τ 2 ) û j m (τ ) 0 . (A14)
Wick's theorem [START_REF] Rammer | Quantum Field Theory of Non-equilibrium States[END_REF][START_REF] Mahan | Many-Particle Physics[END_REF] will be applied for the decomposition of the expectation of the product of eight displacement operators, where only the connected diagrams shown in Fig. 15 are considered. These represent the physically feasible three-phonon anharmonic scattering process. There are 3 × 3 × 2 pairing combinations for each case in Fig. 15, and totally 36 equivalent pairing combinations such that Eq. (A14) becomes

G (2)i j nm (τ, τ ) = (-i) 3 2 C dτ 1 C dτ 2 36 3! • 3! n 1 m 1 l 1 i 1 j 1 k 1 n 2 m 2 l 2 i 2 j 2 k 2 i 1 j 1 k 1 n 1 m 1 l 1 i 2 j 2 k 2 n 2 m 2 l 2 × T C ûi n (τ ) ûi 1 n 1 (τ 1 ) 0 T C û j 1 m 1 (τ 1 ) û j 2 m 2 (τ 2 ) 0 T C ûk 1 l 1 (τ 1 ) ûi 2 n 2 (τ 2 ) 0 T C ûk 2 l 2 (τ 2 ) û j m (τ ) 0 , (A15)
which can be rewritten into

G (2)i j nm (τ, τ ) = i 2 C dτ 1 C dτ 2 n 1 m 1 l 1 i 1 j 1 k 1 n 2 m 2 l 2 i 2 j 2 k 2 i 1 j 1 k 1 n 1 m 1 l 1 i 2 j 2 k 2 n 2 m 2 l 2 G (0)ii 1 nn 1 (τ, τ 1 )G (0) j 1 j 2 m 1 m 2 (τ 1 , τ 2 )G (0)k 1 i 2 l 1 n 2 (τ 1 , τ 2 )G (0)k 2 j l 2 m (τ 2 , τ ). ( A16
)
Putting Eqs. (A11) and (A16) into Eq. (A10), we obtain the contour-order phonon Green's function within second order as

G i j nm (τ, τ ) = G (0)i j nm (τ, τ ) + i 2 C dτ 1 C dτ 2 n 1 m 1 l 1 i 1 j 1 k 1 n 2 m 2 l 2 i 2 j 2 k 2 i 1 j 1 k 1 n 1 m 1 l 1 i 2 j 2 k 2 n 2 m 2 l 2 × G (0)ii 1 nn 1 (τ, τ 1 )G (0) j 1 j 2 m 1 m 2 (τ 1 , τ 2 )G (0)k 1 i 2 l 1 n 2 (τ 1 , τ 2 )G (0)k 2 j l 2 m (τ 2 , τ ). ( A17
)
Comparing Eq. (A17) to Dyson's equation it follows [START_REF] Rammer | Quantum Field Theory of Non-equilibrium States[END_REF][START_REF] Haug | Quantum Kinetics in Transport and Optics of Semiconductors[END_REF] that

G i j nm (τ, τ ) = G (0)i j nm (τ, τ ) + C dτ 1 C dτ 2 n 1 l 2 i 1 k 2 G (0)ii 1 nn 1 (τ, τ 1 ) i 1 k 2 s, n 1 l 2 (τ 1 , τ 2 )G k 2 j l 2 m (τ 2 , τ ), (A18)
and we get the expression of anharmonic phonon scattering self-energy as

i j s, nm (τ 1 , τ 2 ) = i 2 m 1 l 1 n 2 m 2 j 1 k 1 i 2 j 2 i j 1 k 1 nm 1 l 1 ji 2 j 2 mn 2 m 2 G j 1 j 2 m 1 m 2 (τ 1 , τ 2 )G k 1 i 2 l 1 n 2 (τ 1 , τ 2 ), (A19)
where the unperturbed Green's function has been replaced by the full Green's function based on the self-consistent Born approximation. Note that the factor of 4 in the scattering self-energy expression in Ref. [START_REF] Luisier | Atomistic modeling of anharmonic phonon-phonon scattering in nanowires[END_REF] with respect to the present expression comes from the repeated count of the combinations in Eq. (A15) because the middle two pairings such as

T C [û j 1 m 1 (τ 1 )û j 2 m 2 (τ 2 )] 0 are distinguished from T C [û j 2 m 2 (τ 2 )û j 1 m 1 (τ 1 )] 0 [92]
. Such two terms are actually the same since the definition of the contour-order phonon Green's function is irrelevant to the order of atomic displacement operators. Regarding the factor of π /9 in the scattering self-energy in Ref. [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF] with respect to the present one, it is difficult to clarify the reason for this difference since we cannot find any details of the derivation in the literature [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF].

Through the analytic continuation process based on the Langreth theorem [START_REF] Haug | Quantum Kinetics in Transport and Optics of Semiconductors[END_REF], we get the greater/lesser scattering self-energy in real time from Eq. (A19):

>, <i j s, nm (t 1 , t 2 ) = i 2 m 1 l 1 n 2 m 2 j 1 k 1 i 2 j 2 i j 1 k 1 nm 1 l 1 ji 2 j 2 mn 2 m 2 G >, < j 1 j 2 m 1 m 2 (t 1 , t 2 )G >, <k 1 i 2 l 1 n 2 (t 1 , t 2 ). (A20)
For the stationary state heat transport considered in this work, the dependence on t 1 and t 2 will be reduced to the dependence on (t 1 -t 2 ). For convenience, the following Fourier transform and its inverse transform are introduced [START_REF] Mahan | Many-Particle Physics[END_REF]:

>, <i j s, nm (t 1 , t 2 ) = >, <i j s, nm (t 1 -t 2 ) = ∞ -∞ dω 2π exp [-i(t 1 -t 2 )ω] >, <i j s, nm (ω), (A21) >, <i j s, nm (ω) = ∞ -∞ dt 1 exp [iω(t 1 -t 2 )] >, <i j s, nm (t 1 -t 2 ). ( A22
)
Fourier transforms of Eqs. (A21) and (A22) are also applicable to the Green's function in Eq. (A20). After the Fourier transform, Eq. (A20) becomes

>, <i j s, nm (ω) = i 2 m 1 l 1 n 2 m 2 j 1 k 1 i 2 j 2 ∞ -∞ dω 2π i j 1 k 1 nm 1 l 1 ji 2 j 2 mn 2 m 2 G >, < j 1 j 2 m 1 m 2 (ω )G >, <k 1 i 2 l 1 n 2 (ω -ω ), (A23)
which is exactly Eq. ( 8) in the main text. It is also seen from Eq. (A23) that the energy conservation is automatically satisfied in the three-phonon scattering process: ω = ω + (ω -ω ).

APPENDIX B: FOURIER'S REPRESENTATION OF THE ANHARMONIC PHONON-PHONON SCATTERING SELF-ENERGY MATRIX

The starting point of the derivation is the expression of anharmonic phonon scattering self-energy matrix Eq. ( 8) for 1D nanostructures:

>, <i j s, ll (ω) = 1 2 i h l 1 l 2 l 3 l 4 j 1 j 2 j 3 j 4 ∞ -∞ dω 2π i j 1 j 2 ll 1 l 2 j j 3 j 4 l l 3 l 4 G >, < j 1 j 4 l 1 l 4 (ω )G >, < j 2 j 3 l 2 l 3 (ω -ω ). (B1)
The inverse transform of the Fourier's representation, Eq. ( 16), for the anharmonic FC matrix is

i j 1 j 2 ll 1 l 2 = 1 N 2 q ⊥ q ⊥ ˜ i j 1 j 2 l x l 1x l 2x (q ⊥ , q ⊥ ) exp (iq ⊥ • R ⊥ ) exp (iq ⊥ • R ⊥ ). (B2)
For convenience of later derivation, Eq. (B2) is slightly rewritten into

i j 1 j 2 ll 1 l 2 = 1 N 2 q ⊥ q ⊥ ˜ i j 1 j 2 l x l 1x l 2x (q ⊥ , q ⊥ ) exp [i(q ⊥ + q ⊥ ) • R ⊥ ] exp iq ⊥ • R ⊥ , (B3) with R ⊥ = R ⊥ -R ⊥ = (l 1y -l 2y )a 2 + (l 1z -l 2z )a 3 .
The Fourier's representation of the scattering self-energy, Green's function, and third-order FC matrix in Eq. (B1) can be expressed as follows:

<, >i j s, ll (ω) = 1 N q exp (iq • R) <, >i j s, l x l x (ω; q), (B4) G <, > j 1 j 4 l 1 l 4 (ω ) = 1 N q 1 exp (iq 1 • R 1 )G <, > j 1 j 4 l 1x l 4x (ω ; q 1 ), (B5) G <, > j 2 j 3 l 2 l 3 (ω -ω ) = 1 N q 2 exp (iq 2 • R 2 )G <, > j 2 j 3 l 2x l 3x (ω -ω ; q 2 ), (B6) i j 1 j 2 ll 1 l 2 = 1 N 2 q 3 q 4 ˜ i j 1 j 2 l x l 1x l 2x (q 3 , q 4 ) exp [i(q 3 + q 4 ) • R 3 ] exp (iq 4 • R 4 ), ( B7
) j j 3 j 4 l l 3 l 4 = 1 N 2 q 5 q 6 ˜ j j 3 j 4 l x l 3x l 4x (-q 5 , -q 6 ) exp [i(q 5 + q 6 ) • R 5 ] exp (iq 6 • R 6 ), (B8)
where all the subscripts "⊥" are omitted to keep a simple notation, and the relative displacements are defined as follows: R = (l yl y )a 2 + (l zl z )a 3 , (B9)

R 1 = (l 1y -l 4y )a 2 + (l 1z -l 4z )a 3 , (B10) R 2 = (l 2y -l 3y )a 2 + (l 2z -l 3z )a 3 , (B11) R 3 = (l y -l 1y )a 2 + (l z -l 1z )a 3 , (B12) R 4 = (l 1y -l 2y )a 2 + (l 1z -l 2z )a 3 , (B13) R 5 = (l 3y -l y )a 2 + (l 3z -l z )a 3 , (B14) R 6 = (l 4y -l 3y )a 2 + (l 4z -l 3z )a 3 . (B15)
Substitution of Eqs. (B4)-(B8) into Eq. (B1) gives rise to

1 N q exp (iq • R) <, >i j s, l x l x (ω; q) = i h 2 l 1 l 2 l 3 l 4 j 1 j 2 j 3 j 4 1 N 2 q 3 q 4 ˜ i j 1 j 2 l x l 1x l 2x (q 3 , q 4 ) exp [i(q 3 + q 4 ) • R 3 ] exp (iq 4 • R 4 ) × 1 N 2 q 5 q 6 ˜ j j 3 j 4 l x l 3x l 4x (-q 5 , -q 6 ) exp [i(q 5 + q 6 ) • R 5 ] exp (iq 6 • R 6 ) × ∞ -∞ dω 2π 1 N q 1 exp (iq 1 • R 1 )G <, > j 1 j 4 l 1x l 4x (ω ; q 1 ) 1 N q 2 exp (iq 2 • R 2 )G <, > j 2 j 3 l 2x l 3x (ω -ω ; q 2 ). (B16)
Since R = R 3 + R 4 + R 2 + R 5 , multiplying exp(-iq • R) on both sides of Eq. (B16) and summing over (l y , l z ), we re-form the left-hand side term as

1 N q R exp [i(q -q ) • R] <,>i j s, l x x (ω; q) = <,>i j s, l x l x (ω; q ), ( B17
)
where a classical relation in lattice dynamics has been used [START_REF] Ziman | Electrons and Phonons: The Theory of Transport Phenomena in Solids[END_REF][START_REF] Maradudin | Theory of Lattice Dynamics in the Harmonic Approximation[END_REF]:

R exp [i(q -q ) • R] = N (q -q ). (B18)
The right-hand side term becomes

i h 2 l 1 l 2 l 3 l 4 j 1 j 2 j 3 j 4 1 N 2 q 3 q 4 ˜ i j 1 j 2 l x l 1x l 2x (q 3 , q 4 ) R 3 exp [i(q 3 + q 4 -q ) • R 3 ] exp [i(q 4 -q ) • R 4 ] × 1 N 2 q 5 q 6 ˜ j j 3 j 4 l x l 3x l 4x (-q 5 , -q 6 ) exp [i(q 5 + q 6 -q ) • R 5 ] exp (iq 6 • R 6 ) × ∞ -∞ dω 2π 1 N q 1 exp (iq 1 • R 1 )G <, > j 1 j 4 l 1x l 4x
(ω ; q 1 ) 1 N q 2 exp [i(q 2q ) • R 2 ]G <, > j 2 j 3 l 2x l 3x (ω -ω ; q 2 ) = i h 2 l 1 l 2 l 3 l 4 j 1 j 2 j 3 j 4 1 N q 4 ˜ i j 1 j 2 l x l 1x l 2x (qq 4 , q 4 ) exp [i(q 4q ) • R 4 ] × 1 N 2 q 5 q 6 ˜ j j 3 j 4 l x l 3x l 4x (-q 5 , -q 6 ) exp [i(q 5 + q 6q

) • R 5 ] exp (iq 6 • R 6 ) × ∞ -∞ dω 2π 1 N q 1 exp (iq 1 • R 1 )G <, > j 1 j 4 l 1x l 4x
(ω ; q 1 ) 1 N q 2 exp [i(q 2q ) • R 2 ]G <, > j 2 j 3 l 2x l 3x (ω -ω ; q 2 ). (B19)

In the derivation of Eq. (B19), we have identified the sum over (l y , l z ) as a sum over R 3 , and also used the relation Eq. (B18). Therefore, the scattering self-energy is <, >i j s, l x l x (ω; q ) = i h 2 l 1 l 2 l 3 l 4 j 1 j 2 j 3 j 4 1 N q 4 ˜ i j 1 j 2 l x l 1x l 2x (qq 4 , q 4 ) exp [i(q 4q ) • R 4 ] × 1 N 2 q 5 q 6 ˜ j j 3 j 4 l x l 3x l 4x (-q 5 , -q 6 ) exp [i(q 5 + q 6q

) • R 5 ] exp (iq 6 • R 6 ) × ∞ -∞ dω 2π 1 N q 1 exp (iq 1 • R 1 )G <, > j 1 j 4 l 1x l 4x
(ω ; q 1 ) 1 N q 2 exp [i(q 2q ) • R 2 ]G <, > j 2 j 3 l 2x l 3x (ω -ω ; q 2 ).

(B20)

Since only R 5 depends on (l y , l z ), summing both sides of Eq. (B20) over (l y , l z ) (identified as a sum over R 5 ) and using again the relation Eq. (B18), we could reduce Eq. (B20) to N <, >i j s, l x l x (ω; q ) = i h 2 l 1 l 2 l 3 l 4 j 1 j 2 j 3 j 4 1 N q 4 ˜ i j 1 j 2 l x l 1x l 2x (qq 4 , q 4 ) exp [i(q 4q ) • R 4 ]

1 N q 6 ˜ j j 3 j 4 l x l 3x l 4x (q 6q , -q 6 ) exp (iq 6 • R 6 )

× ∞ -∞ dω 2π 1 N q 1 exp (iq 1 • R 1 )G <, > j 1 j 4 l 1x l 4x
(ω ; q 1 ) 1 N q 2 exp [i(q 2q ) • R 2 ]G <, > j 2 j 3 l 2x l 3x (ω -ω ; q 2 ). (B21)

Since R 4 = R 1 -R 2 + R 6 , Eq. ( B21) is reformulated into N <, >i j s, l x l x (ω; q ) = i h 2 l 1 l 2 l 3 l 4 j 1 j 2 j 3 j 4 1 N q 4 ˜ i j 1 j 2 l x l 1x l 2x (qq 4 , q 4 ) 1 N q 6 ˜ j j 3 j 4 l x l 3x l 4x (q 6q , -q 6 ) exp [i(q 6 + q 4q

) • R 6 ] × ∞ -∞ dω 2π
1 N q 1 exp [i(q 1 + q 4q ) • R 1 ]G <, > j 1 j 4 l 1x l 4x

(ω ; q 1 ) × 1 N q 2 exp [i(q 2q 4 ) • R 2 ]G <, > j 2 j 3 l 2x l 3x (ω -ω ; q 2 ). (B22)

In this way, the sum over (l 1y , l 1z ) (identified as a sum over R 1 ) and the sum over (l 2y , l 2z ) (identified as a sum over R 2 ) become independent, and the relation Eq. (B18) is used to remove the dependence on R 1 and R 2 . Then the sum over (l 3y , l 3z ) (identified as a sum over R 6 ) becomes independent and the relation Eq. (B18) is used again. Finally, Eq. (B22) will be reduced to the following expression: <, >i j s, l x l x (ω; q ) = i h 2 l 1x l 2x l 3x l 4x j 1 j 2 j 3 j 4 1 N q 4 ˜ i j 1 j 2 l x l 1x l 2x (qq 4 , q 4 ) ˜ j j 3 j 4 l x l 3x l 4x (-q 4 , q 4q ) × ∞ -∞ dω 2π G <, > j 1 j 4 l 1x l 4x (ω ; qq 4 )G <, > j 2 j 3 l 2x l 3x (ω -ω ; q 4 ). (B23)

Introducing a variable change, q = qq 4 , and rearranging the notation (q → q, q → q ), we get the final expression of anharmonic phonon-phonon scattering self-energy in the Fourier's representation as <, >i j s, l x l x (ω; q) = 1 2 i h l 1x l 2x l 3x l 4x j 1 j 2 j 3 j 4

1 N q ˜ i j 1 j 2 l x l 1x l 2x (q , qq ) ˜ j j 3 j 4 l x l 3x l 4x (qq, -q ) × ∞ -∞ dω 2π G <, > j 1 j 4 l 1x l 4x (ω ; q )G <, > j 2 j 3 l 2x l 3x (ω -ω ; qq ). (B24) FIG. 17. Reconstruction of the harmonic FC matrices as input into the phonon NEGF simulation: (a) schematic of the supercell in DFT calculations, and the device simulation domain, where the solid cuboid represents the central device region whereas the dashed cuboid represents the transverse periodic units that have interaction with the device region with the interaction FC matrix denoted by 0 ∼ 8 ; (b) reconstruction of the harmonic FC matrix 4 within the central region from the elementary harmonic interaction matrices in the DFT supercell; (c) reconstruction of the representative harmonic FC matrix 0 (between the central region and the transverse periodic unit) from the elementary harmonic interaction matrices in the DFT supercell. Other harmonic FC matrices can be reconstructed similarly. Each cube represents a conventional unit cell. harmonic FC matrices within the central device region, and between the device region and the transverse periodic units 0 -8 are reconstructed from the elementary interaction FC matrices within the supercell. The harmonic FC matrix in Eq. ( 11) is then computed based on Eq. [START_REF] Maire | Heat conduction tuning by wave nature of phonons[END_REF] as

˜ (q ⊥ ) = 4 + i=0, 8 i exp [∓i(a 2 + a 3 ) • q ⊥ ] + i=1, 7 i exp (∓ia 2 • q ⊥ ) + i=3, 5 i exp (∓ia 3 • q ⊥ ) + i=2, 6 i exp [∓i(a 2 -a 3 ) • q ⊥ ], (E1)
where the "-" and "+" signs, respectively, correspond to the first and second index of the FC matrix in the summation. The schematic procedure of the reconstruction of the third-order FC matrix is shown in Fig. 18. The treatment is slightly different because the atomic interaction range for third-order FC is very short. From the elementary interaction matrix in the supercell, we extract the third-order FCs for each atom in a unit cell and its neighboring atoms, i.e., a basic third-order FC matrix ˜ i j 1 j 2 l x l 1x l 2x , 0 (q ⊥ , q ⊥ ) for one unit cell (one slab) is obtained. When any of the neighboring atoms lies within the transverse periodic unit cell, a phase factor shall be added based on the Fourier's representation of a third-order FC matrix in Eq. ( 16). The third-order FC matrix for the whole device is constructed simply by repeating ˜ i j 1 j 2 l x l 1x l 2x , 0 (q ⊥ , q ⊥ ) for times equal to the number of slabs. Note that for heterogeneous material properties, the procedure of the reconstruction of device FC matrices shall be slightly adapted to consider the local variation of atomic interaction, for instance, in a recent DFT-based ballistic phonon NEGF modeling of disordered lithiated molybdenum disulfide (MoS 2 ) nanostructures [START_REF] Bunjaku | Thermal properties of disordered Li x MoS 2 : An ab initio study[END_REF]. FIG. 18. Reconstruction of the third-order anharmonic FC matrix as input into the phonon NEGF simulation: (a) schematic of the supercell in DFT calculation, and the device simulation domain with the corresponding third-order FC matrix ˜ i j 1 j 2 lx l 1x l 2x (q ⊥ , q ⊥ ); (b) reconstruction of the third-order FC matrix from the elementary third-order FC matrix for one unit cell ˜ i j 1 j 2 lx l 1x l 2x , 0 (q ⊥ , q ⊥ ) (as explained in the main text) extracted from the DFT supercell. Each cube represents a conventional unit cell, and the parallelpiped denotes the device region.

FIG. 1 .

 1 FIG. 1. Schematic of the physical model in phonon NEGF formalism: (a) 1D nanostructure; (b) 3D nanostructure with transverse periodicity.

FIG. 2 .

 2 FIG. 2. Feynman diagram for three-phonon anharmonic scattering in 3D nanostructures with transverse periodicity.

FIG. 4 .

 4 FIG. 4. Schematic of the discretization of (a) frequency and (b) transverse wave vector in phonon NEGF simulation. The frequency interval is computed as ω = 2ω m /(N m -1), with ω m the maximum crystal phonon frequency and N m the number of frequency points. The transverse reciprocal lattice vectors are defined as b i • a j = 2πδ i j , where i, j = 2, 3.

FIG. 5 .

 5 FIG. 5. Phonon NEGF simulation of ballistic heat transport through Si/Ge interface: (a) schematic of numerical model; blue (green) spheres represent Si (Ge); (b) frequency-dependent interfacial transmission; (c) frequency-dependent interfacial transmissivity (interfacial transmission divided by the transmission of pure Si). The circles and diamonds represent the results of ballistic phonon NEGF with DFT input from Refs. [54,57], respectively, whereas the solid line represents the result of the present work.

FIG. 6 .

 6 FIG.6. Schematic of the numerical model for phonon NEGF simulation of anharmonic heat transport across a silicon thin film with a thickness d. The cross-section period is a conventional unit cell of silicon with eight atoms, and the total number of atoms in the device region is eight times the number of unit cells along the transport direction. Here both the cross section and its two transverse periodic images are shown.

FIG. 7 .

 7 FIG. 7. Heat flow across a silicon thin film with a thickness of 20 uc (10.8 nm) at 300 K by phonon NEGF: (a) position-dependent spectral heat flow by anharmonic NEGF with DFT input; (b) position-dependent heat flow by ballistic NEGF with DFT input (dashed line) and by anharmonic NEGF with DFT input (solid line with circles).

FIG. 9 .

 9 FIG.9. Thickness-dependent nondimensional cross-plane thermal conductivity of silicon thin film at 300 K: experimental data (filled squares with error bar)[START_REF] Jiang | Role of low-energy phonons with mean-free-paths >0.8 μm in heat conduction in silicon[END_REF], previous Monte Carlo (MC) with DFT input (squares)[START_REF] Vermeersch | Cross-plane heat conduction in thin films with ab-initio phonon dispersions and scattering rates[END_REF], nonequilibrium molecular dynamics (NEMD) and Boltzmann transport equation (BTE) with the same Tersoff potential input (triangles)[START_REF] Hu | Unification of nonequilibrium molecular dynamics and the mode-resolved phonon Boltzmann equation for thermal transport simulations[END_REF], present MC with DFT input considering third (circles) and first (diamonds) nearest-neighboring shell in the third-order FC, the present anharmonic NEGF with DFT input, considering only first nearest-neighboring shell in the thirdorder FC (triangles with line).

FIG. 10 .

 10 FIG. 10. Temperature distribution in cross-plane heat transport through a silicon thin film with different thicknesses under a contact temperature difference of 10 K around 300 K: (a) d = 20 uc (10.8 nm); (b) d = 24 uc (13 nm). The squares represent the result by the present anharmonic phonon NEGF with DFT input, the solid line represents the result by the present Monte Carlo (MC) with the same DFT input considering only the first nearest-neighboring shell in the third-order FC, whereas the dashed line represents the uniform distribution in purely ballistic limit.

FIG. 11 .

 11 FIG. 11. Spectral heat flow in cross-plane heat transport through a silicon thin film with a thickness of 20 uc (10.8 nm) at 300 K: The dashed line denotes the result by ballistic phonon NEGF, the solid line represents the result by anharmonic phonon NEGF. The spectral heat flow has been averaged over the 20 slabs of the device for both cases.

FIG. 12 .

 12 FIG. 12. Temperature-dependent cross-plane heat transport through silicon thin film with a thickness of 20 uc (10.8 nm): (a) thermal conductivity and (b) dimensionless temperature distribution by the present anharmonic NEGF and by the Monte Carlo (MC) with consistent DFT input considering only the first nearest-neighbor shell in the third-order FC; the dashed lines denote the ballistic solution of the Boltzmann transport equation (BTE) with the same DFT input. The dimensionless temperature is defined as (T -T c )/(T h -T c ), with T h and T c the temperatures of two contacts.

FIG. 13 .

 13 FIG. 13. Spectral heat flow in cross-plane heat transport through silicon thin film with different thicknesses at 300 K: (a) 3 uc, (b) 5 uc. The dashed line and solid line represent, respectively, the result of the present ballistic and anharmonic NEGF with DFT input; the dash-dotted line represents the result of the anharmonic NEGF considering only 3 × 3 diagonal blocks of the anharmonic phonon scattering self-energy matrix.

  FIG.[START_REF] Maire | Heat conduction tuning by wave nature of phonons[END_REF]. Connected Feynman diagrams for three-phonon anharmonic scattering process.

TABLE III .

 III Mesh independence verification for phonon NEGF simulation of anharmonic heat transport across a silicon thin film with a thickness d = 5 uc (1 uc = 5.4018 Å) at 300 K.

	Frequency mesh	Transverse wave vector mesh	Ballistic thermal conductance (MW/m 2 K)	Thermal conductance (MW/m 2 K)
	61	4 × 4	1130.18	937.73
	81	4 × 4	1090.59	928.05
	81	6 × 6	1069.35	894.55
	101	4 × 4	1074.60	894.32
	101	6 × 6	1040.03	890.97
	101	8 × 8	1053.53	872.05
	121	6 × 6	1061.24	883.33
	121	8 × 8	1065.81	893.08

TABLE IV .

 IV The cross-plane thermal conductivity of silicon thin film at 300 K predicted by Monte Carlo (MC) and by the present anharmonic phonon NEGF with consistent DFT input.

	Film thickness d	Thermal conductivity by MC + DFT (W/m K)	Thermal conductivity by the present NEGF + DFT (W/m K)
	20 uc (10.8 nm)	8.4640	7.4489
	24 uc (13 nm)	9.7662	8.4025
	by 26% compared to the present result (689.45 MW/m 2 K)
	for the simulated case. In contrast, the scattering self-energy
	in Ref. [64] underestimates the thermal conductance (744.76
	and 642.16 MW/m 2 K) by 21% and 28% separately com-
	pared to the present results (941.72 and 890.97 MW/m 2 K)

ACKNOWLEDGMENTS Y.G. appreciates helpful discussions with S. Fiore from ETH Zürich on the reconstruction of the harmonic FC matrix from DFT calculation. This work was supported by the Postdoctoral Fellowship of Japan Society for the Promotion of Science (Grant No. P19353), and CREST Japan Science and Technology Agency (Grants No. JPMJCR19I1 and No. JP- MJCR19Q3). This research used the computational resources of the Oakforest-PACS supercomputer system, The University of Tokyo.

APPENDIX A: DERIVATION OF ANHARMONIC SCATTERING SELF-ENERGY THROUGH DIAGRAMMATIC PERTURBATION EXPANSION APPENDIX C: PROOF OF A GENERAL RELATION BETWEEN NONEQUILIBRIUM PHONON GREEN'S FUNCTIONS FOR 3D NANOSTRUCTURES

Based on the definition of Fourier's representation in Eq. [START_REF] Maldovan | Sound and heat revolutions in phononics[END_REF], the greater/lesser phonon Green's functions are expressed, respectively, as

The element of the transpose of the greater Green's function matrix is

with R ⊥ = (l yl y )a 2 + (l zl z )a 3 here. The following symmetry relation is valid between the greater/lesser phonon Green's function for 1D nanostructures without any periodicity [START_REF] Wang | Nonequilibrium Green's function method for thermal transport in junctions[END_REF]:

Substituting Eq. (C4) into Eq. (C3), we obtain

Since we have R ⊥ = -R ⊥ , Eq. (C5) can be rewritten into

Based on the definition in Eq. (C2), Eq. (C6) becomes exactly

which can be reformulated into matrix form as

APPENDIX D: MPI SCHEME FOR PARALLELIZED CALCULATION OF ANHARMONIC PHONON SCATTERING SELF-ENERGY

We follow the basic procedure in the MPI scheme for the parallelized calculation of the electron-phonon scattering self-energy matrix in electron NEGF [START_REF] Luisier | Proceedings of the 2010 ACM/IEEE International Conference for High Performance Computing, Networking, Storage and Analysis[END_REF]. However, the present situation is slightly more complicated since for each mode (ω; q ⊥ ) we have to consider the Green's functions of two other modes, (ω ; q ⊥ ) and (ω -ω ; q ⊥q ⊥ ), due to the three-phonon anharmonic scattering process shown in Fig. 2. As a first step, for each mode (ω; q ⊥ ) in a CPU, we build a connection table which stores the information of all the possible connected modes (ω -ω ; q ⊥q ⊥ ) via all possible modes (ω ; q ⊥ ). The data exchange is then conducted based on the connection table. Since each mode (ω; q ⊥ ) requires the data of mode (ω ; q ⊥ ) and mode (ω -ω ; q ⊥q ⊥ ), the The solid line with circles represents the present anharmonic phonon NEGF simulation of heat transport across a silicon thin film with a thickness of 2 uc at 300 K; a frequency mesh of N m = 81 and a tranverse wave vector mesh of 4 × 4 is adopted for the test. The dashed line represents the ideal scaling.

data exchange consists of a sending substep and a receiving substep successively: (I) for all the CPUs, send the local data of mode (ω; q ⊥ ) to the CPUs corresponding to all the possible modes (ω ; q ⊥ ), and also to the CPUs corresponding to all the possible modes (ω + ω ; q ⊥ + q ⊥ ); (II) for all the CPUs, receive the data of all the possible modes (ω ; q ⊥ ) from the corresponding CPUs, and also the data of all the possible modes (ω -ω ; q ⊥q ⊥ ) from the corresponding CPUs.

Another important issue is the storage of a third-order FC matrix in the Fourier's representation ˜ i j 1 j 2 l x l 1x l 2x (q ⊥ , q ⊥ ), which will be an extremely large matrix for long nanostructures due to its dependence on two transverse wave vectors. We reduce the memory cost by storing its dependence on only the first wave vector in a local CPU attributed to the full parallelization of the transverse wave vector (the local wave vector acting as the second wave vector). In this way, we also need to exchange data of the third-order FC matrix when computing the anharmonic scattering self-energy Eq. ( 14). The idea and procedure are very similar to that for the data exchange of the Green's function and are not repeated here for elegance. Once the data exchange of both the Green's function and the third-order FC matrix is accomplished, the scattering self-energy for each mode (ω; q ⊥ ) can then be calculated based on Eq. ( 14).

The scalability of the present MPI parallelization scheme is demonstrated in Fig. 16, which shows a scaling of computational time cost versus number of CPUs very close to the ideal scaling limit.

APPENDIX E: RECONSTRUCTION OF HARMONIC AND THIRD-ORDER FORCE CONSTANT MATRICES FOR NANOSTRUCTURES

The basic idea and procedure of the reconstruction of the harmonic FC matrix is shown in Fig. 17. Firstly, all the