*3.1. Auxiliary-Mode Approach*

An auxiliary-mode approach was used to expand the self-energies S<sup>&</sup>lt;,<sup>&</sup>gt; *α* (*t*, *t*) in exponential functions to achieve a numerically efficient implementation to calculate the time evolution of <sup>Π</sup>*α*(*t*). The approach has been previously implemented for electrons [55–57] and for vibrations [150,151]. Since the cos and sin functions in Equation (36) can be easily expressed in terms of exponential functions, the only term which has to be considered is the hyperbolic cotangent. Different schemes have been proposed to obtain a suitable pole decomposition of this term [150]. To bypass the slow convergence of the so-called Matsubara decomposition, a more advanced pole decomposition was suggested by Croy and Saalmann [55], based on a partial fraction decomposition method that displays a faster convergence. However, one needs in this method high-precision arithmetic to compute the poles correctly, so that it is of advantage to have a pole decomposition of the coth function with purely imaginary poles [80]. A Pade decomposition method was recently proposed by Hu et al. [152], which shows very rapid convergence: the hyperbolic cotangent can be written in terms of simple poles as [150]:

$$\coth(\mathbf{x}) \approx \frac{1}{\mathbf{x}} + \sum\_{p=1}^{N^p} \eta\_p \left( \frac{1}{\mathbf{x} - \tilde{\xi}\_p} + \frac{1}{\mathbf{x} - \tilde{\xi}\_p^\*} \right) \tag{38}$$

where *ηp* are residues and *ξ p* (with Im *ξ p* > 0) are the poles. Here, *NP* denotes the number of poles.

Following this auxiliary-mode approach, a spectral density with a Drude regularization (i.e., adding a cut-off frequency *ω*c) of the form <sup>L</sup>*α*(*ω*)=(*ω*2c*ω*)/(*ω*<sup>2</sup> + *ω*2c )L(0) *α* was considered [148]. L(0) *α* ≡ diag(Λ(0) *α* , **0**) contains the coupling matrix Λ(0) *α* between the central domain and the leads. This quantity contains the wide-band limit as a special case for *ω*c → ∞. Using a linear combination of Lorentzians, any spectral density can in principle be approximated. The Drude spectral density is inserted into the expression for the lesser/greater self-energies (see Equation (36)), and for *τ* = *t* − *t* > 0, it turns into [80,147]:

$$\mathcal{S}\_{\mathfrak{a}}^{\varepsilon\_{\varepsilon}^{\varepsilon}}(t, t') = \mathcal{L}\_{\mathfrak{a}}^{(0)}\left[\mathbb{C}\_{\mathfrak{a}}(\tau) \pm \text{sgn}(\tau) \frac{\omega\_{\varepsilon}^{2}}{2} e^{-\omega\_{\varepsilon}|\tau|}\right]. \tag{39}$$

with:

$$\mathbf{C}\_{a}(\tau) = -\mathrm{i} \int\_{-\infty}^{\infty} \frac{d\omega}{2\pi} \left( \frac{\omega\_{c}^{2} \omega}{\omega^{2} + \omega\_{c}^{2}} \right) \coth \frac{\omega}{2k\_{B}T\_{a}} \mathbf{e}^{\mathrm{i}\omega \tau},\tag{40}$$

As mentioned above, the goal of the auxiliary-mode approach is to expand the self-energies and *<sup>C</sup>α*(*τ*) in terms of exponentials. Using the Pade decomposition of the hyperbolic cotangent in Equation (40), one obtains [80,147]:

$$\mathcal{L}\_a(\tau) = -ik\_B T\_{\mathfrak{a}} \omega\_{\mathfrak{c}} e^{-\omega\_{\mathfrak{c}} \tau} - i \sum\_{p=1}^{N\_P} R\_{\mathfrak{a},p} \left( \omega\_{\mathfrak{c}} e^{-\omega\_{\mathfrak{c}} \tau} - \chi\_{\mathfrak{a},p} e^{-\chi\_{\mathfrak{a},p} \tau} \right),$$

with *<sup>R</sup><sup>α</sup>*,*<sup>p</sup>* = <sup>2</sup>*kBTαω*2*c <sup>ω</sup>*2*c* − *<sup>χ</sup>*2*α*,*<sup>p</sup> ηp* and *χα*,*<sup>p</sup>* = −i2*kBTαξ p*. Hence:

$$\mathcal{S}\_{a}^{<\_{s}>}(t,t') = \frac{\partial}{\partial t} \sum\_{p=0}^{N\_{\rm P}} a\_{a,p}^{<\_{s}>} e^{-b\_{\rm r,p}(t-t')} \mathcal{L}\_{a}^{(0)} = \frac{\partial}{\partial t} \mathcal{N}^{<\_{s}>}(t,t') \tag{41}$$

where the coefficients *a*<,<sup>&</sup>gt; *<sup>α</sup>*,*p* and *<sup>b</sup><sup>α</sup>*,*<sup>p</sup>* are related to the auxiliary-mode decomposition. For *τ* < 0, the coefficients are *a*<sup>∗</sup>,<,<sup>&</sup>gt; *<sup>α</sup>*,*p* and *<sup>b</sup>*<sup>∗</sup>*α*,*<sup>p</sup>* [147].

The phonon current matrices <sup>Π</sup>*α*(*t*) (see Equation (35)) can be written as [80,147]:

$$\Pi\_{\mathfrak{a}}(t) = \sum\_{p=0}^{N\_P} \left[ \Phi\_{\mathfrak{a}}^p(t) + a\_{\mathfrak{a},p}^{\*, < <} \mathcal{Q} \cdot \mathcal{L}\_{\mathfrak{a}}^{(0)} \right], \tag{42}$$

with:

$$\Phi\_{a}^{p}(t) = \int\_{t\_{0}}^{t} dt' \,\left(\mathcal{G}^{<}(t, t') \cdot \mathbb{K}\_{\mathrm{eff}}^{T} \cdot \mathcal{N}\_{a, p}^{> \circ}(t', t) - \mathcal{G}^{>}(t, t') \cdot \mathbb{K}\_{\mathrm{eff}}^{T} \cdot \mathcal{N}\_{a, p}^{<}(t', t)\right),$$

The EOM for <sup>Φ</sup>*pα*(*t*) is given by [80,147]:

$$\frac{\partial}{\partial t} \Phi\_{\mathfrak{a}}^{p}(t) = \mathcal{A}\_{\mathfrak{a}}^{p} \cdot \Phi\_{\mathfrak{a}}^{p}(t, t) - \mathcal{B}\_{\mathfrak{a}}^{p} \cdot \mathcal{K}\_{\mathrm{eff}}^{T} \cdot \mathcal{L}\_{\mathfrak{a}}^{(0)} + \mathcal{Q} \cdot \Omega\_{\mathfrak{a}}^{p}(t), \tag{43}$$

with A*pα* = K − *<sup>b</sup>*<sup>∗</sup>*α*,*p*I, B*pα* = <sup>i</sup>*ωcδp*,0*<sup>σ</sup>*(*t*) + *<sup>a</sup>*<sup>∗</sup>,<*α*,*<sup>p</sup>* Q. The functions <sup>Ω</sup>*pα*(*t*) = ∑*α* ∑*<sup>N</sup>*P *p*=<sup>0</sup> Ω*p p αα* (*t*) correlate the main features of both heat baths and their values are obtained by performing the time derivative of Ω*p p αα*(*t*), which is expressed as [147]:

$$\begin{split} \frac{\partial}{\partial t} \Omega\_{\boldsymbol{a}'\boldsymbol{a}}^{p'p}(t) &= \mathcal{C}\_{\boldsymbol{a}'\boldsymbol{a}}^{p'p} \mathcal{L}\_{\boldsymbol{a}'}^{(0)} \cdot \mathcal{Q} \cdot \mathcal{K}\_{\text{eff}}^{T} \cdot \mathcal{L}\_{\boldsymbol{a}}^{(0)} - D\_{\boldsymbol{a}'\boldsymbol{a}}^{p'p} \Omega\_{\boldsymbol{a}'\boldsymbol{a}}^{p'p}(t) \\ &- \omega\_{\boldsymbol{c}} \left( \delta\_{p',0} \mathcal{L}\_{\boldsymbol{a}'}^{(0)} \cdot \mathcal{K}\_{\text{eff}} \cdot \Phi\_{\mathbf{a}}^{p}(t) - \delta\_{p,0} \Phi\_{\mathbf{a}'}^{p',\dagger}(t) \cdot \mathcal{K}\_{\text{eff}}^{T} \cdot \mathcal{L}\_{\mathbf{a}}^{(0)} \right), \end{split} \tag{44}$$

where *Cp p αα* = *a*<sup>&</sup>lt;*<sup>α</sup>*,*p <sup>a</sup>*<sup>∗</sup>,>*α*,*<sup>p</sup>* − *a*<sup>&</sup>gt;*<sup>α</sup>*,*p <sup>a</sup>*<sup>∗</sup>,<*α*,*<sup>p</sup>* and *Dp p αα* = *<sup>b</sup>α*,*p* + *<sup>b</sup>*<sup>∗</sup>*α*,*<sup>p</sup>* are real numbers. Thus, the initial integro-differential equation for the reduced density matrix has been mapped onto a closed set of ordinary differential equations (Equations (37), (43), and (44)), which can be solved using, e.g., a fourth-order Runge–Kutta method [80,147].

The thermal current is calculated as [80,147]:

$$J(t) = -\frac{\mathrm{i}}{2} \mathrm{Tr} \left\{ \mathcal{K}\_{\mathrm{eff}}^T \cdot \mathcal{Q} \cdot \sum\_a \left[ \Pi\_a(t) \cdot \mathcal{Q}^T - h.c. \right] \right\}. \tag{45}$$
