{'fraction_non_alphanumeric': 0.03641758834637173, 'fraction_numerical': 0.01807391421634745, 'mean_word_length': 5.044767269493033, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum mechanical modeling of anharmonic phonon-phonon scattering in nanostructures

The coherent quantum effect becomes increasingly important in the heat dissipation bottleneck of semiconductor nanoelectronics with the characteristic size shrinking down to few nano-meters scale nowadays. 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 non-equilibrium Green function (NEGF) formalism for both 1D and 3D nanostructures, through a diagrammatic perturbation expansion and an introduction of Fourier 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 silicon thin film. The phonon-phonon scattering is shown to be appreciable and to introduce about 20% reduction of thermal conductivity at room temperature even for a film thickness around 10 nm. The present methodology provides a robust platform for the device quantum thermal modeling, as well as the study on the transition from coherent to incoherent heat transport in nano-phononic crystals. This work thus paves the way to understand and to manipulate heat conduction via the wave nature of phonons.


Introduction
The manipulation and control of heat transport in 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 [1] and thermal barrier coating materials [2], whereas high phonon thermal conductivity is desired for heat dissipation problems in nanoelectronics [3]. 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 [4][5][6][7]. These efforts are mainly based on the particle nature of phonons [8,9], since the characteristic size of the system is usually much larger than the dominant width (or called coherent length) of phonon wave-packets. With the development of nano-fabrication and manufacturing technology in recent years, heat conduction tuning via the wave nature of phonons based on nano-phononic crystals (PnC) becomes also possible [10][11][12][13]. Distinctive behaviors of phonon heat transport have been experimentally demonstrated in this wave regime: the period-controlled lattice thermal conductivity minimum of superlattices (SLs) [14], the order of nano-PnC slowing down the heat propagation [15], nanodot induced Anderson localization of thermal phonons in SLs [16], the aperiodicity of SLs reducing thermal conductivity via multiple localization [17], and so on [18][19][20].
The modeling of phonon heat transport in this wave-like coherent regime to particle-like incoherent regime remains, however, a challenging task. The widely adopted semi-classical phonon Boltzmann equation becomes no longer valid [21,22] as it describes only the evolution of phonon-particle population whereas ignoring 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 [15,23,24], 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 [25][26][27][28]. As MD is a classical method, some 3 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.
Non-equilibrium Green's function (NEGF) method [29][30][31][32] 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, NEGF method describes the evolution of both population and coherences of phonon modes [33] such that it will be able to capture the wave-like to particle-like behaviors of phonon heat transport. The NEGF formalism was originally proposed around the last mid-century [34][35][36] during the development of quantum field theory [37] mainly for electrons, and has been relatively well established and widely applied in quantum transport modeling of nano-electronic devices [38][39][40][41]. Later the NEGF method was introduced to model phonon transport in the early years of this century attributed to several pioneering works [42][43][44][45]. Due to the large computational cost, the phonon NEGF has been often applied for ballistic heat transport through relatively simple structures like lowdimensional nanostructures (ultra-thin nanowire [42,46,47], nanotubes [45,48], molecular junctions [49], 2D structures [50][51][52], etc.), interfaces [53][54][55][56][57][58] and SLs [17,[59][60][61]. Very few works [43,44,[62][63][64][65][66] directly take into account the anharmonic phonon-phonon scattering, yet often considering few-atom systems like atomic chain or junction [43,44,62,63], and single-unit-cell interface [66]. An indirect treatment of incoherent phonon scattering was also proposed by the Büttiker probe approach [67][68][69] 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 Section 2, the crucial selfenergy expressions for phonon-phonon scattering remains still diverse in several prototypical literature [43,44,64]. As a result, there still lacks a solid quantitative 4 validation of the anharmonic phonon NEGF formalism to the authors' best knowledge.
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 first-principle-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 [64] by one of the co-authors will be elucidated to ensure feasible yet still accurate large-scale simulation. The remaining of this article will be arranged as follows: the theoretical formulations and numerical implementation of the phonon NEGF method will be described in Section 2; the validation results and related discussions will be shown in Section 3; concluding remarks are made in Section 4.

Mathematical and Numerical Models
In this section, we will present a clarified anharmonic phonon NEGF formalism for each relevant physical system. Section 2.1 will firstly focus on 1D nanostructures without periodicity. Then a Fourier's representation is introduced to the anharmonic formalism in 5 Section 2.2 for 3D nanostructures with transverse periodicity, which is a paradigm for quantum heat transport across thin films, interfaces, SLs and multi-layer nanostructures. In Section 2.3, the numerical implementation of the phonon NEGF formalism will be elucidated in detail, including the computational strategies and techniques for large-scale simulation. In Section 2.4, the input of harmonic and anharmonic FC matrices from firstprinciple (DFT) calculation will be introduced.

Anharmonic phonon NEGF formalism for 1D nanostructures
In this sub-section, we firstly consider quantum heat transport through a 1D nanostructure without any translational periodicity, as it is shown in Figure 1(a). The total Hamiltonian operator of the system in the Heisenberg representation is expressed as [8]: where H0(t) and V(t) are the exactly solvable harmonic part and the third-order anharmonic perturbation part respectively: where n, m, l denotes atomic index, whereas i, j, k denotes Cartesian coordinate index (x, y, z where E denotes the atomic interaction potential, and the subscript '0' denotes equilibrium position of crystal lattice.
The governing equation for the retarded Green's function of the device in a matrix notation is [43,64]: where I denotes the unity matrix, Φ denotes the second-order FC matrix with its component given in Eq. (4), 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 one. ω represents the frequency of the phonon modes through the system. The retarded self-energy matrix includes the contribution from both the two contacts and the anharmonic phonon-phonon scattering within the device: The retarded contact self-energy matrix in Eq. (7) is related to the uncoupled retarded Green's function of the contacts ( ( ) R 1(2)  g ) as [43,53]: with τ1 (or τ2) denoting the harmonic interaction FC matrix between the device and contacts. For brevity, the subscript '1(2)' is introduced to refer to contact 1 or contact 2 hereafter. The contact self-energy matrix is related to the surface Green's function of the contacts, which are calculated by the decimation technique [32,71]. 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. (6) to get the transmission through the system using 7 the Caroli formula [29,42,45]. However, for anharmonic heat transport, we also have to compute the greater/lesser Green's function as [43,64]: where the advanced Green's function of the device is the Hermitian conjugate of the retarded Green's function defined as: GG. The greater/lesser self-energy matrix also includes the contribution from both the two contacts and the anharmonic phonon-phonon scattering in the device: The uncoupled greater/lesser Green's function of the contacts could be derived from the fluctuation-dissipation theorem at equilibrium [62] and are computed as: , fT in Eq. (10). The retarded scattering selfenergy matrix in Eq. (7) is then related to the greater/lesser scattering self-energy matrix as [64]: The second term on the right-hand side of Eq. (14) is the Cauchy principal integral, and is often neglected for computational simplicity [64]. 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 8 similar expressions as available in the literature yet with different coefficients [43,44,64], as it is summarized in Table 1. For a clear comparison, we keep the notation of third-order FC in the original literature whereas 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 [43,44,64]. 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. Table 1. Summary of expressions for greater/lesser anharmonic phonon scattering self-energy.
Third-order term in Hamiltonian Scattering self-energy References   9 with the detailed derivation process given in Appendix A. In comparison to the expressions of anharmonic scattering self-energy in Ref. [43] and Ref. [64], there exists a factor difference of π/9 and 4 respectively with respect to the present result of Eq. (15).
Considering that the third-order FC ( ijk T in Table 1) is the half of the present one ( ijk nml  ), the part of self-energy corresponding to the three-phonon anharmonic scattering in Ref. [44] is actually consistent with the present expression. A further theoretical corroboration of Eq.
(15) will be shown in Section 2.2 when we compare it to the expression of anharmonic phonon scattering self-energy in bulk materials. Up to now, Eq. (6), Eq. (9), Eq. (14) and Eq. (15) constitute a closed set of coupled matrix equations, which have to be solved through an iterative process called self-consistent Born approximation (SCBA) [41,64].

Anharmonic phonon NEGF formalism for 3D nanostructures
In this sub-section, we consider quantum heat transport through 3D nanostructures with transverse periodicity as it is shown in Figure 1(b). Attributed to this periodicity in the transverse direction, the phonon NEGF formalism in Section 2.1 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: where N is the number of transverse wave vectors, and ( ) ( ) and Eq. (15) become: ( , ) ; 2 jj j ij j ij l l l l l l l l l l l l j j j j In Eq. (18), the Fourier's representation of the harmonic FC matrix is defined as: which is consistent with the definition in the previous ballistic NEGF formalism [53,58]. In Eq. (21), the Fourier's representation of the third-order anharmonic FC matrix is defined as: ( , ) ( , ) ; ; where we have kept most of the notations therein. In spite of a similar mathematical form, our result of Eq. (21) is different from Eq. (24) 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. [43], as is already shown in Table 1  The expression of anharmonic phonon-phonon scattering self-energy in Eq. (21) 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 similar diagrammatic perturbation expansion of the phonon Green's function as that in Appendix A, one could obtain the self-energy for three-phonon anharmonic scattering in bulk material as follows [72]: where q is the 3D wave vector in bulk material, ( ) ( ) ( ) ,, F q q q is the third-order anharmonic dynamic matrix contributing to the anharmonic term in the Hamiltonian as [72]: is the phonon normal coordinate operator, with † , aa qq the phonon destruction and creation operator respectively. Due to the translational symmetry of the 3D lattice, ( ,, F q q q includes a delta function as ( ) q q q which will be nonvanishing only when 1 2 3 0 + + = q q q or a reciprocal lattice vector [72,73]. Thus Eq. (25) can be rewritten into: Considering the different notations in the third-order term of the Hamiltonian between Eq.
Finally, some general symmetry relations between the non-equilibrium phonon Green's functions for 3D nanostructures with transvese peridocity are demonstrated: The subscript 'T' in Eq. (28) denotes the transpose of a matrix. We provide a detailed proof of Eq. (28) in Appendix C. The other relations could 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 [62] is summarized in Table 2. Table 2. General relations between the non-equilibrium phonon Green's functions for 1D and 3D nanostructures.

Numerical implementation
The numerical implementation of the anharmonic phonon NEGF formalism by the SCBA iterative solution of Eqs.  and 'RC' denote left contact (contact 1) and right contact (contact 2) respectively, and the device is divided into Nd slabs.

Recursive algorithm
The recursive algorithm was originally proposed in electron NEGF codes for largescale nano-electronic device modeling [74], and has been also applied in anharmonic phonon NEGF modeling of heat transport in nanowires [64]. Similar variant of recursive algorithm was also developed in ballistic phonon NEGF [55]. 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 slab such that the memory cost and time cost will be linearly proportional to the system size along the transport direction as (Nx) 2 Nd and (Nx) 3 Nd respectively [39]. In this work, we adopt the recursive algorithm [39,74] to compute the 14 Nx×Nx 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 Figure 3 is considered for harmonic FC; (II) only the atomic nearest-neighbor interaction is considered for third-order anharmonic FC; (III) only the Nx×Nx 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 [64] 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 it will be shown later in Section 2.4. In principle, more neighbors should be considered in the third-order FC for strongly anharmonic materials like oxides, yet it remains computationally very challenging for large-scale simulations. Concerning approximation (III), an even simpler approximation is assumed in Ref. [64] where only the 3×3 diagonal blocks of anharmonic scattering selfenergy matrix are computed. We will show in Section 3.3 that the approximation in Ref.
[64] 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 [64].

Parallelization scheme
With multiple-CPU (central processing unit) computer facility, the iterative solution of Eqs. (18)-(21) can be parallelized based on the message-passing-interface (MPI) standard [75]. The MPI parallelization scheme has been widely adopted in electron NEGF for large-scale device simulations [41,76], and is also used in anharmonic phonon NEGF in a previous work [64] by one of the co-authors. In comparison to the parallelization scheme in Ref. [64] for nanowire (1D nanostructures) simulations which require data exchange between different frequency points as shown in Eq. (15), here we also have to exchange the data between different transverse wave vectors for simulation of 3D nanostructures as 15 shown in Eq. (21). Thus, we build a dual-level parallelization scheme for the transverse wave vector and frequency points, both of which are uniformly discretized as it is shown in

Macroscopic variable calculation
The local heat flow from the q th slab to (q+1) th slab shown in Figure 3 is related to the first off-diagonal blocks of the lesser phonon Green's function [64] and can be computed as: Eq. (32) can be written as: with the heat flow contributed by each mode expressed as: 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. (32). 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: where Rn denotes the spatial position of the q th slab, and the diagonal blocks of the lesser phonon Green's function [77] and is computed as: with 'Tr' denoting the trace of a square matrix, and the subscript 'nn' represents the subblock of G < corresponding to the q th slab. The local density of states (LDOS) in Eq. (34) In a previous work on coupled electron-phonon NEGF modeling [77], a different definition of local temperature was adopted based on the local phonon number conservation condition.
This approach amounts to removing one  in the integration on both sides of Eq. (34).
Here we keep the convention of using local energy density to characterize the local temperature in non-equilibrium transport as currently done in the heat transport community [21,[67][68][69].  19 where the '−' and '+' signs respectively correspond to the first and second index of FC matrix in the summation. The schematic procedure of the reconstruction of third-order FC matrix is shown in Figure 6. 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 12 12 ,0 ( , )

Results and Discussions
In this section, we mainly aim to demonstrate the validity of the phonon NEGF formalism and computational framework introduced in Section 2. Firstly, a simple case of ballistic heat transport through a Si/Ge interface is considered in Section 3.1. The quantitative validation of the anharmonic phonon NEGF formalism remains challenging due to its limited computational capability to extremely small structures and/or the usual empirical anharmonic FC input in previous works. In Section 3.2, we tackle this challenge by modeling anharmonic heat transport across a silicon thin film with a thickness larger than 10 nm through our DFT-based large-scale NEGF simulation. Finally, we discuss some crucial approximations in the treatment of anharmonic phonon scattering self-energy in Section 3.3.

Validation: ballistic heat transport
We consider ballistic heat transport through a Si/Ge interface as shown in Figure   7(a), with only one unit cell in the device region. The lattice constant and harmonic FC of 22 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 [54,57]. This approximate treatment is rather reasonable since the lattice structure of Ge is indeed very similar to that of Si.
After numerical convergence test, a frequency mesh of Nm=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 homogenous Si film to be discussed in sub-section 3.2. This is due to the very different cutoff frequencies (almost two-fold difference) in the phonon bandstructures of Si and Ge. The ballistic transmission is calculated based on the Caroli formula [29,54]. The spectral transmission and transmissivity through the Si/Ge interface are shown in Figure 7(b) and (c) respectively, which demonstrates an overall good agreement with respect to previous studies [54,57]. The minor difference may arise from the slightly different DFT harmonic FCs calculated by different groups. The circles and diamonds represent the results of ballistic phonon NEGF with DFT input from Ref. [54] and Ref. [57] respectively, whereas the solid line represents the result of present work.

Validation: anharmonic heat transport
In this sub-section, 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 silicon thin film around 300K is considered, as shown in Figure 8. The temperatures of left contact (contact 1) and right contact (contact 2) are set at 305K and 295K respectively (with a temperature difference 10K T = ). The heat transport under such a temperature difference lies still within a linear regime based on our numerical test under a tiny temperature difference of 0.1K. A frequency mesh of 101 and transverse wave vector mesh 6×6 is adopted after numerical convergence test as it is summarized in Table 3.
Denser meshes produce thermal conductance within about 2% variation comparing to the present mesh. We consider a series of thickness for the silicon thin film up to 24uc (13 nm)  is a very good model for silicon with weak normal scattering [9,80]. As a validation of the present MC code, the thickness-dependent cross-plane thermal conductivity of silicon thin film is shown in Figure 10(a) when the 3 rd nearest neighboring shell is considered in the third-order DFT FC. The MC result shows a very good agreement with i) previous MC study considering DFT input [84]; ii) a recent experimental measurement [85].
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 silicon thin film with a thickness of 20uc (10.8 nm) and 24uc (13 nm) is shown in Table 4. It is seen that the result predicted by the present NEGF simulation is very close to the corresponding MC result. The about 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 Section 2.3.1. Although these terms are small, including more terms in the scattering self-energy is expected to increase the thermal conductivity a bit, as it will be inferred later in Section 3.3. Taking into account such effect, this comparison can act as a quantitative validation of the present anharmonic phonon NEGF formalism and numerical framework.  Figure 9. Heat flow across silicon thin film with a thickness of 20uc (10. It is also interesting to note that the results of the phonon Boltzmann equation and the present phonon NEGF formalism agree also well for a film thickness smaller than 10 nm, as shown in Figure 10(b). This indicates that the phonon Boltzmann equation seems to still work for cross-plane heat transport through a silicon thin film of few nanometers. The underlying physical interpretation remains to be investigated in the near future. One possible explanation might be that the coherence of phonons still contributes negligibly to heat transport in this situation. Quantifying the coherence of phonons and its contribution to heat transport remains a 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 Figure 11. At small thickness, the present anharmonic phonon NEGF result well agrees with the general trend of non-dimensional thermal conductivity from different reference data.
We also provide the results of local temperature distribution across the silicon thin film with a thickness of 20uc (10.8 nm) and 24uc (13 nm) in Figure 12. The results by the present anharmonic phonon NEGF simulation show a good agreement with those of the 27 MC simulation with the same DFT input. A large temperature jump near the two contacts is obtained due to the strong non-equilibrium effect between the thin film and contacts, which is well known in the heat transport community and has been already reported [6,7]. 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 Figure 10(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 Figure 13.  Figure 11. Thickness-dependent non-dimensional cross-plane thermal conductivity of silicon thin film at 300K: experimental data (filled squares with error bar) [85], previous Monte Carlo (MC) with DFT input (squares) [84], non-equilibrium molecular dynamics (NEMD) and Boltzmann transport equation (BTE) with the same Tersoff potential input (triangles) [82], present MC with DFT input considering 3 rd (circles) and 1 st (diamonds) nearest neighboring shell in the third-order FC, the present anharmonic NEGF with DFT input, considering only 1 st nearest neighboring shell in the third-order FC (triangles with line).
(a) (b) Figure 12. Local temperature distribution in cross-plane heat transport through silicon thin film with different thicknesses under a contact temperature difference of 10K around 300K: (a) d=20uc (10.8 nm); (b) d=24uc (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 1 st nearest neighboring shell in the third-order FC, whereas the dashed line represents the uniform distribution in purely ballistic limit. 29 Figure 13. Spectral heat flow in cross-plane heat transport through silicon thin film with a thickness of 20uc (10.8 nm) at 300K: 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.

Assessment of approximations in numerical implementation
Finally, we briefly 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 Section 2.3.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 co-authors [64]. 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 Nx×Nx diagonal blocks (c.f. Section 2.3.1). The heat transport across silicon thin film with two thicknesses of 3uc and 5uc at room temperature is simulated. For the thickness of 3uc, 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 5uc, 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 3uc and 5uc respectively. The underestimation will increase with increasing film thickness due to the stronger phonon- 30 phonon scattering rate. The underlying reason can be understood from the spectral heat flow results shown in Figure 14. 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 retrieving the underestimation of thermal conductivity shown in Table 4 and Figure 10(b). However, since the dominant terms have been already considered, the present treatment constitutes a good balance between the computational accuracy and efficiency. Note that a further approximation was made in the product of greater/lesser Green's function matrix in the phonon scattering self-energy of Eq.
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 [41,77] 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 phonon Green's function, which is a more complicated non-linear 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 [64]. 31 (a) (b) Figure 14. Spectral heat flow in cross-plane heat transport through silicon thin film with different thicknesses at 300K: (a) 3uc, (b) 5uc. 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.

Conclusions
In summary, we present a non-equilibrium 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 Boltzmann equation in the particle transport regime based on a classical cross-plane heat transport in 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 scattering rates. 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 32 incoherent heat transport in nanophononic crystals like superlattices, which is pending in our future work. This study thus opens new opportunities for the understanding and tuning of phonon heat transport in the coherent quantum wave regime.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request. 33

Appendix A. Derivation of anharmonic scattering self-energy through diagrammatic perturbation expansion
The starting point of the theoretical derivation is the definition of contour-ordered phonon Green's function for heat transport in 1D nanostructures [43,44]: where < >0 denotes the expectation over the exactly solvable harmonic part at equilibrium with the density matrix defined as [37,40]: where the third-order perturbation part of the Hamiltonian in the interaction representation is: The evolution operator is the basis for the diagrammatic perturbation expansion.
Considering the infinite series in Eq. (A3) within second-order and substituting into Eq.
(A2), we obtain: , , , , where the zeroth-, first-and second-order terms are respectively: The first-order term in Eq. (A7) will be vanishing since it involves odd number of atomic displacement operators. The full expression of the second-order term will be: (2) 12 The Wick's theorem [37,86]  (2) 12 which can be rewritten into:  35 Putting Eq. (A6) and Eq. (A11) into Eq. (A5), we obtain the contour-order phonon Green's function within second-order as: Comparing Eq. (A12) to the Dyson's equation it follows [37,40]: we get the expression of anharmonic phonon scattering self-energy as: where the unperturbed Green's function has been replaced by the full Green's function based on the self-consistent Born approximation.
For the stationary state heat transport considered in this work, the dependence on t1 and t2 will be reduced to the dependence on (t1-t2). For convenience, the following Fourier transform and its inverse transform is introduced [86]: which is exactly Eq. (15) in the main text. It is also seen from Eq. (A18) 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. (15) for 1D nanostructures: where all the subscripts ' ⊥ ' are omitted to keep a simple notation, and the relative displacements are defined as below: ( ) ( ) where a classical relation in lattice dynamics has been used [8,73]: The right-hand side term becomes:    39 In the derivation of Eq. (B19), we have identified the sum over ( ) , yz ll as a sum over 3 R , and also used the relation Eq. (B18). Therefore, the scattering self-energy is:

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 electron-phonon scattering self-energy matrix in electron NEGF [76]. However, the present situation is slightly more complicated since for each mode ( ) Another important issue is the storage of third-order FC matrix in the Fourier's representation 12 12 ( , ) x x x ij j l l l ⊥⊥   qq , 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 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 The scalability of the present MPI parallelization scheme is demonsrated in Fig. D1, which shows a scaling of computational time cost versus number of CPUs very close to the ideal scaling limit.