3.2.1. Proof-of-Principle: One-Dimensional Atomic Chain

To illustrate this approach, a one-dimensional (1D) chain with *N* atoms interacting via force constants with strength *λ* = 1.0 eV/μÅ was considered (see Figure 9a) for *N* = 4 [80,147]. The atomic masses of the atoms was set to 1.0 μ, and only nearest-neighbor interactions were considered.

**Figure 9.** (**a**) Scheme of the one-dimensional atomic chain studied in this work for the case of *N* = 4 atoms, with the filled atoms representing the beginning of the heat baths. *λ* and *η* are the spring force constants between the atoms and the coupling of the central region to the baths. Time-dependent NEGF approach: (**b**) Variation of the total energy of a dimer at *T*0 = 300 K after increasing the number of atoms in the one-dimensional atomic chain for different cut-off frequency (**left**) and *η* parameter (**right**). For comparison, we also plotted the energy values corresponding to the ideal harmonic oscillator case (dashed lines). Panel (**b**) is reproduced with permission from Ref. [147]. Copyright 2018 American Chemical Society.

First, a benchmark of the influence of spectral density parameters on the steady state properties was carried out. In Figure 9b, the dependence of the system energy at 300 K on *N* for different cut-off frequency *ωc* and *η* parameter is displayed [80,147]. The energy is compared to that obtained for an ideal harmonic oscillator in equilibrium, *E*C = ∑ *N i*=1 *ωi* (*n*(*<sup>ω</sup>i*) + 1/2), with *n*(*ω*) being the Bose–Einstein distribution function and *ωi* are the frequencies of the isolated central system. For fixed *η* = *λ*, increasing *ωc* leads to an increase of the system energy, since the spectral density includes more vibrational states. Contrarily, the energy gets closer to the harmonic oscillator values after reducing *η* for fixed *ωc* = 400 THz. This phenomena is expected because, for *η* → 0, the chain can be considered as an isolated harmonic system with no heat injection from the bath and *E*1Dchain = *EC*.

A similar effect is also found by analyzing the phonon transmission *<sup>τ</sup>p<sup>h</sup>*(*ω*) calculated along the lines of Section 2.1. Using the spectral density <sup>Λ</sup>(*ω*)=( *ω*<sup>2</sup> c*<sup>ω</sup>*)/(*ω*<sup>2</sup> + *ω*<sup>2</sup> c )Λ(0), one gets [147]:

$$\Sigma\_{\rm af}^{\varepsilon\_a}(t, t') = -\mathbf{i} \int\_{-\infty}^{\infty} \frac{d\omega}{2\pi} \left(\frac{\omega\_c^2 \omega}{\omega^2 + \omega\_c^2}\right) \Lambda\_a^{(0)} \left[\coth \frac{\omega}{2k\_B T\_a} \pm 1\right] e^{i\omega \left(t - t'\right)}.\tag{46}$$

Hence, the retarded self-energy in time domain is written as:

$$\Sigma^r\_{\mathfrak{a}}(t, t') = \Theta(t - t') 2\mathrm{i}\Lambda^{(0)}\_{\mathfrak{a}} \int\_{-\infty}^{\infty} \frac{d\omega}{2\pi} \left(\frac{\omega\_c^2 \omega}{\omega^2 + \omega\_c^2}\right) e^{i\omega \cdot (t - t')}.\tag{47}$$

By performing a Fourier transform of the previous equation and considering the counter-term, the self-energies are given by:

$$\Sigma^r\_{\mathfrak{a}}(\omega) = \left[ \mathbf{i} \frac{\omega\_c^2 \omega}{\omega^2 + \omega\_c^2} \right] \Lambda\_{\mathfrak{a}}^{(0)} = \left[ \Sigma^a\_{\mathfrak{a}}(\omega) \right]^\dagger. \tag{48}$$

The retarded Green-function finally reads:

$$\mathbf{G}^{r}(\omega) = \left[\omega^{2}\mathbf{I} - \mathbf{K} - \Sigma\_{L}^{r}(\omega) - \Sigma\_{R}^{r}(\omega)\right]^{-1},\tag{49}$$

where **K** is the force constant matrix of the one-dimensional chain. The transmission function *<sup>τ</sup>p<sup>h</sup>*(*ω*) is computed with Equation (23) and the steady heat flux is calculated by using the Landauer approach (see Equation (22)).

Figure 10a shows *<sup>τ</sup>p<sup>h</sup>*(*ω*) for different *N* with a cut-off at 400 THz. New transmission peaks appear for larger *N* due to emergence of new vibrational modes [80]. In addition, it was found that the maximum frequency *ωmax* with non-zero transmission depends on *N*. Thus, for *N* > 8, *ωmax* remains constant ( ∼195 THz). For *N* → <sup>∞</sup>, *<sup>τ</sup>p<sup>h</sup>* is constant and is zero for *ω* > *ωmax*, i.e., all modes have the same transmission probability *<sup>τ</sup>p<sup>h</sup>* = 1.0 (see also [33]). Figure 10b shows the influence of *η* for a dimer with *ωc* = 400 THz. *η* has a considerable influence on the phonon transmission. For *η* ≥ *λ*, *<sup>τ</sup>p<sup>h</sup>* is similar to a Gaussian and the dimer cannot be understood as a weakly coupled system anymore. For *η* ≤ 0.8 *λ*, on the contrary, two main transmission peaks corresponding to the two dimer modes can be resolved. For even weaker coupling, *<sup>τ</sup>p<sup>h</sup>* will yield two delta functions at these frequencies, correctly describing the vibrations of the system. This result confirms the analysis carried out for the system energy and shown in Figure 9b. The influence of the cut-off frequency on *<sup>τ</sup>p<sup>h</sup>* is weak compared to *η*: increasing *ωc* ten times only leads to a slight reduction of the phonon transmission at high frequencies and the frequency spectrum becomes wider (see Figure 10c) [80,147].

Based on the Landauer formalism, the temperature of the heat baths only appears in the Bose–Einstein distribution (see Equation (22)). However, in the TD-NEGF approach, the temperatures appear in the auxiliary-mode expansion of the self-energy. Once the system is in thermal equilibrium, a symmetric temperature bias Δ *T* = *TL* − *TR* = 2*ξT*0 (*ξ* > 0) is applied, with *T*0 being the mean temperature at which the system was previously equilibrated [80,147]. The left and right baths temperatures are expressed as *TL* = (1 + *ξ*)*<sup>T</sup>*0 and *TR* = (1 − *ξ*)*<sup>T</sup>*0. Figure 11a shows the steady heat flux as a function of *N* for different *η* (*<sup>ω</sup>c* = 400 THz). The heat flux values were obtained for a mean temperature of *T*0 = 300 K and *ξ* = 0.1. The values in both methods become closer after reducing *η* in agreemen<sup>t</sup> with the above discussed results in Figures 9b and 10b. Additionally, the heat flux converges for a given *N* for all *η*s [80,147].

**Figure 10.** Landauer approach: (**a**) Variation of the phonon transmission function *<sup>τ</sup>p<sup>h</sup>* of an one-dimensional atomic chain as a function of the number of atoms. (**b**) Influence of the coupling parameter on the phonon transmission function *<sup>τ</sup>p<sup>h</sup>* of an atomic dimer. (**c**) Variation of *<sup>τ</sup>p<sup>h</sup>* with respect to the cut-off frequency for *N* = 2 (**top**) and *N* = 4 (**bottom**). The figure is reproduced with permission from Ref. [147]. Copyright 2018 American Chemical Society.

**Figure 11.** Landauer approach: (**a**) Steady heat flux as a function of the number of atoms in the one-dimensional atomic chain for different *η* values. (**b**) Cut-off frequency *ωc* dependence of the steady heat flux for the atomic dimer at various temperatures bias Δ*T*. For comparison, we also plotted the values obtained using Landauer approach. The figure is reproduced with permission from Ref. [147]. Copyright 2018 American Chemical Society.

The influence of *ωc* on the steady heat flux was studied for *η* = *λ* (see Figure 11b). One sees that each approach displays a different behavior, i.e., the heat flux increases and decreases after increasing *ωc* for the TD-NEFG method and the Landauer approach, respectively [80,147]. This effect can be tuned by the value of Δ *T*. However, independently of Δ *T*, the heat flux for both approaches become closer with increasing *ωc*.
