Phononic heat transport in nanomechanical structures: steady-state and pumping

TL;DR

This paper develops a phase space Green's function approach to analyze phononic heat transport in nanomechanical systems, capturing atomic details and non-equilibrium effects, including steady-state and driven conditions.

Contribution

It introduces a novel formalism for phononic heat transport that accurately describes small molecular junctions and accounts for non-stationary and quantum thermodynamic properties.

Findings

Correctly recovers Landauer-like heat conductance in steady state

Derives expressions for non-stationary heat currents under periodic driving

Applies formalism to diatomic molecular junctions

Abstract

We study the heat transport due to phonons in nanomechanical structures using a phase space representation of non-equilibrium Green's functions. This representation accounts for the atomic degrees of freedom making it particularly suited for the description of small (molecular) junctions systems. We show that for the steady state limit our formalism correctly recovers the heuristic Landauer-like heat conductance for a quantum coherent molecular system coupled to thermal reservoirs. We find general expressions for the non-stationary heat current due to an external periodic drive. In both cases we discuss the quantum thermodynamic properties of the systems. We apply our formalism to the case of a diatomic molecular junction.

Equations

$$H(t)=\frac{1}{2}\,\vec{p}^{\;\text{T}}\cdot\vec{p}+\frac{1}{2}\,\vec{u}^{\;\text{T}}\cdot\hat{K}(t)\cdot\vec{u}\equiv\frac{1}{2}\,\mbox{\boldmath$\zeta$}^{\,\text{T}}\cdot\check{\mathcal{M}}(t)\cdot\mbox{\boldmath$\zeta$},$$

$$\mbox{\boldmath$\zeta$}=\begin{pmatrix}\vec{u}\\ \vec{p}\end{pmatrix}\qquad\text{and}\qquad\check{\mathcal{M}}(t)=\begin{pmatrix}\hat{K}(t)&\hat{0}\\ \hat{0}&\hat{I}\end{pmatrix},$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\mbox{\boldmath$\zeta$}=\check{\mathcal{Q}}\cdot\frac{\partial}{\partial\mbox{\boldmath$\zeta$}}H=-\check{\mathcal{K}}(t)\cdot\mbox{\boldmath$\zeta$},$$

$$\check{\mathcal{Q}}=\left(\begin{array}{cc}\hat{0}& \hat{I}\\ -\hat{I}& \hat{0}\end{array}\right)\text{and}\check{\mathcal{K}}(t)\equiv -\check{\mathcal{Q}}\cdot \check{\mathcal{M}}(t)=\left(\begin{array}{cc}\hat{0}& -\hat{I}\\ \hat{K}(t)& \hat{0}\end{array}\right).$$

$$\check{C}(\tau,\tau^{\prime})\equiv\frac{1}{\imath\,\hbar}\left\langle\mathbb{T}_{\mathcal{C}}\,\mbox{\boldmath$\zeta$}(\tau)\otimes\mbox{\boldmath$\zeta$}(\tau^{\prime})\right\rangle\equiv\begin{pmatrix}\hat{C}^{(uu)}&\hat{C}^{(up)}\\ \hat{C}^{(pu)}&\hat{C}^{(pp)}\end{pmatrix}(\tau,\tau^{\prime}),$$

$${\check{C}}^{>}(t,t^{\prime })$$

$${\check{C}}^{<}(t,t^{\prime })$$

$${\check{C}}^{\mathbb{T}}(t,t^{\prime })$$

$${\check{C}}^{\overline{\mathbb{T}}}(t,t^{\prime })$$

$${\check{C}}^r(t,t^{\prime })$$

$$=\theta (t-t^{\prime })\left({\check{C}}^{>}-{\check{C}}^{<}\right)(t,t^{\prime }),$$

$${\check{C}}^a(t,t^{\prime })$$

$$=\theta (t^{\prime }-t)\left({\check{C}}^{<}-{\check{C}}^{>}\right)(t,t^{\prime }),$$

$${\check{C}}^K(t,t^{\prime })$$

$$=\left({\check{C}}^{>}+{\check{C}}^{<}\right)(t,t^{\prime }).$$

$$\left(\check{\mathcal{I}}\frac{\partial }{\partial t}+\check{\mathcal{K}}(t)\right)\cdot {\check{C}}^{\gtrless }(t,t^{\prime })=0,$$

$$\left(\check{\mathcal{I}}\frac{\partial }{\partial t}+\check{\mathcal{K}}(t)\right)\cdot {\check{C}}^{\mathbb{T},\overline{\mathbb{T}}}(t,t^{\prime })=\pm \delta (t-t^{\prime })\check{\mathcal{Q}}.$$

$$\left(\check{\mathcal{I}}\frac{\partial }{\partial t}+\check{\mathcal{K}}(t)\right)\cdot {\check{C}}^K(t,t^{\prime })=0,$$

$$\left(\check{\mathcal{I}}\frac{\partial }{\partial t}+\check{\mathcal{K}}(t)\right)\cdot {\check{C}}^{r,a}(t,t^{\prime })=\delta (t-t^{\prime })\check{\mathcal{Q}},$$

$$\breve{\mathcal{C}}(t,t^{\prime })=$$

$$\equiv$$

$$\left(\check{\mathcal{I}}\frac{\partial }{\partial t}+\check{\mathcal{K}}(t)\right)\cdot \check{\mathcal{G}}(t,t^{\prime })=\delta (t-t^{\prime })\check{\mathcal{Q}},$$

$$\check{\mathcal{G}}(t,t^{\prime })\cdot \left(\check{\mathcal{I}}\overleftarrow{\frac{\partial }{\partial t^{\prime }}}+{\check{\mathcal{K}}}^{\text{T}}(t^{\prime })\right)=-\delta (t-t^{\prime })\check{\mathcal{Q}}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\check{\mathcal{G}}(t,t)$$

$$=-\check{\mathcal{K}}(t)\cdot \check{\mathcal{G}}(t,t)-\check{\mathcal{G}}(t,t)\cdot {\check{\mathcal{K}}}^{\text{T}}(t).$$

$$\breve{\mathcal{P}}\cdot \breve{\mathcal{C}}(t,t^{\prime })\cdot {\breve{\mathcal{P}}}^{\text{T}}=$$

$$\equiv$$

$$\underline{\check{\mathcal{G}}}[\omega ]={\int }_{-\infty }^{\infty }\mathrm{d}(t-t^{\prime }){\text{e}}^{ı\omega (t-t^{\prime })}\underline{\check{\mathcal{G}}}(t-t^{\prime }),$$

$$\underline{\check{\mathcal{G}}}[\omega ]=$$

$$=$$

$$\underline{\check{\mathcal{G}}}[\omega ]\equiv \left(\begin{array}{cc}{\check{\mathcal{G}}}^{(uu)}[\omega ]& {\check{\mathcal{G}}}^{(up)}[\omega ]\\ {\check{\mathcal{G}}}^{(pu)}[\omega ]& {\check{\mathcal{G}}}^{(pp)}[\omega ]\end{array}\right)=\left(\begin{array}{cc}\hat{G}[\omega ]& ı\omega \hat{G}[\omega ]\\ -ı\omega \hat{G}[\omega ]& \hat{G}[\omega ]\cdot \hat{K}\end{array}\right),$$

$$\hat{G}[\omega ]=({\omega }^2\hat{I}-\hat{K}{)}^{-1}.$$

$$\hat{G}[\omega ]=\frac{1}{2}{\int }_{-\infty }^{\infty }\frac{\mathrm{d}\bar{\omega }}{2\pi }\hat{J}(\bar{\omega })\left(\frac{1}{\omega -\bar{\omega }}-\frac{1}{\omega +\bar{\omega }}\right),$$

$$\hat{J}(\omega )=2\pi \sum _j\frac{1}{{\omega }_j}\delta (\omega -{\omega }_j)|j⟩⟨j|.$$

$${\hat{G}}^{r,a}[\omega ]$$

$$={\left[(\omega \pm ı0^{+}{)}^2\hat{I}-\hat{K}\right]}^{-1},$$

$${\hat{G}}^{\mathbb{T},\overline{\mathbb{T}}}[\omega ]$$

$$={\left[{\omega }^2\hat{I}-(\sqrt{\hat{K}}\mp ı0^{+}\hat{I}{)}^2\right]}^{-1}.$$

$$\begin{gathered} {\underline{\check{\mathcal{G}}}}^r(t-t^{\prime })=\theta (t-t^{\prime }) \\ \times \left(\begin{array}{cc}-\frac{\sin \left[\sqrt{\hat{K}}(t-t^{\prime })\right]}{\sqrt{\hat{K}}}& \cos \left[\sqrt{\hat{K}}(t-t^{\prime })\right]\\ -\cos \left[\sqrt{\hat{K}}(t-t^{\prime })\right]& -\sqrt{\hat{K}}\sin \left[\sqrt{\hat{K}}(t-t^{\prime })\right]\end{array}\right) \\ +\text{solution of homogeneous equation}, \end{gathered}$$

Full text

I Introduction

Significant progress has been recently achieved on the understanding of phononic heat transfer at the molecular level [1, 2, 3]. In addition to the investigation of fundamental aspects of the problem [1, 4], several authors have realized that phonons, usually regarded as an energy waste, can be manipulated and controlled to carry and process information. Exploring analogies with electrons and photons, theoretical proposals have been put forward aiming the fabrication of devices such as thermal diodes [5], thermal transistors [6, 7], and thermal logic gates [8], some of them already experimentally verified [9, 10, 11]. These ideas have given rise to the emerging field of phononics [2, 12].

The presence of an external time-dependent drive, such as an external force or time-varying thermal bath temperature, gives another interesting twist to the problem, making possible to explore non-equilibrium phenomena such as directed heat pumping and cooling [13, 14, 15, 2, 16, 17, 18, 19].

Early reports on the measurement of quantized thermal conductance in suspended nanostructures [20, 21] attracted attention to the field. More recently, ballistc thermal conductance has been experimentally studied in carbon nanotubes [22, 23, 24], silicon nanowires [25, 26], as well as molecular and atomic contacts [27, 28]. The experimental advances in these studies are remarkable and pose important challenges to the quantum theory of thermal conductance [24, 28].

One of the fundamental tools for the theoretical study of non-equilibrium properties of quantum systems is the non-equilibrium Green’s functions (NEGF) theory [29, 30]. This approach, originally developed for fermionic systems [31, 32], has been nicely adapted to describe the heat transfer in small junctions systems [33, 34, 35, 4, 36, 37, 38]. Despite its success, the implementation of the NEGF to calculate phonon heat currents driven by a temperature difference between source and drain still has some caveats, like the need to symmetrize the heat current to obtain the standard Landauer-like transmission formula [36, 37, 38]. The relevance of NEGF for phononics calls for a deeper and careful analysis of the formalism.

The purpose of the paper is twofold. First, we present a rigorous method for the description of quantum thermal transport properties due to phonon or atomic degrees of freedom using nonequilibrium Green’s function in phase space. We show that our formal developments solve the problems of the previous works [36, 37, 38] and recover the well known Landauer-like formula for the stationary heat current in the ballistic regime [39, 40, 41, 42, 43, 44]. Second, we extend the formalism to address systems under the influence of a time dependent drive. As an example, we derive general expressions for the heat current pumped by an external time-dependent periodic potential for a system coupled to two thermal reservoirs at the same temperature.

We show how to apply our method by analyzing the steady-state heat transport properties of a diatomic molecule coupled to thermal reservoirs by semi-infinite linear harmonic chains. Next, we study the heat current pumped through the system due to a time-dependent driving force and discuss its thermodynamic properties.

The paper is organized as follows: In Sec. II we introduce the phase space representation of the Green’s functions on which our derivations are built. We begin Sec. III by presenting the model Hamiltonian addressed in this study. We then use the Green’s function formalism to derive expressions for the thermal current due to a source-drain temperature difference and the heat current pumped by an external periodical drive of the system atomic degrees of freedom. In Sec. IV, we apply our results to the simple model of a diatomic molecular junction. Finally, we present our conclusions in Sec. V.

II Green’s functions in phase space

In this the section we use a phase space representation of non-equilibrium Green’s functions [45, 46]. We show that this representation is very convenient for a canonical quantization of the displacements $\vec{u}\equiv(u_{1},\ldots,u_{n})$ and their canonical conjugated momenta $\vec{p}\equiv(p_{1},\ldots,p_{n})$ in a $2n$-dimensional phase space.

Let us consider a quadratic Hamiltonian expressed in terms of space phase variables $(\vec{u},\vec{p})$ representing a system of coupled oscillators. The model Hamiltonian reads

[TABLE]

where, for the sake of compactness, we assume that the masses are identical and have unit value. $\hat{K}(t)$ is the force constant matrix that represents the couplings of the oscillators network. The dynamic variable $\zeta$ and the matrix $\check{\mathcal{M}}$ have the symplectic structure

[TABLE]

where $\hat{I}$ is the identity matrix.

The equation of motion for ${\bm{\zeta}}$ reads

[TABLE]

We define the phase space correlation functions $\hat{C}(\tau,\tau^{\prime})$ on the Keldysh contour [29] as

[TABLE]

where $\imath\hbar\,\hat{C}^{(\alpha\beta)}\equiv\langle\mathbb{T}_{\mathcal{C}}\,\vec{\alpha}(\tau)\otimes\vec{\beta}(\tau^{\prime})\rangle$. The correlation functions $\hat{C}^{(\alpha\beta)}(\tau,\tau^{\prime})$ are a straightforward phase space generalization of standard Green’s functions [29, 30], as we discuss below.

As standard [29], the greater, lesser, time-ordered, and anti-time-ordered correlations functions read

[TABLE]

where $(\check{C}^{\,\mathbb{T}}+\check{C}^{\,\overline{\mathbb{T}}}-\check{C}^{>}-\check{C}^{<}\,)(t,t^{\prime})=0$.

Alternatively, the correlation functions can be represented by their retarded $\check{C}^{r}$, advanced $\check{C}^{a}$, and Keldysh $\check{C}^{K}$ components, namely

[TABLE]

Using Eqs. (3a) and (5) we obtain the equations of motion for $\check{C}^{\gtrless}(t,t^{\prime})$ and $\check{C}^{\mathbb{T},\overline{\mathbb{T}}}(t,t^{\prime})$, namely

[TABLE]

Similarly, using Eqs. (3a) and (7), we show that $\check{C}^{\text{K}}(t,t^{\prime})$ and $\check{C}^{r,a}(t,t^{\prime})$ satisfy

[TABLE]

where $\check{\mathcal{I}}$ is the $2n\times 2n$ identity matrix. To obtain Eq. (8a), we use the identity $\check{C}^{>}(t,t)-\check{C}^{<}(t,t)=\check{\mathcal{Q}}$, that follows from the canonical commutations relations.

To make the notation compact, we write the correlation function in a block structure as (Keldysh space)$\,\otimes\,$(symplectic space) in its irreducible representation, namely

[TABLE]

where $\sigma_{1}$ is the first Pauli matrix. Note that $\check{\mathcal{G}}(t,t^{\prime})$ has also a symplectic structure and satisfies (by inspection) the equation of motion

[TABLE]

with a self-adjoint equation

[TABLE]

Using Eqs. (10) and (11) we obtain the following identity

[TABLE]

Performing the Keldysh rotation [30] in Eq. (II), we obtain a reducible representation of the correlation function in terms of the quantities defined in Eq. (6) as

[TABLE]

where $\breve{\mathcal{P}}=\frac{1}{\sqrt{2}}\left(I_{2}+\imath\,\sigma_{2}\right)\otimes\check{\mathcal{I}}$ and $\sigma_{2}$ is the second matrix of Pauli.

Let us now introduce the frequency representation of the correlation functions. Assuming time translational invariance, i.e., that the matrix $\check{\mathcal{K}}$ does not depend on time, one defines $\underline{\check{\mathcal{G}}}[\omega]$ in terms of the Fourier transform

[TABLE]

for $\underline{\check{\mathcal{G}}}(t-t^{\prime})=\check{\mathcal{G}}(t,t^{\prime})$. We study the time-dependent problem in Sec. III.2.

By inserting Eq. (14) in (10) [or in (11)], we write

[TABLE]

where

[TABLE]

with

[TABLE]

Equation (16) has been obtained in Ref. 36 by directly taking the Fourier transform of the displacement $\{u_{i}\}$ and the canonically conjugate momentum operators $\{p_{i}\}$. We note that despite being very appealing, this straightforward procedure is formally problematic, since the canonical commutation relations $\left[u_{i}(t),p_{j}(t)\right]$ can not be consistently defined in the frequency domain (see Appendix A for more details). This problem can be circumvented [38] by performing the Fourier transform of the phase space correlation functions, as described above.

The Green’s function $\hat{G}[\omega]$ can be represented as

[TABLE]

where the spectral operator $\hat{J}(\bar{\omega})$ is

[TABLE]

Here we have used that $\hat{K}$ is a positive-semidefinite matrix [47], which satisfies $\hat{K}|j\rangle=\omega_{j}^{2}\,|j\rangle$ with $\omega_{j}\geqslant 0$ (recall that $\langle j|j^{\prime}\rangle=\delta_{j,j^{\prime}}$ and $\sum_{j}|j\rangle\langle j|=\hat{I}$).

The general expression (18) does not distinguish the retarded, advanced, ordered, and anti-ordered components of $\hat{G}[\omega]$. A proper representation of the components requires a regularization around the poles $\omega=\pm\bar{\omega}$ of Eq. (18), namely

[TABLE]

The Green’s functions $\hat{G}^{r,a}(t,t^{\prime})$ and $\hat{G}^{\mathbb{T},\overline{\mathbb{T}}}(t,t^{\prime})$ are obtained by the inverse Fourier transform of Eqs. (20) and are consistent with Eqs. (5) and (6), as they should.

Substituting Eqs. (20a) and (16) in the inverse Fourier transform Eq. (14), we write the retarded component of $\underline{\mathcal{G}}(t-t^{\prime})$ as

[TABLE]

Similarly, the ordered and anti-ordered components read

[TABLE]

which satisfy $\check{\mathcal{G}}^{\mathbb{T}}(0^{\pm})-\check{\mathcal{G}}^{\overline{\mathbb{T}}}(0^{\pm})=\pm\check{\mathcal{Q}}$.

The Keldysh component of the correlation function is, in general, more demanding to obtain. As standard, the exception is the equilibrium case. In this limit, the fluctuation-dissipation theorem [48] relates the Keldysh component of the correlation function of a bosonic system to its retarded and advanced components as

[TABLE]

where $f(\omega)=\left(\text{e}^{\beta\hbar\omega}-1\right)^{-1}$ is the Bose-Einstein distribution function. One can also write

[TABLE]

where $\sigma=\pm 1$ and

[TABLE]

As a result, the equilibrium lesser and greater Green’s functions are given by

[TABLE]

III Model Hamiltonian

In this section, we describe the heat transport properties of a molecular junction modeled by a central region $C$ representing a nanostructure coupled by multiple leads connected to reservoirs in thermal equilibrium [36, 37]. We recall that we only consider thermal transport due vibrational degrees of freedom, which is the dominant mechanism in insulator systems.

This partition scheme allows one to write the general Hamiltonian of Eq. (1) as

[TABLE]

where

[TABLE]

correspond to the Hamiltonian of the $\alpha$-lead, central region and tunneling, respectively. We define the decoupled Hamiltonian $H_{a}^{0}$ corresponding to the $a$-partition as

[TABLE]

The force constant matrix in Eq. (1) is decomposed as $\hat{K}(t)=\hat{K}^{0}+\hat{V}(t)$, where $\hat{K}^{0}$ gives the dynamical matrix of the decoupled partitions

[TABLE]

and $\hat{V}(t)$ corresponds to the coupling between different partitions, namely

[TABLE]

These definitions allow us to write the tunneling Hamiltonian $H_{T}(t)$ as

[TABLE]

where $\vec{u}\equiv\left[\bigoplus_{\alpha}\vec{u}_{\alpha}\right]\oplus\vec{u}_{C}$. Note that $\hat{V}=\hat{V}^{\text{T}}$ and therefore $V_{\alpha C}=V_{C\alpha}^{\text{T}}$ for all $\alpha$ terminals.

The model Hamiltonian in Eq. (26) includes $V_{aa}$ ($a=\alpha,C$) terms that have not been explicitly accounted for by previous works [36, 37, 38]. Neglecting $V_{aa}$ can be problematic for the consistency of NEGF. This can be seen using the adiabatic switch-on picture, the standard implementation of NEGF in the steady-state regime (A discussion of different implementation schemes can be found, for instance, in Ref. [49]). The absorption of $V_{aa}$ into $K^{0}_{\alpha\alpha}$ modifies the free Green’s functions making their calculation troublesome. This issue becomes clear in the formal development below as well as in the applications discussed in Sec. IV.

To discuss the thermodynamic properties of the system it is convenient to describe the molecular junction as formed by reservoirs coupled to an extended central region, which we refer to as “molecule”. Accordingly, we write Eq. (26) as

[TABLE]

where the molecule Hamiltonian reads

[TABLE]

The energy of the extended molecule is defined as $E_{M}(t)\equiv\left\langle H_{M}(t)\right\rangle$, namely

[TABLE]

where the components of lesser functions are explicit given by

[TABLE]

with $a,b=\{C,\alpha\}$, in line with Eq. (4).

One can define the thermal current flowing through an open molecule connected to multiple reservoirs by comparing its energy variation

[TABLE]

with the energy continuity equation, expressed as

[TABLE]

where one associates $J_{\alpha}(t)$ to the thermal current from $\alpha$-reservoir into the molecule and $\Phi(t)$ is power developed by the ac sources (or drives) in the molecule. Hence, by inspection one infers that

[TABLE]

and

[TABLE]

Using the equation-of-motion method [48], we write the thermal current from $\alpha$-reservoir into the molecule in terms of the correlation functions as

[TABLE]

while the power developed by the external time-dependent drives reads

[TABLE]

In the following subsections we study separately the steady-state transport ($\dot{V}_{ab}=0$) and the heat transport due to pumping by an external drive ($\dot{V}_{ab}\neq 0$) for $a,b=\{\alpha,C\}$.

III.1 Steady-state transport

Let us now calculate the steady-state thermal current flowing from the $\alpha$-lead due to a temperature difference in the reservoirs. Here, we consider the heat current expression (38) for a time-independent coupling matrix $\hat{V}$.

Since the Hamiltonian does not explicitly depends on time, it is convenient to work in the frequency representation. The Fourier transform of $C^{<(up)}_{C\alpha}(t,t^{\prime})$ is

[TABLE]

where $C_{C\alpha}^{<(up)}[\omega]=\imath\,\omega\,G^{<}_{C\alpha}[\omega]$. Substituting Eq. (42) into Eq. (38), we cast the steady-state heat current as

[TABLE]

The system Green’s function \hat{G}[\omega]=\big{(}\omega^{2}\,\hat{I}-\hat{K}\big{)}^{-1} satisfies the Dyson equation

[TABLE]

where $\hat{K}=\hat{K}^{0}+\hat{V}$ and \hat{g}[\omega]=\big{(}\omega^{2}\,\hat{I}-\hat{K}^{0}\big{)}^{-1}. Note that the free Green’s function $\hat{g}[\omega]$ is block diagonal in the partitions.

From Eq. (III.1) we obtain

[TABLE]

where, for notational convenience, we introduce an effective embedding self-energy

[TABLE]

and an effective free Green’s function

[TABLE]

where $g_{a}[\omega]=(\omega^{2}I_{a}-K^{0}_{a})^{-1}$.

In Sec. IV.1 and in Appendix B we discuss the importance of including $V_{aa}$ in the surface Green’s function.

For $a=\alpha$, it corresponds to Green’s function in thermal equilibrium with the $\alpha$-reservoir at a temperature $T_{\alpha}$. Hence, using Eqs. (20a) and (25) we write obtain

[TABLE]

where $\imath A_{\alpha}(\omega)\equiv g^{r}_{\alpha}[\omega]-g^{a}_{\alpha}[\omega]$ is the $\alpha$-lead “free” spectral function and $f_{\alpha}(\omega)=\left(\text{e}^{\beta_{\alpha}\hbar\omega}-1\right)^{-1}$ with $\beta_{\alpha}=1/k_{B}T_{\alpha}$.

In general, the retarded and advanced surface Green’s functions ${g}_{\alpha}^{r,a}[\omega]$ are computed by decimation techniques [50, 36].

The lesser components of $G_{\alpha C}$ and $G_{C\alpha}$ are obtained by applying the Langreth rules [29, 48] to Eq. (45). By inserting the result in Eq. (43), we obtain

[TABLE]

The self-energies are given in terms of

[TABLE]

where $K_{\alpha}=K_{\alpha}^{0}\,+\,V_{\alpha\alpha}$ and $\imath\,\tilde{A}_{\alpha}(\omega)=\tilde{g}^{r}_{\alpha}[\omega]-\tilde{g}^{a}_{\alpha}[\omega]$. Hence,

[TABLE]

By expressing the self-energies in terms of the line width functions, we write the heat current as

[TABLE]

Applying the Langreth rules to Eq. (45c) and using Eq. (51), we obtain

[TABLE]

that are inserted in Eq. (52) to finally arrive at the steady-state heat current

[TABLE]

where

[TABLE]

rigorously obtaining the Landauer heat conductance that has been phenomenologically put forward [40] and adopted by several authors, see for instance, Refs. 43, 44, 51. As a consequence, the numerical implementation of the heat current $J_{\alpha}^{(S)}$ given by Eq. (54), is obviously the same as the one using the scattering matrix [43, 51].

The explicitly symmetric tunneling Hamiltonian $H_{T}(t)$, Eq. (27c), leads to an expression for the heat current $J_{\alpha}(t)$ with terms depending on both $V_{C\alpha}$ and $V_{\alpha C}$. This ensures that $J_{\alpha}(t)$ accounts for processes corresponding to the heat flow from the central region $C$ to the $\alpha$-lead as well as from $\alpha$ to $C$. Our result differs from the heat current derived by Wang and collaborators [36, 37, 38]. These authors derive the heat current using the Hamiltonian without explicitly taking into account processes associated to $V_{LC}$ (corresponding to $\alpha=L$). The obtained expression for heat current depends only on the hybrid Green’s function $G^{<}_{CL}$. Furthermore, the absence of $V_{\alpha C}$ (or $V_{C\alpha}$) in their Hamiltonian implies that the self-energy $\Sigma_{L}=V_{CL}\cdot g_{L}\cdot V_{LC}$ has to be introduced in a somewhat arbitrary manner. Moreover, Refs. [36, 37, 38] need the ad hoc symmetrization, $J=(J_{L}+J_{L}^{*}-J_{R}-J_{R}^{*})/4$, to obtain the well known Caroli formula for the transmission since the integrand of Eq. (43) is not purely real in the absence of $V_{\alpha C}$ (or $V_{C\alpha}$).

The transmission coefficient $\mathcal{T}_{\alpha\beta}(\omega)$ is interpreted as the probability of an energy $\hbar\omega$ to be transmitted from the reservoir $\alpha$ to the reservoir $\beta$ and has the same structure of the Meir-Wingreen formula [32] that describes the electronic conductance of fully coherent systems of non-interacting electrons.

It is straightforward to verify that $\mathcal{T}_{\alpha\beta}(\omega)=\mathcal{T}_{\beta\alpha}(\omega)$, which implies that in steady-state $J^{(S)}\equiv J_{L}^{(S)}=-J_{R}^{(S)}$. Hence, $\mathrm{d}E_{M}/\mathrm{d}t=0$ and, as expected, the molecule energy does not change in time.

III.2 Pumping transport

Let us now study the heat current in nanoscopic systems due to a time-dependent external drive, as motivated in the introduction. As in the stationary case, we employ the NEGF theory, since more standard approaches, like the Kubo-Greenwood one, are only suitable for bulk systems.

The analysis of heat currents in time-dependent systems is far more involved for bosonic degrees of freedom than for the electronic ones. In the latter case, the Fermi energy (and the corresponding Fermi velocity) establishes a characteristic time scale for the electronic dynamics. In experiments [52] the external driving is slow with respect to the electronic dynamics, which allows to approach the problem using the adiabatic approximation [53, 54, 55, 56, 57]. In the bosonic case there is no internal characteristic time scale and analytical progress has to resort on the assumption that the driving force is small to employ perturbation theory.

As an example of time-dependent transport, we study the case of periodically driven system in time. We assume that the coupling between regions depends on time as $\hat{V}(t)=\hat{V}+\varepsilon\,\hat{v}(t)$, where $\varepsilon$ is a dimensionless parameter.

The initial state is the fully connected molecule-leads system in equilibrium.

Defining an auxiliary matrix $\check{\mathcal{V}}(t)$ as

[TABLE]

we can write $\check{\mathcal{K}}(t)=\check{\mathcal{K}}-\check{\mathcal{Q}}\cdot\check{\mathcal{V}}(t)$ or, equivalently, $\check{\mathcal{M}}(t)=\check{\mathcal{M}}+\check{\mathcal{V}}(t)$. It follows from Eq. (10) that the Dyson’s equation reads

[TABLE]

where $\underline{\check{\mathcal{G}}}(t-t^{\prime})$ denotes the steady-state Green’s function transport, given by Eqs. (14) to (17). We consider $\varepsilon\ll 1$ and treat the problem using pertubation theory. This is an alternative approach to the Floquet analysis used in Refs. 17, 19.

We note that the Floquet method is extremely efficient, irrespective of coupling strength, provided the ratio between the band width and the driving frequency is not large, a condition that keeps the size of the Hilbert space computationally manageable. The opposite limit of small $\Omega$ is in general computationally prohibitive for this method. For electronic systems, however, it has been argued that if the characteristic single particle dwell time $\tau_{d}$ (evaluated at the Fermi energy) in the scattering region is much smaller than $1/\Omega$ only few harmonics of the perturbation are coupled. This allows for an effective truncation of the Hilbert space. The dwell time $\tau_{d}$ depends on the spectral density and on the strength of its coupling to the leads [58]. Since these quantities typically show a strong energy dependence, one has to verify if $\tau_{d}\Omega\ll 1$ is indeed fulfilled. In general the latter condition rules out the application of the Floquet approach for small $\Omega$ to a potentially large number of systems.

The Green’s function deviation from steady-state, $\delta\check{\mathcal{G}}(t,t^{\prime})\equiv\check{\mathcal{G}}(t,t^{\prime})-\underline{\check{\mathcal{G}}}(t-t^{\prime})$, is conveniently represented by

[TABLE]

where

[TABLE]

and the set $\left\{\hat{\Lambda}_{n}[\omega,\omega^{\prime}]\right\}$ is defined by the recurrence relation

[TABLE]

where \hat{G}[\omega]=\big{(}\omega^{2}\,\hat{I}-\hat{K}\big{)}^{-1} has been discussed in the previous section and

[TABLE]

We model the coupling terms as

[TABLE]

where $\phi_{\alpha}(t)$ is a dimensionless function that describes the pumping time-dependence of the $\alpha$-lead. For a periodic pumping, i.e., $\phi_{\alpha}(t+\tau)=\phi_{\alpha}(t)$ the pumping function can be expressed by a Fourier series in harmonic form as

[TABLE]

By construction \big{\langle}\phi(t)\big{\rangle}_{\tau}=0, where $\langle\cdots\rangle_{\tau}\equiv\frac{1}{\tau}\int_{0}^{\tau}\mathrm{d}t\,(\ldots)$ stands for the time average over a period. We assume that $|\phi_{\alpha}(t)|_{\text{m{a}x}}=1$.

Expanding the Dyson equation, Eq. (57), in a power series in $\varepsilon$, we write the energy $E_{M}(t)$ of the extended molecule as

[TABLE]

The explicitly expression for $E_{M}^{(n)}(t)$ are rather lengthy and are given in Appendix C. For a periodic pumping we show that $E^{(0)}_{M}$ does not depend on time and $E_{M}^{(n)}(t)=E_{M}^{(n)}(t+\tau)$ for $n=1,2,\ldots$ (see Appendix C).

We express the variation of the extended molecule energy between $t$ and $t+\Delta t$ in the form of a first law of thermodynamics, namely, $\Delta E_{M}^{(\Delta t)}\equiv\sum_{\alpha}Q_{\alpha}^{(\Delta t)}+W^{(\Delta t)}$. Note that $-Q_{\alpha}^{(\Delta t)}$ corresponds the heat transferred from the molecule to the $\alpha$-reservoir, while $W^{(\Delta t)}$ is the energy transferred to the molecule that does not come from reservoirs, namely,

[TABLE]

where $J_{\alpha}(t)$ and $\Phi(t)$ are, respectively, the thermal current flowing from $\alpha$-reservoir into the molecule and the power developed by the ac sources.

For a periodic process after a cycle of period $\Delta t=\tau$, we finding that $\Delta E_{M}^{(\tau)}=0$, so that

[TABLE]

where we define

[TABLE]

$\left\langle J_{\alpha}(t)\right\rangle_{\tau}$ and $\left\langle\Phi(t)\right\rangle_{\tau}$ can be evaluated by using a perturbative expansion

[TABLE]

where $J_{\alpha}^{(P)}$ and $\Phi^{(P)}$ are discussed in Appendix C and can be cast as

[TABLE]

where the quantities $A_{\beta\gamma}^{\alpha}(n)$, $B_{\beta\gamma}^{\alpha}(n)$, $D_{\beta\gamma}(n)$, and $E_{\beta\gamma}(n)$ are given by intricate expressions involving combinations of equilibrium Green’s functions. The latter are explicitly given by Eq. (114).

Note that in Eq. (68) the first order contributions in $\varepsilon$ vanish. The second order terms $J_{\alpha}^{(P)}$ and $\Phi^{(P)}$ depend explicitly on the periodic profile $\phi_{\alpha}(t)$.

The perturbative approach we put forward allows us to write the pumping currents and power order by order in terms of products and sums of steady-state Green’s functions, which are represented by square matrices of the order of the number of degrees of freedom of the system. Hence, here the numerical bottleneck for addressing realistic systems is the same as in the steady-state, namely, the calculation of the equilibrium Green’s functions as a function of the frequency. Having obtained these objects by any standard method, one needs only to insert the corresponding quantities in the expressions given in the App. C. We note that the non-perturbative regime requires a calculation of the system Green’s functions by directly solving the corresponding differential equations, that is in general a very challenging task.

IV Application: Molecular junction

We investigate the consequences of our findings using the molecular junction model presented in Sec. III. We consider a one-dimensional system where a central region with $N$ atoms is attached to two semi-infinite linear chains acting as leads, as depicted in Fig. 1.

For the sake of clarity, we consider the simplest non trivial case of a diatomic molecule, namely, $N=2$. The force constant between the atoms in the leads and its first neighbors is $k$. The force constant between the atoms in the central region is $k_{C}$ while the left (right) atom connects to the left (right) lead though a coupling $k_{L}$ ($k_{R}$).

In this model the inter-partition and central coupling reduced matrices are

[TABLE]

The retarded and advanced components of the modified Green’s functions are

[TABLE]

for $\alpha=L,R$. Note that the property $\tilde{g}_{\alpha}^{r}[-\omega]=\tilde{g}_{\alpha}^{a}[\omega]$ is satisfied according to the Eq. (20a). The derivation of Eqs. (71) is presented in App. C.

IV.1 Steady-state

Equations (70) and (71), allow us to calculate the retarded and advanced self-energies $\tilde{\Sigma}_{L(R)}^{r,a}[\omega]$ defined in Eq. (46), the level-width functions $\tilde{\Gamma}_{L(R)}[\omega]$ given by Eqs. (51a) and (96), and the central region Green’s functions $G_{CC}^{r,a}[\omega]$. The local density of states (LDOS) at the site $j=A,B$ in the central region reads

[TABLE]

The factor $2\omega$ is present to convert the value coming directly from the imaginary part of $G^{r}_{CC}[\omega]$ into the DOS per unit of $\omega$, ensuring that $\int\text{DOS}(\omega)\,\mathrm{d}\omega$ equals the number of propagating channels in the system.

For the equal force constant case we can calculate the LDOS and the transmission analytically, namely

[TABLE]

Figure 2 shows the DOS at one of the sites in the central region for $k_{L}=k_{R}=k_{C}=k$. Our formalism recovers the standard DOS for a linear chain. The singularity at $\omega=\sqrt{4k}$ agrees with the frequency in which the dispersion relation of a linear chain $\omega=\sqrt{4k}\sin(k_{x}a/2)$ becomes flat, i.e., at the edge of the first Brillouin zone. Here $k_{x}$ is the longitudinal momentum and $a$ is the lattice parameter. Also, the transmission coefficient $\mathcal{T}(\omega)$ corresponds to a perfect transmission inside the frequency band of the leads $|\omega|<\sqrt{4k}$ and it is zero otherwise.

In the limit of small temperatures and small temperature differences, namely, $T_{L/R}=T\pm\Delta T/2$ with $\Delta T\ll T$ for $T\to 0$, the thermal current for steady-state can be written as $J^{(S)}_{L,R}=\pm\,\sigma(T)\;\Delta T$, where $\sigma(T)$ corresponds to the thermal conductance defined by

[TABLE]

where $\omega_{c}\equiv\sqrt{4k}$. From Eq. (55) it is possible to verify that $\mathcal{T}(T\rightarrow 0)=1$.

The low temperature limit of $\sigma(T)$ is

[TABLE]

as theoretically predicted [39, 40] and experimentally observed [21].

Thus, at low temperatures the thermal conductance $\sigma(T)\propto T$

vanishes for $T\to 0$, as required by the third law of thermodynamics.

Let us now study situations where the force constants are different. In the weak coupling limit, $k_{R},k_{L}\ll k,k_{C}$, the central region is nearly disconnected from the outside world having only one resonant level at $\omega_{C}=\sqrt{2k_{C}}$. Thus, the conductance is only expected to be significant at the vicinity of $\omega_{C}$. Instead, Fig. 3 shows one peak at $\omega\approx\omega_{C}$ and two additional strong peaks, one at zero frequency and another intermediate peak at $0<\omega<\omega_{C}$.

The first peak at $\omega=0$ corresponds to the acoustic mode that has an infinite long wavelength so that the short ranged “defects” introduced by $k_{L},k_{R},k_{C}\neq k$ do not affect the transport across the system. This picture is reenforced by noticing that the zero frequency peak is robust against changes in the value of $k_{L}$, $k_{C}$ and $k_{R}$ in the weak coupling regime, see Fig. 4.

By coupling the diatomic molecule to leads, the resonance level at $\omega_{C}$ is shifted and acquires broadening, as described by the self-energy $\tilde{\Sigma}^{r}[\omega]$. Hence the peak near $\omega_{C}$ is very sensitive to variations in $k_{L}$ and $k_{R}$. These features, are illustrated in Fig. 3, by inspecting a set of transmission curves where we keep $k_{C}$ constant and increase $k_{L}=k_{R}$.

On the other hand, the remaining peak at $0<\omega<\omega_{C}$ depends only on the values of $k_{L}$ and $k_{R}$. In the weak coupling regime, a semi-classical picture explains this additional transmission peak. The natural interfaces frequencies $\omega_{\alpha}\propto\sqrt{k_{\alpha}}$, with $\alpha=L,R$, are much smaller then $\omega_{C}$. The large separation in frequencies suggest that the resonance close to $\omega_{C}$ is dominated by the isolated molecule mode, while the other corresponds to an oscillation of a frozen central region. The Green’s functions of such a system gives resonances at the frequencies

[TABLE]

For $k_{L}=k_{R}$, $\omega_{1}=\sqrt{2k_{C}+k_{L}}$ and $\omega_{2}=\sqrt{k_{L}}$ that are plotted in Fig. 3 as vertical dotted lines matching the peaks positions. For $k_{L}\neq k_{R}$, the symmetry is broken and the maximum transmission at all the peaks, except for the one with zero frequency, is no longer perfect.

We note that our results are qualitatively similar to those in Ref. [36], that analyze the steady-state transport through a benzene ring. There is an important difference though: Taking into account $V_{aa}$ in $\tilde{g}_{\alpha}$ guarantees that ${\cal T}(\omega\rightarrow 0)\rightarrow 1$, which is a necessary condition to obtain the quantum of thermal conductance for $T\rightarrow 0$. In distinction, by using $g_{\alpha}$ as the surface Green’s function, as done in Refs. [36, 37, 38], one obtains ${\cal T}(\omega\rightarrow 0)\rightarrow 0$.

IV.2 Pumping

For simplicity, let us analyze a pumping process between reservoirs at the same temperature. In this case, the steady-state current from the $\alpha$-reservoir is $J_{\alpha}^{(S)}=0$. Hence, $J_{\alpha}^{(P)}$ gives the leading contribution to the heat flow.

We consider the case of pumping functions with a phase difference $\varphi$, namely, $\phi_{L}(t)=\phi_{R}(t-\varphi/\Omega)$, which implies that $\varphi_{n}^{(L)}-\varphi_{n}^{(R)}=n\,\varphi$ and $a_{n}^{(L)}=a_{n}^{(R)}\equiv a_{n}$ for $n\geqslant 1$. According to Eq. (69), we can express the $\alpha$-thermal pumped current as

[TABLE]

where $\mathcal{A}^{\alpha}_{\text{homo}}(n\Omega)\equiv A^{\alpha}_{LL}(n)+A^{\alpha}_{RR}(n)$, $\mathcal{A}^{\alpha}_{\text{hete}}(n\Omega)\equiv A^{\alpha}_{LR}(n)+A^{\alpha}_{RL}(n)$, $\mathcal{B}^{(\alpha)}(n\Omega)\equiv B^{\alpha}_{LR}(n)-B^{\alpha}_{RL}(n)$ for $\alpha=L,R$. For the symmetric coupling case, i.e., $k_{L}=k_{R}\neq k_{C}$, we can show that $\mathcal{A}^{L}_{\text{homo/hete}}(n\Omega)=\mathcal{A}^{R}_{\text{homo/hete}}(n\Omega)$ and $\mathcal{B}^{L}(n\Omega)=-\mathcal{B}^{R}(n\Omega)$.

Similarly, the pumped power reads

[TABLE]

where $\mathcal{D}_{\text{homo}}(n\Omega)\equiv D_{LL}(n)+D_{RR}(n)$, $\mathcal{D}_{\text{hete}}(n\Omega)\equiv D_{LR}(n)+D_{RL}(n)$ and $\mathcal{E}(n\Omega)\equiv E_{LR}(n)-E_{RL}(n)$ defined for $0<n\,\Omega<4\sqrt{k}$ and zero otherwise. For further details see Appendix C. For a symmetric setup (i.e. $k_{L}=k_{R}\neq k_{C}$), we can show that $\mathcal{E}(n\Omega)\equiv 0$. Note that $\Phi^{(P)}$ satisfies the condition $\Phi^{(P)}>0$,

as exemplified for the diatomic molecule in the Fig. 5, which corresponds to positive rate of work performed on the system. Therefore, for reservoirs at the same temperature, the entropy production per cycle cast as $(\Delta S)_{\text{cycle}}/\tau=\varepsilon^{2}\,\dot{\mathcal{S}}^{(P)}+\mathcal{O}(\varepsilon^{4})$ satisfies

[TABLE]

as expected from the second law of thermodynamics. Note that the overall partition scheme, the definitions of the heat currents, and the pumped power are consistent with general thermodynamic properties.

The pumping function is determined by the choice of the parameter set $\{a_{n}\}$ and phase difference $\varphi$. We study four examples of pumping functions: single-mode represented by $a_{n}=\delta_{n,1}/2$; square oscillation represented by $a_{n}=2/(n\pi)$ for $n$ odd (and zero otherwise); triangle oscillation by $a_{n}=\frac{4}{n^{2}\pi^{2}}(-1)^{\frac{n-1}{2}}$ for $n$ odd (and zero otherwise); sawtooth oscillation by $a_{n}=-1/(\pi n)$ for all $n$.

The thermal current absorbed by the $\alpha$-reservoir, $-J_{\alpha}^{(P)}>0$, as a function of pumping frequency $\Omega$ for different pumping profiles and phase difference is shown in Fig. 5. We find a suppression of thermal current according to the type of pumping in the following order: square, single-mode, triangle and sawtooth. The pumping peak occurs in the frequency window $2\sqrt{k}<\Omega<3\sqrt{k}$ and a sub-peak within $0<\Omega<\sqrt{k}$, accompanied by a weak suppression in the domain $\sqrt{k}<\Omega<2\sqrt{k}$ and the strong suppression for $\Omega\gtrsim 4\sqrt{k}$. As the phase difference is increased, the suppression in $\sqrt{k}<\Omega<2\sqrt{k}$ is intensified. Note that the unperturbed proposed setup is symmetric ($k_{L}=k_{R}\neq k_{C}$). The phase shift $\varphi$ causes the difference between $J_{L}^{(P)}$ and $J_{R}^{(P)}$, see Fig. 5.

The effective heat flux between the two reservoirs $\Delta J^{(P)}(\Omega)\equiv J^{(P)}_{R}(\Omega)-J^{(P)}_{L}(\Omega)$ for the symmetric coupling case (i.e. $k_{L}=k_{R}\neq k_{C}$) reads

[TABLE]

where $\Delta J^{(P)}>0$ corresponds to the heat flux from right-reservoir to left-reservoir and $\Delta J^{(P)}<0$ to reverse direction. Note that for $\varphi=0,\pm 2\pi,\pm 4\pi,\pm 6\pi,\ldots$, we obtain $\Delta J^{(P)}=0$, as expected. For the single-mode case we verify that the maximum value of $|\Delta J^{(P)}|$ occurs for $\varphi=\pm\pi/2,\pm 3\pi/2,\pm 5\pi/2,\ldots$. $\Delta J^{(P)}$ versus the pumping frequency $\Omega$ is represented in the Figs. 5(g)-(i).

The time-dependent drive breaks the system symmetry. Hence, one can engineer configurations of $\Omega$ and $\phi_{\alpha}(t)$ that direct the pumped heat either to the left or to the right. Equation (80) shows that by replacing $\varphi\rightarrow-\varphi$ the directionality is reversed, namely, $\Delta J^{(P)}(\varphi)=-\Delta J^{(P)}(-\varphi)$.

For all considered pumping profiles and phase differences, we find $\Delta J^{(P)}=0$ at $\Omega=\Omega_{1}\approx 1.371\sqrt{k}$ and $\Omega=\Omega_{2}\approx 2.029\sqrt{k}$, see Fig. 5. The largest negative peak of $\Delta J^{(P)}$ occurs between $\Omega_{1}$ and $\Omega_{2}$. Other two positive peaks occur on $\Omega<\Omega_{1}$ and $\Omega>\Omega_{2}$.

The external drives contribute to the energy transfer by exciting the original unperturbed propagating energy $\hbar\omega$ to $\hbar(\omega+n\Omega)$ where $n=1,2,3,\ldots$. In our model, the energy transfer between leads is only possible if the injected energy $\hbar\omega$ and the excited energy $\hbar(\omega+n\Omega)$ satisfy the conditions $|\omega|\leq\omega_{c}$ and $|\omega+n\Omega|\leq\omega_{c}$, respectively, where $\omega_{c}\equiv 2\sqrt{k}$. Outside this frequency window the leads linewidths vanish and no energy transfer is allowed. Thus, only modes with $n\leq 2\omega_{c}/\Omega$ will contribute to the energy transfer between reservoirs. As $\Omega$ increases a smaller number of modes contribute and the overall energy transport decreases, as we see in Fig. 5. For $\Omega>2\omega_{c}=4\sqrt{k}$ the external drives do not induce energy transport to the reservoirs as no positive integer can satisfy $n<1$. For small $\Omega$, many modes $n$ with $n\leq 2\omega_{c}/\Omega$ compete and contribute to the energy transfer.

From Eqs. (67) and (68), we define the (cooling) efficiency $\kappa$ of the heat pump that operates between two reservoirs left ($L$) and right ($R$) at equal temperatures, considering a full period as

[TABLE]

$|\kappa|$ is the ratio between the net current and the total heat current driven by the external ac source per cycle. Figure 6 shows that $\kappa$ can be positive or negative depending on $\Omega$ and $\phi_{L,R}(t)$. $\kappa=\pm 1$ correspond to situations where the heat currents have opposite signs, independent of their magnitudes.

For $\Omega\gtrsim 0.5\sqrt{k}$, the thermal energy is always absorbed by the reservoirs for all the studied pumping profiles (see Fig. 5). Thus, the denominator of Eq. (81), $|Q_{R}^{(\tau)}|+|Q_{L}^{(\tau)}|$, corresponds to realized work $W>0$. In contrast, when $\Omega\lesssim 0.5\sqrt{k}$ for the single-mode and the triangle profiles, we find a plateau $\kappa=-1$, indicated in the Fig.6. In this cases, the external drive pumps heat from the left ($-Q_{L}^{(\tau)}<0$) to the right reservoir ($-Q_{R}^{(\tau)}>0$) and $\kappa=-1$ indicates all the thermal energy extracted from the left reservoir by the external drive is transferred to the right one. Unfortunatly, the corresponding heat current is rather small. Outside of this regime, we find that in order to optimize the pumped heat $\Delta J^{(P)}$ one should: (i) tune the pumping frequency to $\Omega\approx\sqrt{k}$ or $\Omega\approx 1.75\sqrt{k}$ producing positive ($\Delta J^{(P)}>0$) or negative ($\Delta J^{(P)}<0$), respectively, heat transfer (see Fig. 5) and (ii) tune the phase shift to $\pi/2$ (see Eq. (80)).

The high frequency, $\Omega\gg\sqrt{k}$, asymptotic behavior of $\kappa$ (see Fig. 6) is trivial due to the contribution of a smaller number of modes $n$ as discussed before. We also note a surprisingly similar behavior of the efficiency curves for single-mode and triangle pumping profiles. This can be explained by inspecting Eq. (80) and recalling the rapid decrease of $a_{n}^{2}\propto 1/n^{4}$ with $n$ for the latter case. In other words, the pumped current, Eq. (77), for the triangle profile is dominated by the first mode $n=1$, resulting in a frequency dependence similar to the single-mode profile.

V Conclusions

In this paper we have presented a rigorous description of quantum thermal transport properties due to phonons in molecular and nanomechanical systems using the non-equilibrium Green’s function theory. We approached the problem using a phase-space representation on the quantum correlations functions in the Keldysh contour, a convenient generalization of the standard Green’s functions technique [29, 30].

We have shown that in the stationary regime our approach recovers a Landauer-like transmission formula, as expected. Our derivation solves some inconsistencies of previous theoretical works based on NEGF [36, 37, 38]. For instance, the use of phase-space correlation functions avoids the necessity of taking the Fourier transform of $\vec{u}$ and canonically conjugate $\vec{p}$ operators, that is troublesome for the commutation relations $[\vec{u},\vec{p}]$. The partition we put forward in Sec. III, avoids conceptual difficulties with the adiabatic switch-on picture on picture on which the formalism is based.

Finally, starting with a symmetrized Hamiltonian, our formalism avoids the necessity of imposing the ad hoc symetrization $J=\left(J_{L}+J^{\ast}_{L}-J_{R}-J^{\ast}_{R}\right)/4$ used in Refs. [36, 37, 38].

We extend the formalism to study the heat transport in systems subjected to a time-dependent external drive which opens the possibility of addressing situations of interest for applications in phononics. In distinction to the electronic case where the Fermi velocity and the system size give a characteristic time scale for the dynamics, the absence of such time scale in bosonic systems leads us to develop a new perturbation theory scheme assuming that the external drive is weak.

We apply our results to a model of a diatomic molecule coupled to semi-infinite linear chains in equilibrium with thermal reservoirs. The simplicity of the model allows for an amenable computation and to understand its main physical features using simple analytical considerations. This gives us confidence on the method and we expect it to be used to treat more realistic systems.

Acknowledgements.

This work is supported by the Brazilian funding agencies CAPES, CNPq, FAPERJ and FAPEAL. The authors thank the hospitality of the International Institute of Physics (IIP) in Natal (Brazil), where this work was concluded.

Appendix A Canonical commutation relations and Fourier transform in frequency space

The canonical quantization procedure of a classical Hamiltonian expressed in terms of the set of independent variables represented by the displacements $\{u_{i}\}$ and by canonically conjugated momenta $\{p_{j}\}$, renders the commutation relations $\left[u_{i}(t),p_{j}(t)\right]=\imath\hbar\,\delta_{i,j}$ and $\left[u_{i}(t),u_{j}(t)\right]=0=\left[p_{i}(t),p_{j}(t)\right]$.

The standard approach [36] is not consistent with the above relations, as we show below.

Let us consider the Fourier transform of $u_{i}(t)$ and $p_{i}(t)\equiv\dot{u}_{i}(t)$ as [36, 37, 38]

[TABLE]

where the condition $(u_{i}[-\omega])^{\dagger}=u_{i}[\omega]$ must be satisfied as a result of $\left(u_{i}(t)\right)^{\dagger}=u_{i}(t)$ (and reciprocally $\left(p_{i}(t)\right)^{\dagger}=p_{i}(t)$).

Hence, the canonical commutation relations $\left[u_{i}(t),p_{j}(t)\right]=\imath\hbar\,\delta_{i,j}$ and $\left[u_{i}(t),u_{j}(t)\right]=0$ can be written as

[TABLE]

Note that it is not possible to obtain a consistent result for $\left[u_{i}[\omega],u_{j}[\omega^{\prime}]\right]$ with the (83a) e (83b), simultaneously. This results from the fact that canonicity relations involve operators at equal times, which is not compatible with the transformation (82b). Thus, the frequency Fourier transform (82b) is not a canonical transformation.

In our construction we obtain the equations of motion for the NEGFs. The Fourier transform to frequency space is performed in the Green’s functions arguments and not in the displacement and momentum operators, as standard. Hence, circumventing potential problems with the commutation relations.

Appendix B Surface Green’s functions for semi-infinite lattices

In this Appendix we present a novel direct analytical calculation of $\tilde{g}_{\alpha}^{r,a}[\omega]$ for the non-ideal coupling case, namely, $k_{\alpha}\neq k$. (The results of Refs. [36, 37] are recovered by taking $k_{\alpha}=k$.) Next, we discuss the importance of the term $V_{aa}$ which originates the difference between ${g}_{\alpha}^{r,a}[\omega]$ and $\tilde{g}_{\alpha}^{r,a}[\omega]$.

According to (47) and (48c), we write

[TABLE]

where $|e_{1}\rangle=(1,0,\cdots,0)^{\text{T}}$ represents the surface site and

[TABLE]

with $\omega_{\pm}\equiv\omega\pm\imath\,0^{+}$ and $n\to\infty$.

Applying the method of co-factors in (84) we write

[TABLE]

where $d_{n}^{\pm}\equiv(-1)^{n}\det\mathcal{D}_{n}^{\pm}$. The Laplace’s method gives the following recurrence equation

[TABLE]

The discriminant $\Delta=\omega_{\pm}^{2}\,\left(\omega_{\pm}^{2}-4k\right)$ of the associated characteristic equation has non trivial roots $|\omega|=\sqrt{4k}$. Hence, we split the solution of the recurrence equation in two frequency domains: (i) $|\omega|\leqslant\sqrt{4k}\,$ and (ii) $|\omega|>\sqrt{4k}$.

(i) For $|\omega|\leqslant\sqrt{4k}$ we introduce the parametrization $\omega_{\pm}\equiv\omega\pm\imath\,0^{+}=\sqrt{4k}\,\sin\left(\theta_{\pm}/2\right)$ with $\theta_{\pm}=\theta\pm\imath\eta$ and $\eta=0^{+}$ for $\theta\in[-\pi,\pi]$. Substituting the latter in (87), we find

[TABLE]

Since $\tan\left[m(\theta\pm\imath\,\eta)\right]\sim\pm\imath$ for $m\gg 1$ and $\eta>0$,

[TABLE]

that, with the help of the identities $2k\cos\theta=2k-\omega^{2}$ and $2k\sin\theta=\omega\sqrt{4k-\omega^{2}}$, leads to

[TABLE]

Note that (90) satisfies the property $\tilde{g}_{\alpha}^{r}[-\omega]=\tilde{g}_{\alpha}^{a}[\omega]$ in line with (20a).

(ii) For $|\omega|>\sqrt{4k}$, we parametrize $\omega_{\pm}\equiv\omega\pm\imath\,0^{+}=\sqrt{4k}\,\cosh\left(\theta_{\pm}/2\right)\,\text{sgn}(\theta)$ with $\theta_{\pm}=\theta\pm\imath\eta$ and $\eta=0^{+}$ for $\theta\in\mathbb{R}^{\ast}$, where $\text{sgn}(\theta)$ is the sign function of $\theta$. Substituting this parametrization in (87) we find

[TABLE]

Since $\coth\left[m(\theta\pm\imath\,\eta)\right]\sim\text{sgn}(\theta)$ for $m\gg 1$ and $\eta>0$,

[TABLE]

Using $2k\,\cosh\theta=\omega^{2}-2k$ and $2k\,\text{sgn}(\theta)\,\sin\theta=\sqrt{\omega^{2}(\omega^{2}-4k)}$, we obtain

[TABLE]

Note that $\tilde{g}^{r,a}_{\alpha}[\omega]\sim 1/\omega^{2}\to 0$ for $|\omega|\gg\sqrt{4k}$, which guarantees convergence in the integrations.

We can write $\tilde{g}^{r,a}_{\alpha}[\omega]$ in a convenient form as

[TABLE]

where the real auxiliary functions $\tilde{\gamma}_{\alpha}[\omega]$ and $\tilde{\mu}_{\alpha}[\omega]$ are

[TABLE]

It is straightforward to verify that $\tilde{g}_{\alpha}^{r}[\omega]$ satisfies the Kramers-Kronig relations, as it should [59, 60].

The $\alpha$-contact line width function, $\tilde{\Gamma}_{\alpha}[\omega]$, Eq. (51a), becomes

[TABLE]

Let us now analyze the role of term $V_{aa}$ in the transmission $\mathcal{T}(\omega\rightarrow 0)$, given by Eq. (55). We consider a system where the central region is composed by a dimer as shown in Fig. 1. According to Eqs. (46), (70) and (90) the low frequency limit of $\tilde{\Gamma}_{\alpha}[\omega]$ and $\tilde{\Sigma}^{r,a}[\omega]$ are

[TABLE]

where $P_{L}=\begin{pmatrix}1\;&\;0\\ 0\;&\;0\end{pmatrix}$, $P_{R}=\begin{pmatrix}0\;&\;0\\ 0\;&\;1\end{pmatrix}$, $V_{CC}=\sum_{\alpha}k_{\alpha}\,P_{\alpha}$ and $\alpha=L,R$. The central region Green’s function is

[TABLE]

where $I_{2}$ is a $2\times 2$ identity matrix.

Using Eq. (98) in Eq. (99), we obtain $G^{r,a}_{CC}[\omega]=[-K_{CC}^{0}\,\pm\,\imath\,\sqrt{k}\,\omega\,(P_{L}+P_{R})+\mathcal{O}(\omega^{2})]^{-1}$, where the spring-constant matrix of the decoupled central region $K_{CC}^{0}$ is singular and gives rise to the expansion

[TABLE]

Substituting Eqs. (97) and (100) into Eq.(55), we obtain $\mathcal{T}(\omega)=1+\mathcal{O}(\omega)$, that is, $\mathcal{T}(\omega\rightarrow 0)=1$.

In summary, the term $V_{aa}$ in Eq. (70) leads to $\text{Re}\{\tilde{\Sigma}^{r,a}[0]\}=-V_{CC}$. The latter cancels out the term $-V_{CC}$ in Eq. (99), leading to Eq. (100), that results in unit transmission for $\omega\rightarrow 0$. We conclude that for the general non-ideal coupling case, the use of $\tilde{g}^{r,a}_{\alpha}[\omega]$ is key to obtain the correct transmission low-frequency behavior.

Appendix C Perturbative weak pumping regime

In this Appendix, we derive the perturbation expansion for the energy $E_{M}(t)$ of the extended molecule and analyze the periodic behavior for pumped-induced heat transport. Next, we obtain the perturbation expansion in $\varepsilon$ for thermal current $J_{\alpha}(t)$ and the power $\Phi(t)$.

C.1 Energy $E_{M}(t)$

We expand the Dyson equation, Eq. (57), in a power series in $\check{\mathcal{V}}(t)$, to obtain (after a lengthy but straightforward calculation) the energy $E_{M}(t)$ as

[TABLE]

where $\underline{\hat{K}}\equiv\hat{K}^{0}+\hat{V}$ and \big{(}\hat{G}[\omega]\cdot\underline{\hat{K}}\big{)}_{CC}=\sum_{\alpha}\big{(}G_{C\alpha}^{<}[\omega]\cdot V_{\alpha C}+G_{CC}^{<}[\omega]\cdot V_{CC}^{(\alpha)}\big{)}.

Here, $\phi_{\alpha}[\omega]$ is the Fourier’s transform of the pumping function, Eqs. (61) and (63), given by

[TABLE]

$\Xi_{ab,\beta}^{<}[\omega,\omega^{\prime}]$ and $\Xi_{ab,\beta\gamma}^{<}[\omega,\nu,\omega^{\prime}]$ are lesser components of

[TABLE]

respectively, with latin letters corresponding to reservoirs or $C$ and greek letters corresponding to reservoirs only.

Note that $\Lambda_{1}^{<}$ and $\Lambda_{2}^{<}$ of Eq. (60) are related to $\Xi_{ab,\beta}^{<}$ and $\Xi_{ab,\beta\gamma}^{<}$ by

[TABLE]

Note that $E_{M}^{(0)}$ is constant. Substituting (C.1) in Eqs. (101b) and (101c), we can see that $E_{M}^{(n)}(t)=E_{M}^{(n)}(t+\tau)$ for $n=1,2,\ldots$.

C.2 Current $J_{\alpha}(t)$ and power developed by the ac sources $\Phi(t)$

Substituting the results of Eqs. (57)-(62) into (40)-(41), we obtain the current $J_{\alpha}(t)$ from $\alpha$-lead and the power developed by the ac sources $\Phi(t)$ in the form of a perturbative series in $\varepsilon$ as

[TABLE]

where $J^{(n)}_{\alpha}(t)$ and $\Phi^{(n)}(t)$ are $n$-order contribution of the series of $J_{\alpha}(t)$ and $\Phi(t)$, respectively. Below we give the explicit expressions for the first and second-order contributions.

C.2.1 First-order contribution

The coefficients $J_{\alpha}^{(1)}(t)$ and $\Phi^{(1)}(t)$ read

[TABLE]

Using Eqs. (61) and (C.1), we find

[TABLE]

where $\langle\cdots\rangle_{\tau}=\frac{1}{\tau}\int_{0}^{\tau}\mathrm{d}t\,(\cdots)$. Thus,

[TABLE]

C.2.2 Second-order contribution

The coefficients $J_{\alpha}^{(2)}(t)$ and $\Phi^{(2)}(t)$ read

[TABLE]

Using the Eqs. (63) and (C.1), we obtain

[TABLE]

Hence,

[TABLE]

where we introduced the following integrals

[TABLE]

Defining J^{(P)}_{\alpha}\equiv\big{\langle}J_{\alpha}^{(2)}(t)\big{\rangle}_{\tau} and \Phi^{(P)}\equiv\big{\langle}\Phi^{(2)}(t)\big{\rangle}_{T}, we get

[TABLE]

where

[TABLE]

Equations (113) and (114) and the energy conservation $\sum_{\alpha}J_{\alpha}^{(P)}+\Phi^{(P)}=0$ (according Eqs. (66)-(68)), lead to the following conditions

[TABLE]

where we introduced symmetrization $O_{(\beta\gamma)}\equiv\frac{1}{2}\left(O_{\beta\gamma}+O_{\gamma\beta}\right)$ and anti-symmetrization $O_{[\beta\gamma]}\equiv\frac{1}{2}\left(O_{\beta\gamma}-O_{\gamma\beta}\right)$ shorthand notations.

For the calculation of $\mathcal{J}_{\beta\gamma}^{(\alpha)}(\sigma\Omega_{n})$ and $\mathcal{F}_{\alpha\beta}(\sigma\Omega_{n})$, we use $V_{C\alpha}\cdot\Xi_{\alpha C,\beta}[\omega,\omega^{\prime}]$, $\Xi_{CC,\beta}[\omega,\omega^{\prime}]$ and $V_{C\alpha}\cdot\Xi_{\alpha C,\beta\gamma}[\omega,\omega^{\prime},\omega]$ of (45) and (103), as

[TABLE]

where we define

[TABLE]

where

[TABLE]

[TABLE]

where we define the lesser component of the set (117) as

[TABLE]

[TABLE]

and the retarded-($r$) and advanced-($a$) component of (117) as

[TABLE]

with $x=r,a$ and where we denote the components of generalized function as

[TABLE]