*2.1. Ballistic Phonon Transport*

One of the most powerful methodologies to study quantum (thermal) transport is the non-equilibrium Green's functions (NEGF) formalism [35,64,65]. The NEGF method has its origin in quantum field theory [66], and has been developed to study many-particle quantum systems under both equilibrium and nonequilibrium conditions. Different formulations were derived during the early 1960s [67–69]. Thus, Keldysh developed a diagrammatic approach, Kadanoff and Baym formulated their approach based on equations of motion. Both methods are suitable for studying a dynamic system in nonequilibrium. For instance, by using the Keldysh formalism, one can obtain formal expressions for the current and electron density [70]. The method has also been successfully used to study electron transport properties in open quantum systems [71,72]. Moreover, NEGF has been recently used on thermal transport investigations not only in the ballistic regime [73–75], but also including phonon–phonon scattering [76–79]. In this section, we describe the NEGF formalism to address ballistic phonon transport in nanostructures. Phonon–phonon interactions require the inclusion of extra self-energy terms, which depend on products of single-phonon Green's functions, so that the whole problem must be solved self-consistently; however, this goes beyond the scope of this review (see [35,64] for more details concerning this point).

The main difference between the NEGF formalism and ordinary equilibrium theory is that all time-dependent functions are defined on the so-called Schwinger-Keldysh contour (see Figure 1). However, a simplification occurs when *t*0 → − ∞ (Keldysh contour). If the interactions are switch on adiabatically, the contribution from the [*<sup>t</sup>*0, *t*0 + *β*] piece vanishes. The information lost by this procedure is related to initial correlations. In many physical situations, such as in steady state transport, it is a plausible approximation to assume that initial correlations have been washed out by the interactions when one reaches the steady state. On the contrary, for the transient response, the role of initial correlations may be important (see Section 3). Here, we consider the Keldysh contour consisting of two branches running from −∞ to ∞ and from ∞ to <sup>−</sup>∞. Therefore, one can introduce a contour-ordered Green's function as [80]:

$$G(\mathbf{r}, \mathbf{r}') = -i \left\langle T\_{\mathbb{C}} A(\mathbf{r}) B^T(\mathbf{r}') \right\rangle,\tag{1}$$

with *TC* as the contour-order operator. Based on it, six real-time Green's functions can be defined [35]:

**Figure 1.** Schwinger/Keldysh contour *C* in the imaginary time plane, *C* = {*t* ∈ C, *t* ∈ [*<sup>t</sup>*0, ∞] *t* ∈ [*<sup>t</sup>*0, <sup>−</sup>*β*]}. For more clarity, the different contour branches are displayed slightly off the axes. Time-ordering: time *t*2 is later on the contour than time *t*, and *t* is larger than *t*1 [80].


*A*(*t*) and *B*(*t*) are operators in the Heisenberg picture and Θ(*t*) is the Heaviside step function. The angular brackets denote trace with the canonical density matrix, i.e., ··· = Tr(*ρ* ···), with *ρ* = *e* <sup>−</sup>*β<sup>H</sup>*/Tr(*e*<sup>−</sup>*β<sup>H</sup>*) and *β* = 1/(*kBT*), and *H* is the Hamiltonian of the system. The notation [*<sup>A</sup>*, *B<sup>T</sup>*] represents a matrix and should be understood as *AB<sup>T</sup>* − *BA<sup>T</sup><sup>T</sup>*.

In equilibrium or non-equilibrium steady state, the Green's functions only depend on the time difference, *t* − *t*. The Fourier transform of *Gr*(*t* − *t*) = *Gr*(*<sup>t</sup>*, *t*) is defined as *<sup>G</sup><sup>r</sup>*[*ω*] = - +∞ −∞ *<sup>G</sup><sup>r</sup>*(*t*)*eiωtdt*. Using the basic definitions, the following linear relations hold in both frequency and time domains [35]:

$$\begin{aligned} G^r - G^a &= G^{>} + G^{<}, \\ G^t + G^{\bar{l}} &= G^{>} + G^{<}, \\ G^t - G^{\bar{l}} &= G^r + G^a. \end{aligned} \tag{2}$$

Only three of the six Green's functions are linearly independent. In systems with time translational invariance, the functions *Gr* and *Ga* are related by *Ga*[*ω*] = (*Gr*[*ω*])†. Hence, under general non-equilibrium steady-state conditions, only two are independent, with a typical choice of working with *Gr* and *G*<sup>&</sup>lt;, although other combinations are possible. Extra relations are defined in the frequency domain for bosons [81]:

$$\begin{aligned} G^{<} [\omega]^{\dagger} &= -G^{<} [\omega], \\ G^{r} [-\omega] &= G^{r} [\omega]^{\*}, \\ G^{<} [-\omega] &= G^{>} [\omega]^{T} = -G^{\left[\omega\right]} ^{\*} + G^{r} [\omega]^{T} - G^{r} [\omega]^{\*}. \end{aligned} \tag{3}$$

Therefore, based on the last two equations, only the positive frequency part of the functions is needed. Equations (2) and (3) are generally valid for non-equilibrium steady states. However, for systems in thermal equilibrium satisfying the fluctuation–dissipation theorem [82], there is an additional equation relating *Gr* and *G*<sup>&</sup>lt;:

$$\mathcal{G}^{<}\left[\omega\right] = f(\omega) \left(\mathcal{G}^{r}[\omega] - \mathcal{G}^{a}[\omega]\right),\tag{4}$$

where *f*(*ω*) = *e h*¯ *ω kBT* − 1 −1 is the Bose–Einstein distribution function at temperature *T*. *kB* is the Boltzmann constant. Indeed, the correlation function *G*<sup>&</sup>lt; contains information of fluctuations, while *Gr* − *Ga* describes dissipation of the system. *<sup>G</sup>*<sup>&</sup>gt;[*ω*] = *eβω <sup>G</sup>*<sup>&</sup>lt;[*ω*] also applies for equilibrium systems and, consequently, there is only one independent Green's function under equilibrium conditions.

In phonon transport calculations, a partitioning scheme is applied, consisting in splitting the whole system in three regions: one central region (also denoted as device region), connected to two thermal baths on the left (L) and right (R) (see Figure 2). At the simplest level, the thermal baths can be considered as a collection of non-interacting harmonic oscillators. All elastic and/or inelastic scattering processes are therefore assumed to be confined to the central (or device) region. Since we focus on thermal transport mediated by the vibrational system, the phonon Hamiltonian of the whole system is given by [80]:

$$H = \sum\_{a=L,\mathcal{C},\mathcal{R}} H\_a + (u^L)^T V^{LC} u^C + (u^C)^T V^{CR} u^R + V\_{n\prime} \tag{5}$$

where *Hα* = 1 2 (*u*˙*<sup>α</sup>*)*Tu*˙*<sup>α</sup>* + 1 2 (*u<sup>α</sup>*)*TK<sup>α</sup>u<sup>α</sup>* represents the Hamiltonian of the region *α*; *α* = *L*, *C*, *R*, for the left, center, and right regions, respectively. *u<sup>α</sup>* is a column vector consisting of all the displacement variables in region *α*, and *u*˙*<sup>α</sup>* is the corresponding conjugate momentum. The following transformation of coordinates has been considered, *uj* = √ *mjxj*, where *xj* is the relative displacement of *j*th degree of freedom. *Kα* is the mass-reduced force constant matrix. This matrix is the mass-weighted second derivative of energy with respect to displacement at the equilibrium positions:

$$\left[\mathcal{K}\_{\gamma\beta}\right]\_{ij} = \mathcal{K}\_{ij}^{\gamma\beta} = \frac{1}{\sqrt{m\_i m\_j}} \frac{\partial^2 E}{\partial u\_i^{\gamma} u\_j^{\beta}}.\tag{6}$$

*VLC* = ( *VCL*)*<sup>T</sup>* is the coupling matrix between the left lead to the central region; and similarly for *VCR*. The last term *Vn* represents possible many-body interactions, such as phonon–phonon interaction [35]. It contains higher order (higher than 2) derivatives of the energy with respect to the displacements, evaluated at the equilibrium positions.

**Figure 2.** Schematic representation of the common partitioning scheme for phonon transport calculation using Green's function technique. The entire system is split into three regions: central region and, left and right heat baths. Each of this region are characterized by their own Hamiltonian *Hα* with *α* = *L*, *C*, *R*. The coupling matrices between heat baths and central region are *VLC* and *VRC*.

The most important quantity to calculate is the heat flux *J*, which is defined as the energy transferred from the heat source to the junction in a unit time, and is equal to the energy transferred from the junction to the heat sink in a unit time. Here, it is assumed that no energy is accumulated in the junction. According to this definition, the heat flux out of the left lead is:

$$J\_L = -\left< H\_L(t) \right> = i \left< [H\_L(t), H] \right> = i \left< \left[ H\_L(t), V^{LC}(t) \right] \right>. \tag{7}$$

In the steady state, energy conservation means that *JL* + *JR* = 0. For simplicity, we set *h*¯ = 1. Using Heisenberg's equation of motion, *JL* can be written as:

$$\begin{split} J\_L &= \left\langle (\dot{\mathbf{u}}^L)^T(t) V^{LC} \mathbf{u}^C(t) \right\rangle \\ &= \lim\_{\mathbf{t'} \to \mathbf{t}} \sum\_{j,k} V^{LC}\_{jk} \left\langle (\dot{\mathbf{u}}^L\_j)^T(t') \mathbf{u}^C\_k(t) \right\rangle. \end{split} \tag{8}$$

Thus, the heat flux depends on the expectation value of (*u*˙ *Lj* )*<sup>T</sup>*(*t*)*uCk* (*t*), which can be written in terms of the lesser Green's function *<sup>G</sup>*<sup>&</sup>lt;*CL*(*<sup>t</sup>*, *t*) = −*i u<sup>L</sup>*(*t*)*u<sup>C</sup>*(*t*)*<sup>T</sup><sup>T</sup>*. Since operators *u* and *u*˙ are related in Fourier space (frequency domain) as *u*˙[*ω*] = <sup>−</sup>*iωu*[*ω*], the derivative is eliminated and one obtains: 

$$J\_L = -\frac{1}{2\pi} \int\_{-\infty}^{\infty} \text{Tr}\left(V^{LC} G\_{\tilde{C}L}^{<}[\omega]\right) \omega \, d\omega \,. \tag{9}$$

The Green's functions of interacting systems can be efficiently obtained by solving their equations of motion (EOM) [64]. In this section, EOMs will only be used to obtain expressions for retarded and lesser GFs of the central region. This topic will be expanded with more details in Section 3. First, we have that the contour-ordered GF *<sup>G</sup>*(*<sup>τ</sup>*, *τ*) = −*i <sup>T</sup>τ<sup>u</sup>*(*τ*)*u*(*τ*)*<sup>T</sup>* satisfies the following equation:

$$-\frac{\partial^2 G(\mathbf{r}, \mathbf{r}')}{\partial \mathbf{r}^2} - \mathcal{K} G(\mathbf{r}, \mathbf{r}') = I \delta(\mathbf{r}, \mathbf{r}') \tag{10}$$

The equation per each region is obtained by partitioning the matrices *G* and *K* into submatrices *G<sup>α</sup>*,*α* and *K<sup>α</sup>*,*α* , *α*, *α* = *L*, *C*, *R*. The free Green's function for the decoupled system *g* is easily obtained by solving:

$$-\frac{\partial^2 g^a(\tau, \tau')}{\partial \tau^2} - \mathcal{K}^a g^a(\tau, \tau') = I \delta(\tau, \tau'). \tag{11}$$

The corresponding free GFs in frequency domain are written as:

$$\mathcal{g}\_a^r[\omega] = \left[ (\omega + i\eta)^2 - \mathcal{K}^a \right]^{-1} \, , \tag{12}$$

where *η* is an infinitesimal positive quantity to single out the correct path around the poles when performing an inverse Fourier transform, such that *gr* = 0 for *t* < 0. Other Green's functions can be obtained using the general relations among them (see Equations (2) and (3)). Hence, the contour-ordered non-equilibrium GF can be written as:

$$\mathcal{G}^{\rm CL}(\boldsymbol{\tau}, \boldsymbol{\tau}') = \int d\boldsymbol{\tau}'' G^{\rm CC}(\boldsymbol{\tau}, \boldsymbol{\tau}'') V^{\rm CL} \mathcal{g}^{\rm L}(\boldsymbol{\tau}'', \boldsymbol{\tau}'), \tag{13}$$

$$\mathcal{G}^{\rm CC}(\boldsymbol{\tau}, \boldsymbol{\tau}') = \mathcal{g}^{\rm C}(\boldsymbol{\tau}, \boldsymbol{\tau}') + \int d\tau\_1 \int d\tau\_2 \mathcal{g}^{\rm C}(\boldsymbol{\tau}, \tau\_1) \Sigma(\tau\_1, \tau\_2) \mathcal{G}^{\rm CC}(\tau\_2, \tau'), \tag{14}$$

with <sup>Σ</sup>(*<sup>τ</sup>*1, *<sup>τ</sup>*2) being the total self-energy including the coupling to the baths and given by:

$$
\Sigma(\tau\_1, \tau\_2) = \Sigma\_L(\tau\_1, \tau\_2) + \Sigma\_R(\tau\_1, \tau\_2) = V^{\mathbb{C}L} \mathcal{g}^L(\tau\_1, \tau\_2) V^{\mathbb{C}\mathbb{C}} + V^{\mathbb{C}R} \mathcal{g}^R(\tau\_1, \tau\_2) V^{\mathbb{R}\mathbb{C}}.\tag{15}
$$

*gL* and *gR* are the GF of the isolated semi-infinite leads. Since the bath-device interaction terms are short-ranged, it is usually only necessary to compute the projection of the bath GF on the layer directly in contact with the device. The resulting GF can be calculated by an iteration method [83] or by decimation techniques [84]. The Dyson equation (see Equation (14)) can be written in the frequency domain as:

$$\begin{split} \mathcal{G}^{r}\_{\mathcal{C}\mathcal{C}}[\omega] &= \left( (\omega + i\eta)^{2} I - \mathcal{K}^{\mathbb{C}} - \Sigma^{r}[\omega] \right)^{-1}, \\ \mathcal{G}^{<}\_{\mathcal{C}\mathcal{C}}[\omega] &= \mathcal{G}^{r}\_{\mathcal{C}}[\omega] \Sigma^{<}[\omega] \mathcal{G}^{a}\_{\mathcal{C}}[\omega]. \end{split} \tag{16}$$

Several physical quantities can be calculated using these relations, e.g., the local density of states (LDOS) is expressed as:

$$\eta\_i(\omega) = -\frac{2\omega}{\pi} \left( \text{Im} G'[\omega] \right)\_{\text{i\'i}} \,\, \, \, \, \, \tag{17}$$

and the total phonon DOS as *η*(*ω*) = ∑*Ni*=<sup>1</sup> *ηi*(*ω*). The DOS and LDOS give the distribution of phonons in frequency as well as in real space [81]. This is very useful to analyze quantum transport processes, as shown below. Although Equation (17) is obtained for the special case of no phonon–phonon interactions, the same formula is valid in the presence of phonon–phonon interactions described by an interaction self-energy such as in Equation (16). This typical approach assumes that the non-crossing approximation applies, allowing to treat the effect of contacts and interactions as two independent additive contributions. Clearly, this is valid in the limit of small interactions acting only within the central region.

Next, it is useful to introduce the Γ function describing the phonon scattering rate into the thermal baths:

$$
\Gamma[\omega] = i \left( \Sigma^\prime[\omega] - \Sigma^d[\omega] \right) = \Gamma\_\prime[\omega] + \Gamma\_\hbar[\omega], \tag{18}
$$

This function has an important relation with the spectral function, *<sup>A</sup>*[*ω*] = *Gr*[*ω*]Γ[*ω*]*Ga*[*ω*]. By applying the Langreth theorem [64] to Equation (13), the lesser GF *<sup>G</sup>*<sup>&</sup>lt;*CL*turns into:

$$\mathcal{G}\_{\mathcal{C}L}^{<\cdot<}[\omega] = \mathcal{G}\_{\mathcal{C}\mathcal{C}}^{r}[\omega]V^{\mathcal{C}L}g\_{L}^{<}[\omega] + \mathcal{G}\_{\mathcal{C}\mathcal{C}}^{<}[\omega]V^{\mathcal{C}L}g\_{L}^{a}[\omega]. \tag{19}$$

Consequently, the heat flux coming from the left lead (see Equation (9)) can be written as:

$$J\_L = -\frac{1}{2\pi} \int\_{-\infty}^{+\infty} d\omega \omega \text{Tr}\left( G^r[\omega] \Sigma\_L^{\le'}[\omega] + G^{<'}[\omega] \Sigma\_L^a[\omega] \right). \tag{20}$$

For simplicity, the subscripts *C* related to the central region have been dropped. The upperletters are used to identify Green's functions on the central region and lowercase letters for the leads. After symmetrizing with respect to the left and right leads, the heat flux becomes:

$$J = \frac{1}{4} \left( J\_L + J\_L^\* - J\_R - J\_R^\* \right) \,. \tag{21}$$

The final expression reads:

$$J = \int\_0^\infty \frac{d\omega}{2\pi} \hbar \omega \tau\_{\rm ph}[\omega] \left(f\_L - f\_{\rm R}\right). \tag{22}$$

This result is formally similar to the Landauer equation obtained for electron transport. Here, however, *fL*,*<sup>R</sup>* are the Bose–Einstein distributions for the left and right leads and *<sup>τ</sup>p<sup>h</sup>*[*ω*] is the phonon transmission function, given by:

$$\tau\_{ph}[\omega] = \text{Tr}\left(G^{\prime}[\omega]\Gamma\_L[\omega]G^{\prime}[\omega]\Gamma\_R[\omega]\right) \,. \tag{23}$$

The retarded GF of the central region connected to the thermal baths is given by:

$$G^r[\omega] = \left[ (\omega + i\eta)^2 I - K^\mathbb{C} - \Sigma^r\_L[\omega] - \Sigma^r\_R[\omega] \right]^{-1} \tag{24}$$

Then, the thermal conductance can then be computed according to *<sup>κ</sup>p<sup>h</sup>* = limΔ*T*→<sup>0</sup> *J* Δ *T* , Δ *T* as the temperature difference between the thermal baths, with *TL* = *T* + Δ *T*/2 and *TR* = *T* − Δ *T*/2, respectively. A linear expansion of the Bose–Einstein distribution in Δ *T* yields [80]:

$$\kappa\_{\rm ph} = \frac{1}{2\pi} \int\_0^\infty d\omega \,\omega \, T[\omega] \frac{\partial f(\omega)}{\partial T}. \tag{25}$$

Notice that the thermal conductance can only be obtained by an integration over the whole frequency range of the phonon transmission. In practice, the derivative of the Bose–Einstein distribution will reduce (depending on the temperature) the real integration range. This is in contrast to the Landauer conductance for electrons, where, strictly speaking, only states near the Fermi energy are playing a role.
