*3.1. Model*

Let us consider the energy transfer via a two-level system between two environmental systems (*L* and *R*) consisting of an infinite number of bosons [28,54–56]. With the definition of the lower (higher) level of the two-level system as |0 (|1), respectively, the Hamiltonian of the model is written as

$$H\_{\mathbb{S}\nu} = \sum\_{m=0,1} \varepsilon\_{m} |m\rangle\langle m|\_{\prime} \quad H\_{\mathbb{E}} = \sum\_{k,\nu} \hbar \omega\_{k,\nu} b\_{k,\nu}^{\dagger} b\_{k,\nu}$$

$$H\_{\mathbb{S}\nu} = \sum\_{\mathbb{V}} X\_{\mathbb{V}} (|0\rangle\langle 1| + |1\rangle\langle 0|) , \qquad (\nu = L, R) , \tag{15}$$

where *Xν* = ∑*k hg*¯ *<sup>k</sup>*,*<sup>ν</sup>*(*b*†*k*,*<sup>ν</sup>* + *bk*,*<sup>ν</sup>*), with *b*†*k*,*ν* and *bk*,*<sup>ν</sup>* the creation and annihilation boson operators of the *k*th mode of the *ν*th environment. The model scheme is shown in Figure 1A.

We consider the energy transfer from out-of-phase temperature modulations of the two environments, corresponding to the bias averaged out during a period T [28]. To discuss the nonadiabaticity for this model, we study energy (boson) transfers under cyclic and piecewise modulations of the environmental temperatures *TL* and *TR* dividing T into *N* intervals, *ti* ≤ *t* ≤ *ti*+<sup>1</sup> (*i* = 1, ··· , *N*) with *t*1 = 0 and *tN*+<sup>1</sup> = T . We need to discretize the temperature modulation because conventional treatments describing relaxation phenomena require the environmental temperature to remain constant. By changing the number of intervals of the temperature modulation, we compared each time interval with the relaxation time of the relevant two-level system and thus we are able to discuss nonadiabaticity explicitly, for example, from the scale between *τ*<sup>−</sup><sup>1</sup> *r* and Ω. In taking the limit *N* → <sup>∞</sup>, we reveal energy transfer features under a continuous modulation. In Figure 1B, we plot the time dependence of the temperature modulations used in this study calculated for a typical number of discrete time intervals *N* = 20.

**Figure 1.** (**A**) Model scheme: a two-level system as an anharmonic junction interacts with two environments (*L* and *R*). (**B**) Temperature modulations, *TL*(*t*) = 200 + 100 cos(*ω<sup>t</sup>* + *<sup>π</sup>*/4), and *TR*(*t*) = 200 + 100 sin(*ω<sup>t</sup>* + *<sup>π</sup>*/4), discretized with *N* = 20.

### *3.2. FCS Formalism Applied to Pumping of Energy Quanta*

We apply the general formalism of the FCS in the former section to this model focusing on weak system–environment coupling and considering long time (Born-Markovian) limits by taking *t* → ∞ in the upper bound of the integral in Equation (10). In this limit, the super-operator *ξ*(*λ*)(*t*) becomes time independent during each interval. We find that the mean value of the transferred quantity between the relevant two-level system and the *ν*th environment in the *i*th time interval, <sup>Δ</sup>*q<sup>ν</sup>i* is written as

$$
\langle \Delta q\_i^{\nu} \rangle = \hbar \omega\_0 \{ A\_i^{\nu} \int\_{t\_i}^{t\_{i+1}} dt' \rho\_{00}(t') - B\_i^{\nu} \delta t \},
\tag{16}
$$

with *ω*0 = (*<sup>ε</sup>*1 − *<sup>ε</sup>*0)/¯*h*, *ρ*00(*t*) = 0|*ρ*(*<sup>λ</sup>*=<sup>0</sup>)(*t*)|0, and *δt* = *ti*+<sup>1</sup> − *ti*. *<sup>A</sup><sup>ν</sup>i* and *Bνi* in Equation (16) are coefficients defined as *<sup>A</sup><sup>ν</sup>i* = <sup>−</sup>(*k<sup>ν</sup>d*,*<sup>i</sup>* + *<sup>k</sup><sup>ν</sup>u*,*<sup>i</sup>*) and *Bνi* = <sup>−</sup>*k<sup>ν</sup>u*,*<sup>i</sup>* with rate constants *<sup>k</sup><sup>ν</sup>u*,*<sup>i</sup>* and *<sup>k</sup><sup>ν</sup>d*,*i*, which govern the time evolution of *ρ*00(*t*) during *ti*−<sup>1</sup> ≤ *t* < *ti*. Their explicit expressions are *<sup>k</sup><sup>ν</sup>d*,*<sup>i</sup>* = <sup>Γ</sup>*νnν*,*<sup>i</sup>* and *<sup>k</sup><sup>ν</sup>u*,*<sup>i</sup>* = <sup>Γ</sup>*ν*(<sup>1</sup> + *<sup>n</sup>ν*,*<sup>i</sup>*), where *<sup>n</sup>ν*,*<sup>i</sup>* = 1/(exp[*h*¯ *βνi <sup>ω</sup>*0] − 1) is the Bose–Einstein distribution for the inverse temperature *βνi* of the *ν*th environment during the *i*th interval and Γ*ν* denotes the feature of the system–environment coupling as Γ*ν* = <sup>2</sup>*πhν*(*<sup>ω</sup>*0) with the coupling spectral density *hν*(*ω*) ≡ ∑*k <sup>g</sup>*<sup>2</sup>*k*,*<sup>ν</sup><sup>δ</sup>*(*<sup>ω</sup>* − *<sup>ω</sup>k*,*<sup>ν</sup>*) = *λω* exp[−*<sup>ω</sup>*/*<sup>ω</sup>c*], where *λ* is the coupling strength and *ωc* is the cutoff frequency. To obtain <sup>Δ</sup>*q<sup>ν</sup>i* , we need *ρ*00(*t*), the time evolution for which is

$$\phi\_{00}(t) = -K\_{d,i}\rho\_{00}(t) + K\_{\mu,j}\rho\_{11}(t)$$

with *Kd*,*<sup>i</sup>* = ∑*ν <sup>k</sup><sup>ν</sup>d*,*<sup>i</sup>* and *Ku*,*<sup>i</sup>* = ∑*ν <sup>k</sup><sup>ν</sup>u*,*i*. The differential equation for *ρ*00(*t*) is solved to give

$$
\rho\_{00}(t) = \rho\_{s,i} + e^{\Lambda\_i t} (\rho\_{00}(t\_{i-1}) - \rho\_{s,i}),
\tag{17}
$$

where we denote *ρ<sup>s</sup>*,*<sup>i</sup>* = <sup>−</sup>*Ku*,*<sup>i</sup>*/Λ*<sup>i</sup>* with Λ*i* = <sup>−</sup>(*Kd*,*<sup>i</sup>* + *Ku*,*<sup>i</sup>*). Using Equations (12), we find that the total transferred energy during the period is calculated to be

$$
\langle \Delta q^{\upsilon} \rangle = \sum\_{i=1}^{N+1} \langle \Delta q\_i^{\upsilon} \rangle = \hbar \omega\_0 (\mathcal{G}\_1^{\upsilon} + \mathcal{G}\_2^{\upsilon} + \mathcal{G}\_3^{\upsilon}),
\tag{18}
$$

where

$$\mathcal{G}\_1^v \quad = \sum\_{i=1}^{N+1} (A\_i^v \rho\_{s,i} - B\_i^v) \delta t\_\prime \tag{19}$$

$$\mathcal{G}\_2^{\vee} = \sum\_{i=1}^{N} \frac{A\_{i+1}^{\vee}}{\Lambda\_{i+1}} (\rho\_{s,i+1} - \rho\_{s,i}),\tag{20}$$

$$\mathcal{G}\_3^{\mathbb{V}} = \sum\_{i=1}^{N+1} \Phi\_{0j}^{\mathbb{V}} + \sum\_{i=2}^{N} (\rho\_{s,i-1} - \rho\_{s,i})\Psi\_i^{\mathbb{V}} + (\rho\_{s,n-1} - \rho\_{s,n})\frac{A\_n^{\mathbb{V}}}{\Lambda\_n} \varepsilon^{\Lambda\_n \delta t},\tag{21}$$

with

$$\boldsymbol{\phi}\_{0,i}^{\boldsymbol{\nu}} = \ (\boldsymbol{\rho}\_{00}(0) - \boldsymbol{\rho}\_{s,1}) \boldsymbol{f}^{\boldsymbol{\nu}}(1, i), \tag{22}$$

$$\Psi\_i^{\nu} = \frac{A\_i^{\nu}}{\Lambda\_i} e^{\Lambda\_i \delta t} + \sum\_{m=i}^{N} f^{\nu}(i, m+1), \tag{23}$$

$$f^{\nu}(p,q) \quad = \frac{A\_q^{\nu}}{\Lambda\_q} \epsilon^{\sum\_{v=p}^{q-1} \Lambda\_v \delta t} (\epsilon^{\Lambda\_q \delta t} - 1). \tag{24}$$

In the next subsection, we show the physical meanings of these obtained terms.
