2.3.3. Selective Molecular-Scale Phonon Filtering

We have recently proposed nano-junctions consisting of two colinear (6,6)-nanotubes (NT) joined by a central molecular structure as a potential molecular-scale phonon filter (see Figure 6a) [121]. The nanotubes act as heat baths kept at the same temperature. The left NT is considered as a reference bath with a broad phonon frequency spectrum, playing the role of a "source" of phonon modes. The molecular system in the central part is as a mode selector, and selected modes are then propagated to the right contact [80]. Besides carbon NTs, boron-nitride (BN) and silicon carbide (SiC) nanotubes as right baths were also considered. The filtering capability of the device thus depends on a mode-specific propagation resulting from the combined effect of molecular vibrations selection rules and the overlap of the contact spectral densities with the molecular region.

The phonon transport problem was treated using the previously described Green's function technique as implemented in the PHONON tool. Figure 6b,c shows the influence of interconnecting chains on the filtering effect for the case where both thermal baths are (6,6)-CNTs with *ωD* ∼ 1800 cm<sup>−</sup>1. As bridging molecular systems four parallel chains of ethylene, benzene, and azobenzene were chosen [121]. In Figure 6b, spectral gaps emerge in the transmission induced by the presence of the molecular chains. The overall transmission of the junctions is reduced by roughly a factor of four when compared with the infinite CNT due additional phonon scattering effects at the interfaces.

The influence of azobenzene chains is stronger comparing to the other monomers. Thus, chains with only two monomers already induce phonon gaps and filter out roughly half of the spectral range (see the lower panel of Figure 6b). Contrary to the benzene case, the transmission at low frequencies was strongly reduced, a result probably related to the lower number of modes and additional scattering at lower frequencies induced by larger structural distortions [141,142]. As a result, azobenzene-based junctions display the lowest thermal conductance, *<sup>κ</sup>p<sup>h</sup>* (see Figure 7b for only CNT-based leads). In brief, it becomes clear that channel selection and phonon filtering can be strongly controlled by the chemical composition of the bridge [121].

**Figure 6.** (**a**) Schematic representation of the nanoscale phonon filter proposed in Ref. [121]. A two-terminal junction is considered, where the role of the thermal baths is played by two semi-infinite (6,6) nanotubes (CNT, BNT, SiCNT) which are bridged by molecular chains consisting of ethylene, benzene, and azobenzene monomers. *ωD* represents the Debye frequency in each nanotube. (**b**) Phonon transmission functions *<sup>τ</sup>p<sup>h</sup>* for benzene-based junctions. We also added the plot corresponding to the phonon transmission function of and infinite CNT (grey) and a single infinite molecular chain of benzene monomers (brown). Highlighted with dashed-line circles are the regions where phonon gaps clearly develop by increasing the chain length. (**c**) Variation of *<sup>τ</sup>KL*(*j* = 4, *N*) as a function of the number of monomers (*N*) in the different studied molecular junctions by considering both thermal baths made of (6,6)-CNTs. Each junction consists of four molecular chains in parallel, *j* = 4. The figure is reproduced with permission from Ref. [121]. Copyright 2019 American Chemical Society.

To quantify the deviation from the "source of modes" distribution (or degree of filtering), quantified in our case by the transmission spectrum of the CNT, a key design magnitude *<sup>τ</sup>*KL(*j*, *N*) was defined as [121]:

$$\tau\_{\rm KL}(j\_\prime N) = \frac{1}{\omega\_{\rm D}^{\rm CNT}} \int\_0^{\omega\_{\rm D}^{\rm CNT}} d\omega \,\tau\_{\rm CNT}(\omega) \ln \frac{\tau\_{\rm CNT}(\omega)}{\tau\_{\rm MJ,j}(\omega)},\tag{26}$$

with the index *j* referring to the number of molecular chains in parallel interconnecting the two thermal baths and *N* the number of monomers in one chain. *<sup>τ</sup>*CNT(*ω*) and *<sup>τ</sup>*MJ,j(*ω*) are the corresponding transmission functions for an infinite CNT and for the CNT-molecule junction containing *j* molecular chains in parallel. A perfect filter (zero transmission) yields *<sup>τ</sup>*KL(*j*, *N*) → <sup>∞</sup>, while no filtering at all yields *<sup>τ</sup>*KL(*j*, *N*) = 0. As shown in Figure 6c, the azobenzene junctions display the highest *τ*KL (i.e., highest filter efficiency) due to the efficient blocking of high-frequency modes above 1000 cm<sup>−</sup>1.

The ethylene-based chains are less efficient, but still display a larger effect than the benzene chains; this is mostly related to two issues: the complete filtering of frequencies larger than 1500 cm<sup>−</sup><sup>1</sup> and the presence of a relatively large phonon gap between 500 cm<sup>−</sup><sup>1</sup> and 750 cm<sup>−</sup>1. *τ*KL for benzene-based junctions also increases with the length of the molecular chain and shows a tendency to saturate for *N* > 12 monomers.

**Figure 7.** Phonon thermal conductance *<sup>κ</sup>p<sup>h</sup>* as a function of the temperature. (**a**) *<sup>κ</sup>p<sup>h</sup>* values of the infinite CNT and BNNT as well as of the CNT-BNNT junctions. The remaining panels show the thermal conductance in the different junction types with a chain length corresponding to four monomers for: (**b**) ethylene; (**c**) benzene; and (**d**) azobenzene. The figure is reproduced with permission from Ref. [121]. Copyright 2019 American Chemical Society.

Figure 7 highlights the influence of changing the material of the right bath on the phonon filtering. In these calculations, all chains consist of four monomers. The thermal conductance of CNT-BNNT junctions is reduced when compared to the perfect tubes due to interface scattering (see Figure 7a [143–146]). By inserting the molecular system, a reduction of the thermal conductance by roughly a factor of 3–4 is produced, the effect being more pronounced for the azobenzene junction [121]. For azobenzene-based molecular junctions, the thermal conductance barely changes when going from CNT-CNT to CNT-BNNT (see Figure 7d). This is a consequence of the strong suppression of high frequency transport channels in the transmission. The saturation of the thermal conductance is determined by the nanotube with the smaller Debye cutoff (going from CNT (*<sup>ω</sup>D* ∼ 1800 cm<sup>−</sup>1) to BNNT (*<sup>ω</sup>D* ∼ 1400 cm<sup>−</sup>1) to SiCNT (*<sup>ω</sup>D* ∼ 880 cm<sup>−</sup>1)), the value of the saturation point is, however, influenced by the specific composition of the molecular chains.

### **3. Atomistic Framework for Time-Dependent Thermal Transport**

A novel atomistic approach able to treat transient phonon transport was recently developed by Medrano Sandonas et al. [147]. The approach is based on the solution of the equation of motion for the phonon density matrix *<sup>σ</sup>*(*t*), calculated within the NEGF formalism, by using an auxiliary-mode approach [80,147]. The latter has been previously used for time-resolved electron transport [55–57]. Unlike recent related approaches with limited application range (in terms of an atomistic treatment of the underlying system) [61–63], our method can be efficiently combined with an atomistic description as implemented in standard first-principle or parameterized approaches.

The basic structure of this approach is shown in Figure 8. Two thermal baths made of non-interacting harmonic oscillators in thermal equilibrium are contacted to a scattering region, whose vibrational features are assumed to be well represented by a quadratic Hamilton operator. The total system is described by the Hamiltonian:

$$H = H\_{\mathbb{C}} + \sum\_{ak} \left( \frac{1}{2} p\_{ak}^2 + \frac{1}{2} \omega\_{ak}^2 u\_{ak}^2 \right) + \sum\_{a,k} \frac{1}{2} \left( \mathbf{u}^T \cdot \mathbf{V}\_{ak} \boldsymbol{\mu}\_{ak} + \boldsymbol{\mu}\_{ak} \mathbf{V}\_{ak}^T \cdot \mathbf{u} \right) \tag{27}$$

**Figure 8.** Schematic representation of the target molecular junctions by using the TD-NEGF approach. A molecular system is connected to two harmonic thermal baths, which are the source for the heat flow in the molecule. The figure is reproduced with permission from Ref. [147]. Copyright 2018 American Chemical Society.

The first term *H*C = (1/2)**p***<sup>T</sup>* · **p** + (1/2)**u***<sup>T</sup>* · **K**eff · **u** is the Hamiltonian of the central domain, **u** is a column vector consisting of all the displacement variables in the region, and **p** contains the corresponding momenta. Both vectors have length *N*, with *N* being the number of degrees of freedom in the central region. We chose renormalized displacements *ui* = √*mixi*, where *mi* is the mass associated to the *i*th vibrational degree of freedom, and *xi* is the actual displacement having the dimension of length [80,147]. The effective force-constant matrix **K**eff = **K** + **K**ct has dimension *N* × *N*, and includes the force constant matrix of the central region, **K**, and a counter-term **K**ct [80,147]. The index *α* ∈ {*<sup>L</sup>*, *R*} labels the left (L) and right (R) heat baths and *k* denotes their vibrational modes with frequency *ωαk*. The second term of Equation (27) is the Hamiltonian of the heat bath. The last term represents the interaction between the central region and the baths, given by coupling vectors **V***α<sup>k</sup>* which are assumed to vanish before time *t*0 → <sup>−</sup>∞. Written in this form, the coupling leads to a renormalization of the bare force-constant matrix, which can be canceled by the above introduced counterterm **K**ct = ∑*<sup>α</sup>*,*<sup>k</sup>*(**<sup>V</sup>***α<sup>k</sup>* · **<sup>V</sup>***Tα<sup>k</sup>*)/*ω*2*α<sup>k</sup>* [80,147]. Hence, the coupling to the thermal baths will introduce dissipation but no shift in the vibrational spectrum [148]. The equations of motion (EOM) for the central region (**u** and **p**) and normal modes of one lead (*uk*) read [80,147]:

$$\frac{\partial}{\partial t} \left( \begin{array}{c} \mathbf{u} \\ \mathbf{p} \end{array} \right) = \mathcal{K}\_{\mathrm{eff}} \cdot \left( \begin{array}{c} \mathbf{u} \\ \mathbf{p} \end{array} \right) - \sum\_{k} \mu\_{k} \left( \begin{array}{c} \mathbf{0} \\ \mathbf{V}\_{k} \end{array} \right). \tag{28}$$

The 2*N* × 2*N* dimensional auxiliary matrices (denoted by caligraphic symbols) are given by the expressions [80,147]:

$$\mathcal{I} \equiv \begin{pmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{pmatrix}, \quad \mathcal{Q} \equiv \begin{pmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{I} & \mathbf{0} \end{pmatrix}, \quad \mathcal{K}\_{\mathrm{eff}} \equiv \begin{pmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{K}\_{\mathrm{eff}} & \mathbf{0} \end{pmatrix} \equiv \begin{pmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{K} & \mathbf{0} \end{pmatrix} + \begin{pmatrix} \mathbf{0} & \mathbf{0} \\ -\mathbf{K}\_{\mathrm{tr}} & \mathbf{0} \end{pmatrix}.$$

The energy of the central region can be written in terms of the phonon density matrix *σ*(*t*) = <sup>i</sup>G<sup>&</sup>lt;(*<sup>t</sup>*, *t*), with G<sup>&</sup>lt;(*<sup>t</sup>*, *t*) being a lesser Green's function [80,147]:

$$E\_{\mathbb{C}}(t) = \frac{1}{2} \text{Tr} \left\{ \mathcal{K}\_{\text{eff}}^T \cdot \mathcal{Q} \cdot \sigma(t) \right\},\tag{29}$$

and the total energy is *ETot*(*t*) = *EC*(*t*) + *Ebath*(*t*). In the absence of external forces, the total energy is conserved, and the heat current coming from the heat baths can be defined as *J*(*t*) = − *∂∂t <sup>E</sup>*C(*t*). Hereafter, - = 1 is taken. The time evolution of the heat flux is related to the lesser GF as [62,63]:

$$\mathcal{G}^{<}(t, t') = -\mathrm{i} \left( \begin{array}{c} \langle \mathbf{u}(t') \mathbf{u}^{T}(t) \rangle \\ \langle \mathbf{p}(t') \mathbf{u}^{T}(t) \rangle \end{array} \begin{array}{c} \langle \mathbf{u}(t') \mathbf{p}^{T}(t) \rangle \\ \langle \mathbf{p}(t') \mathbf{p}^{T}(t) \rangle \end{array} \right) . \tag{30}$$

To obtain the time dependence of the lesser GF, Dyson's equation was derived (see Ref. [80] for details). Using Langreth's rules, the lesser GF can be written in terms of retarded and advanced GFs:

$$\mathcal{G}^{<}(t, t') = \int d\tau\_2 \int d\tau\_3 \, \mathcal{G}^{R}(t, \tau\_2) \cdot \mathcal{S}^{<}(\tau\_2, \tau\_3) \cdot \mathcal{G}^{A}(\tau\_3, t'). \tag{31}$$

For the latter, EOMs are given by [80,147]:

$$\frac{\partial}{\partial t} \mathcal{G}^{\mathbb{R}}(t, t') = \delta(t, t') \mathcal{Q} + \mathcal{K}\_{\mathrm{eff}} \cdot \mathcal{G}^{\mathbb{R}}(t, t') + \mathcal{Q} \cdot \int dt\_2 \, \mathcal{S}^{\mathbb{R}}(t, t\_2) \cdot \mathcal{G}^{\mathbb{R}}(t\_2, t'). \tag{32}$$

$$\frac{\partial}{\partial t'} \mathcal{G}^A(t, t') = \delta(t, t') \mathcal{Q}^T + \mathcal{G}^A(t, t') \cdot \mathcal{K}\_{\mathrm{eff}}^T + \int dt\_2 \, \mathcal{G}^A(t, t\_2) \cdot \mathcal{S}^A(t\_2, t') \cdot \mathcal{Q}^T. \tag{33}$$

S*<sup>R</sup>*,*A*,<sup>&</sup>lt; denotes the respective self-energy. Setting *t* → *t* in Equation (31), the EOM for lesser GF becomes:

$$\begin{split} \frac{\partial}{\partial t} \mathcal{G}^{<}(t,t) &= \mathcal{K}\_{\mathrm{eff}} \cdot \mathcal{G}^{<}(t,t) + \mathcal{G}^{<}(t,t) \cdot \mathcal{K}\_{\mathrm{eff}}^{T} \\ &\quad + \mathcal{Q} \cdot \left[ \int\_{\tau\_{2}}^{t} d\tau\_{2} \left( \mathcal{S}^{>}(t,\tau\_{2}) \mathcal{G}^{<}(\tau\_{2},t) - \mathcal{S}^{<}(t,\tau\_{2}) \mathcal{G}^{>}(\tau\_{2},t) \right) \right] \\ &\quad + \left[ \int\_{\tau\_{2}}^{t} d\tau\_{2} \left( \mathcal{G}^{>}(t,\tau\_{2}) \mathcal{S}^{<}(\tau\_{2},t) - \mathcal{G}^{<}(t,\tau\_{2}) \mathcal{S}^{>}(\tau\_{2},t) \right) \right] \cdot \mathcal{Q}^{T}, \end{split} \tag{34}$$

where the thermal current matrices <sup>Π</sup>*α*(*t*) are defined as:

$$\Pi\_{\mathfrak{a}}(t) = \int\_{t2}^{t} d\tau\_{2} \, \left( \mathcal{G}^{>}(t, \tau\_{2}) \mathcal{S}\_{\mathfrak{a}}^{<}(\tau\_{2}, t) - \mathcal{G}^{<}(t, \tau\_{2}) \mathcal{S}\_{\mathfrak{a}}^{>}(\tau\_{2}, t) \right) . \tag{35}$$

For harmonic baths, S<sup>&</sup>lt;,<sup>&</sup>gt; *α* (*t*, *t*) can be obtained as:

$$\mathcal{S}\_{\mathfrak{a}}^{\mathcal{G}}(t, t') = -i \int\_0^\infty \frac{d\omega}{\pi} \left[ \coth \frac{\omega}{2k\_B T\_\mathfrak{a}} \cos \omega (t - t') \pm i \sin \omega (t - t') \right] \mathcal{L}\_{\mathfrak{a}}(\omega) \,. \tag{36}$$

*Tα* is the bath temperature and L*α*(*ω*) is the spectral density of reservoir *α* [148,149]. The EOM for the phonon density matrix reads [80,147]:

$$\frac{\partial}{\partial t}\boldsymbol{\sigma}(t) = \boldsymbol{\mathcal{K}}\_{\text{eff}} \cdot \boldsymbol{\sigma}(t) + \boldsymbol{\sigma}(t) \cdot \boldsymbol{\mathcal{K}}\_{\text{eff}}^{T} + \mathop{\rm in}\_{a \in \{L, \mathbb{R}\}} \left(\boldsymbol{\Pi}\_{a}(t) \cdot \boldsymbol{\mathcal{Q}}^{T} - h.c.\right). \tag{37}$$
