3.2.2. Atomistic System: Carbon-Based Molecular Junctions

Medrano Sandonas et al. [147] showed that the TD-NEGF method reviewed in this section can be combined with atomistic methodologies for addressing the phonon dynamics in real systems. Due to the larger number of degrees of freedom, the matrix dimensions considerably increase and, hence, the computational cost. As typical examples, the work in Ref. [147] focused on poly-acetylene (PA, 4 atoms) and poly-ethylene (PE, 6 atoms) dimers connected to thermal baths (see Figure 12a for the case of PA dimer). The main difference between the two systems was the presence of double C bonds in PA compared with single bonds in PE. Both structural optimization and force constant calculations were performed with the Gaussian09 code [153]. Λ(0) *α* will thus take values corresponding to realistic bonds

between the central region and the reservoirs (for more details, see [80]). For both junctions, a cut-off *ω*c = 100 THz was used (roughly two times the maximum frequency of the vibrational spectrum). The number of poles in the auxiliary-mode expansion at 100 K, 300 K, and 500 K was 10, 8, and 4, correspondingly [80,147].

**Figure 12.** Time-dependent NEGF approach: (**a**) Energy density plot for molecular junctions made of poly-acetylene (PA) dimer. (**b**) Variation of the steady-state heat flux as a function of *ξ* for the PA dimer at different *T*0. Inset: Dynamics of the heat flux for both leads for PA and PE dimers after applying a temperature bias of Δ*T* = 60 K at *T*0 = 300 K. The figure is reproduced with permission from Ref. [147]. Copyright 2018 American Chemical Society.

To gain a deeper understanding of the thermal properties in the transient state, the energy density *<sup>D</sup>*(*<sup>E</sup>*, *t*) was defined as:

$$D(E, t) = \sum\_{i=1}^{N} (1/\gamma \sqrt{2\pi}) \exp[-\left((E - E\_i)/\gamma \sqrt{2}\right)^2] \tag{50}$$

with *γ* = 0.001 eV and {*Ei*} the set of eigenvalues of <sup>Z</sup>(*t*). In Figure 12a, the results for PA dimer during thermal equilibration at *T*0 = 300 K are displayed. As shown in the figure, all modes of the <sup>Z</sup>(*t*) matrix display very low energy at the beginning of the transient. The lowest lying modes gain then energy and reach a maximum at equilibrium. However, the eigenvalues in PE need a longer time to converge as compared to PA [147]. This difference arises from the different coupling strengths to the leads (related to the matrices Λ(0) *α* ). Consequently, the magnitude of the oscillations in the heat flux during the transient after applying a temperature bias is different, being larger for PA, as shown in the inset of Figure 12b. This difference in covalently bonded configurations also leads to a larger heat flux for the PA dimer [147]. Moreover, in agreemen<sup>t</sup> with linear response, the heat flux behaves nearly linear in *ξ* for different mean temperatures *T*0 (see Figure 12b) [80,147].
