Minimal quantum thermal machine in a bandgap environment: non-Markovian features and anti-Zeno advantage

A minimal model of a quantum thermal machine is analyzed, where a driven two level working medium (WM) is embedded in an environment (reservoir) whose spectrum possesses bandgaps. overlap with hot or cold reservoirs whose spectra are separated by a bandgap. Approximate and exact treatments supported by analytical considerations yield a complete characterization of this thermal machine in the deep quantum domain. For slow to moderate modulation, the spectral response of the reservoirs is close to equilibrium, exhibiting sideband (Floquet) resonances in the heat currents and power output. In contrast, for faster modulation, strong-coupling and non-Markovian features give rise to correlations between the WM and the reservoirs and between the two reservoirs. Power boost of strictly quantum origin ('quantum advantage') is then found for both continuous and segmental fast modulation that leads to the anti-Zeno effect of enhanced spectral reservoir response. Such features cannot be captured by standard Markovian treatments.


Introduction
The description of macroscopic, classical thermal machines is based on the paradigm that their dynamics is represented by an adiabatic sequence of equilibrium states wherein the system and the reservoir are separable [1]. In recent years, it has turned out that in the quantum domain this paradigm does not commonly lead to heat-machine operation that is distinctly quantum mechanical [2,3]. To describe non-equilibrium thermal machines in the quantum domain, the leading approaches [2][3][4][5][6] have resorted to Born-Markov master equations, often under simplified assumptions: for example, individual strokes in an Otto cycle have been treated independently and then "glued" together to a full cycle. This in turn requires complete system-reservoir thermalization and separability in the appropriate strokes, thus remaining within the foregoing paradigm of classical thermodynamics.
A broader dynamical description of thermal machines that is fully reconciled with quantum mechanics should, however, account for the non-separability and correlations of the reservoirs with the system (the working medium -WM): these features are generic at low temperatures [7] and inevitably affect the machine operation [8]. Their treatment requires us to venture beyond the Born-Markov approximation [9]. Considerable progress in this direction has been recently achieved by revealing non-Markovian and correlation effects in heat machines [10][11][12][13][14][15]. In particular, it has been shown [16] that non-Markovianity under fast driving of the WM can lead to a power boost of strictly quantum origin ("quantum advantage") which is a manifestation of the anti-Zeno effect, whereby fast driving can enhance the system-reservoir coupling [17][18][19]. An important finding has been that switching on and off the WM-reservoir coupling in a stroke has to be fully accounted for, as it may strongly influence the power output and efficiency [20,21]. Another development concerns the possibility of maintaining coherence as a resource of the machine in either the WM [22] or the piston mode [23].
On the methodological side, non-perturbative treatments of the non-equilibrium quantum dynamics have been established, such as generalized master-Floquet equations [24] and exact approaches based on the path integral representation of the reduced density operator of the WM [7]. Among the latter treatments one finds specifically the Hierarchical Equations of Motion (HEOM) [25][26][27], the Stochastic Liouville-von Neumann Equation(SLN) [28,29], the multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) approach [30][31][32], the nonequilibrium Green's functions (NEGF) [33][34][35], and the continuous-time quantum Monte Carlo (CT-QMC) algorithm [36,37] that have been recently applied to explore quantum heat phenomena. These approaches provide us with powerful tools for exploring the mostly unfamiliar domain of quantum thermodynamics far from equilibrium, particularly in the presence of low-temperature mesoscopic or microscopic reservoirs.
Current experimental studies of heat machines encompass a variety of quantum platforms: a single trapped ion as WM [38]; superconducting-circuit WM coupled to electronically engineered reservoirs [39][40][41][42][43][44], nitrogen vacancy center WM coupled to nuclear spin baths in diamond [4], phonon reservoirs in solid-state circuits [45], carbon nanotubes acting as mechanical resonators [46]; spin-wave or phonon-reservoirs in systems of trapped ions [47], as well as inelastic spin-exchange collisions giving rise to quantum-controlled heat transfer directionality [48]. Although remarkable in terms of control and fabrication, these machines mostly operate in the domain of classical thermodynamics or only exhibit limited quantum effects. Yet experiment has thus far not implemented theoretical proposals for heat machines with advantageous properties of quantum origin [ 2)]. This modulation gives rise to side-bands k ω s , k ∈ Z (orange). In this way, the TLS realizes a working medium that interacts continuously with two reservoirs at different temperatures T h > T c having localized spectral distributions at low frequencies around Ω s (slow bath) and at high frequencies around Ω f (fast bath) with respective bandwidths ∆ α . minimal designs of continuously operating quantum heat machines [16,23,[54][55][56][57]. Such continuous operation is experimentally advantageous as it does not require on-offswitching of the WM-reservoir thermal couplings and thus circumvents the conceptual difficulties that are inherent in conventional stroke (reciprocal) heat machines, for example, a four stroke quantum Otto engine [20,21,58].
Here we consider, as illustrated in Fig. 1, a minimal design of a continuously operating heat machine: the WM is a two-level system (qubit) whose transition frequency is periodically modulated [8,16,54,55,57] such that it is in alternating spectral overlap with either the hot or the cold reservoir. The innovation in the present design is that the reservoir spectra are separated by a bandgap and are therefore completely disjoint. Bandgap reservoirs arise in photonic crystals [59,60] but can also be realized in superconducting circuits by placing transmission lines or LC-oscillators as interfaces between transmon qubits and ohmic resistors.
To augment quantum effects, we adopt the strong system-reservoir coupling regime that arises near bandgaps [59]. We take the cold reservoir to be at zero temperature (T c = 0). Under such conditions, we reveal highly nontrivial features such as anharmonicity of the WM reponse, non-equilibriation of the WM in finite-time cycles as well as quantum coherences and correlations of the WM-reservoir complex. Those features are revealed by the HEOM approach [25,[61][62][63][64], a non-perturbative simulation technique of open quantum systems which is derived from the formally exact representation of the reduced density matrix in terms of path integrals [7,65]. The HEOM has the additional benefit that, as part of the hierarchy, perturbative approaches such as the Redfield and the Lindblad master equations can be easily computed by the same code. In particular, the first order truncation of HEOM gives a generalized timenonlocal master equation which allows us to find asymptotic Floquet states within a non-Markovian, perturbative method [66,67]. By comparing the results of these various approaches, insight is provided into the applicability of the perturbative treatments and the relevance of higher order quantum correlations.
This paper is arranged as follows. In Sec. 2, we present the heat engine model, whose operation essentially depends on the designed reservoir spectra. In Sec. 3, we outline the exact HEOM approach, its approximated treatment, and its use to extract heat currents. Simulation results are presented in Sec. 4, where HEOM data are shown together with the Redfield-plus approximation and a Markovian treatment. The analysis of sideband (Floquet) resonances follows in Sec. 5. The power boost of the engine that reflects the non-Markovian anti-Zeno quantum advantage is discussed in Secs. 6 and 7. The conclusions are presented in Sec. 8.

Minimal design of thermal machine with bandgap reservoirs
Quantum heat engines (QHEs) can be theoretically studied within the framework of open quantum dynamics [8,20,49,54,[67][68][69][70][71]. Within this framework, the total Hamiltonian consists of three parts, where H 0 (t) pertains to a system that acts as WM subject to external driving [41,72], H B represents two thermal reservoirs (baths) at temperatures, T h > T c , henceforth denoted as hot reservoir and cold reservoir temperatures, respectively, and H I is the coupling between the WM and the reservoirs. Here, we consider a paradigmatic WM, namely, a TLS with periodic frequency modulation, i.e., with σ ± = 1 2 (σ x ± iσ y ) written in terms of the Pauli matrices σ n , n = x, y, z. The modulation of the transition frequency is taken to be of the form whose amplitude λ and driving frequency ω s are the main variable parameters in the following analysis. We shall = k B = 1. In the context of a thermal machine, the situation, where ω(t) < ω 0 corresponds to an expansion of the WM, while ω(t) > ω 0 represents its compression. Within the set-up illustrated in Fig. 2, a conventional heat engine then operates such that the low (high) frequency bath has the lower (higher) temperature T c (T h ). However, we will also consider the reversed situation and demonstrate that it may induce even stronger heat fluxes.
The bare system Hamiltonian is diagonal in the eigenbasis of σ + σ − = (σ z +1)/2, i.e. σ + σ − |q = q|q with q = 0, 1. Consequently, the bare system time evolution operator can easily be expressed as where ω (k) = ω 0 + kω s are the quasi-energies (denoted by orange lines in Fig. 1) and J k (·) are Bessel functions of the first kind. The WM has thus multiple equidistant quasi-energy levels (in agreement with Floquet theory) which, as we will show, has prominent effects on quantum heat transport. For λ/ω s → 0 one has J |k| =0 (λ/ω s ) → 0, J 0 (λ/ω s ) → 1 so that the system effectively reduces to a static system with level spacing ω 0 . As thermal reservoirs we consider quasi-continua of independent harmonic degrees of freedom with localized spectral distributions, namely, a slow bath H s composed of oscillators in the frequency range below ω 0 and a fast bath H f with oscillator frequencies above ω 0 . Accordingly, one can write where m α,i , ω α,i , x α,i and p α,i denote the mass, frequency, coordinate and momentum of the ith oscillator in the αth reservoir. The bilinear interaction between these reservoirs and the WM is given by with X α = i c α,i x α,i denoting a collective coordinate of the αth reservoir and c α,i being the coupling constant of the ith mode with α. The last term guarantees that the reservoirs only act dynamically upon the system without any coupling-induced distortion of the system. We note that the total Hamiltonian in Eq. (1) can be transformed into the form of the conventional spin-boson model [7] by the canonical transformatioñ The effect of the two reservoirs on the WM is completely described by the coupling weighted spectral densities summed over all modes which leads to µ α = (2/π) ∞ 0 dωJ α (ω)/ω. According to our setting of two reservoirs with localized spectral distributions and only negligible overlap among them, the following continuum form of spectra is convenient, i.e., with central frequencies Ω f > ω 0 > Ω s and frequency scales ξ α which determine the widths ∆ α of the spectral distributions. In what follows, all masses are set to unity and all frequencies are scaled with ξ s = ξ f so that the dimensionless bandwidths ∆ α are of order unity. The WM continuously interacts with the reservoirs with varying spectral overlap due to the external driving (see Figs. 1 and 2). At resonance, the effective couplings read As noted above, in order to focus on the quantum regime with strong non-Markovian effects, we consider the cold reservoir to have T c = 0. This immediately implies that our machine cannot be a quantum refrigerator.

Simulation Techniques
The described setting is highly non-trivial, first, because we consider localized, nonoverlapping reservoir (bath) spectra and the machine is operated at low temperatures. Hence, correlations between WM and the thermal reservoirs are expected to be so strong that conventional master equations are not applicable. We thus rely on the exact quantum dynamical simulations within the Hierarchical Equations of Motion (HEOM) approach.

Hierarchical equations of motion
Here, we briefly describe the essence of the HEOM approach and its derivation from the path integral expression of the reduced density matrix. For the sake of simplicity, we consider only a single reservoir for a Hamiltonian of the form (5) and assume factorized initial states at the time zero and ρ m is the density operator of the relevant system (WM). The generalization to correlated initial states has also been discussed [73,74].
In path integral representation [65] the reduced density operator is obtained as Generally, a continuous system coordinate can be discretized using a system-specific discrete variable representation (DVR) [75]. Within the HEOM and for the TLS-WM considered here, a representation in terms of the eigenstates |q ∈ {|0 , |1 } is convenient. The coordinates q + (t) and q − (t) denote forward and backward system paths, respectively, and S ± [q ± (t)] the corresponding actions, These paths q ± (t) directly determine also the σ ± x (t) according to where t + denotes the time slice on the forward and backward paths that follows the time slice t. The effective impact of the reservoir onto the system dynamics is described by the influence functional [65] which reads The derivation of the real-time HEOM starts by first expanding [25,76] or fitting [77,78] the bath correlation function as a sum of exponential terms, i.e.
with the collective bath operator X as in Eq. (6) and the Bose-Einstein distribution (14) denote proper coefficients in an expansion in terms of exponentials with proper coefficients γ k . All characteristics of the reservoirs are expressed by the correlation functions which are depicted in Fig. 2 for the slow and fast reservoirs in Eq. (8), respectively. In contrast to conventional Ohmic-type spectral distributions J(ω) ∝ ω that at higher temperatures lead to an exponential decay of the auto-correlation C(t), we here observe long correlation times τ R 30 (in arbitrary units) that exceed the driving (modulation) periods τ s = 2π/ω s if ω s 2π/τ R ≈ 0.2. Consequently, a separation of time scales on which Markovian perturbative approaches are based, does not exist and memory effects are strong as we will see below. The limited bandwidth of the reservoirs yields beating patterns, particularly for the high frequency bath. By resorting to the following auxiliary density operator (ADO) definition, the HEOM formulation leads to the equation of motion [25,61,62,[79][80][81][82][83] dρ which shows that the WM interacts via σ x with the collective bath force. The ADOs ρ n s are labeled by the subscript n denoting a set of integers {n 1 , ..., n k , ...}, with n k ≥ 0 associated with the kth exponential term in Eq. (14); n + k and n − k denote {n 1 , ..., n k + 1, ...}, and {n 1 , ..., n k − 1, ...}, respectively. The super-operator acting on these ADO is defined as L 0 (t)ρ n = [H 0 (t), ρ n ]. The WM reduced density operator in this notation is ρ m = ρ {0,...0,...} .
Assuming that ρ 0 (t) is of order one, the magnitude of ρ n (t) is proportional to k d n k k , which may be divergent for strong system-bath coupling η as |n| .
According to (18) we expect the reduced density to approach a periodic steady state ρ m (t) → ρ m (t + τ s ) with τ s = 2π/ω s being the driving (modulation) period. This periodicity characterizes all single time-dependent observables (one-point correlations) such as the excited (ground) state population P e (t) = σ + σ − t (P g = 1−P e ) or the heat currents I α (t).

Perturbative treatment
Approximate treatments of open system dynamics have been developed to second order in the system-reservoir coupling. Together with a time scale separation between fast decaying reservoir correlations and relaxation dynamics of the reduced density operator, this leads to the Redfield master equation [9]. Interestingly, an extended Redfield equation is also obtained from the HEOM if it is curtailed at first-order, i.e. resricted to ADOs with k n k = 1. Namely, in the interaction picture this Redfield equation has the form where σ I x (t), ρ I s (t) and ρ I 0 + k (t) denote the system (WM), reduced density matrix and first-order ADOs in the interaction picture, respectively. Note that Eqs. (19a) and (19b) have the same structure as their counterparts in the generalized Floquet theory [66,67,70]. Upon solving Eq. (19b) and inserting it into Eq. (19a), one has In the above expression, the correlation function in Eq. (14) and the Born approximation [9] have been used but not the Markov approximation. Thus, Eq. (20) is an integro-differential equation nonlocal im time which is henceforth denoted Redfield-plus (Redfield + ) to distinguish it from the conventional Redfield formulation. It can be conventiently solved with the help of auxiliary variables [78,91,92].
In the Markov limit one sets in Eq. (20) ρ I m (τ ) → ρ I m (t), leading to the conventional time-local Redfield master equation This equation can be easily solved in the time domain with the help of Eq. (14) through the use of the auxiliary operator This substitution transforms Eq. (21) into The HEOM approach can been seen as an infinite-order extension of the Redfieldplus/Redfield approximation [80,93,94]. This allows us to reveal consistently the impact of higher order system-reservoir correlations which are particularly subtle for heat currents.

Heat current, power, and efficiency
In the framework of the HEOM, effects of the environment on the system dynamics can be obtained from the ADOs [26,63,64,95,96]. Here, we concentrate on the quantum heat current which is linear in the collective bath force X. One starts, in the interaction picture, with the following two equations where ρ I T (t) denotes the total density matrix in the interaction picture. In the next step, relations between first-order ADOs (ρ 0 + k ) in the HEOM and first-order moments of X are constructed according to Higher-order relations can be found in Refs. [64,95].
In the presence of two thermal reservoirs, the above relation is inserted into the definition [34,64] derived from first principal [5,49] to obtain the quantum heat current between WM and reservoirs, i.e.
where ν = c, h indicates cold and hot bath (either slow or fast), respectively. By definition, a positive value of I ν (t) corresponds to energy flowing into system, while a negative value corresponds to its inverse. We note in passing that an alternative definition of heat current as minus the energy change in the reservoir, which is directly related to the system energy changes when the system-reservoir coupling energy is negligible. However, the situation is very different when studying strong coupling setups where the system is periodically driven and in this case changes in the coupling energy must be accounted for [33,34,58]. Deeper insight into the operation of the thermal machine in the quantum regime is given by the normalized correlations between the slow and the fast reservoir modes mediated by the driven two level WM which can be easily obtained from the HEOM, i.e.
with k, l denoting slow and fast baths effective modes, respectively. These correlations are not present in the classical treatment nor in the standard Redfield master equation.
In steady state, we have ρ sst n (t) = ρ sst n (t + τ s ), and also I sst ν (t) = I sst ν (t + τ s ). It is thus convenient to defineĪ as an average of the heat current over one driving period. When representing the current as I st ν (t) = m I ν,m exp(−imω s t), this then impliesĪ ν = I ν,0 . The operation of the thermal machine as a heat engine is characterized by its power output and its efficiency. In steady-state, the power can be calculated directly fromĪ c andĪ h according to Here, ρ 12 0 + k and ρ 21 0 + k denote ADO elements and P e (t) = σ + σ − t denotes the excited state population. From the above definition positive values of P mean that heat is converted into work (heat engine operation), while negative P means that work is dissipated in the reservoirs (dissipator). In steady-state, a formal definition of the efficiency is given by Note that here the efficiency has a physical meaning only for positive values and is then a figure of merit for the heat engine. In addition, when T c = 0, the efficiency is limited to 1 by the first law of thermodynamics, i.e. energy conservation, which forbids an output that exceeds the input [8], otherwise, is not. With the above set of equations at hand, we are in a position to explore heat transfer properties of the quantum heat engine in more detail.

Simulation results
In this section, we present numerical results based on the HEOM and the approximate treatments Redfield + and Born-Markov Redfield, respectively. The total initial density matrix is taken to be a factorized state i.e. ρ T (0) = ρ m (0)⊗ α e −βαHα /Tr[e −βαHα ], whose time evolution is followed until a steady state is approached. In this regime, observables are calculated for various parameter sets, where thermal reservoirs are characterized by their temperature and central frequency (Ω α , T i ), α = s, f and i = h, c while bandwidths ∆ α are kept constant throughout. As a first step, we compare data obtained from the approximate treatment with exact ones from the HEOM. The task is to find parameter domains for which the Redfield + is sufficiently accurate and conversely to identify, where it fails due to very strong WM-reservoirs correlations (strong non-Markovianity).
In the numerical simulations, the spectral function S(ω) = J(ω)n β (ω) is properly fitted to an optimized rational function with tolerance δ S ≤ 10 −6 and then Fourier transformed to correlation functions C(t) written as a sum of exponential terms.  Concurrently, and on-the-fly filtering [85] algorithm is adopted in order to achieve high efficiency. Atomic units (a.u.) are used here in order to treat a variety of regimes.

Perturbative versus exact treatment
Let us first recall the relevant time scales in the periodic steady state, i.e. the external modulation period τ s = 2π/ω s and the typical correlation time (memory time) of the thermal reservoirs τ R . At low temperatures τ R ∼ β so that Markovian treatments fail and particularly τ s τ R . For bandgap environments additional time scales come into play, namely, the central band frequencies Ω α and respective widths ∆ α . Now, in Figs. 3 to 7 we compare the performance of the perturbative Redfield + with the exact HEOM for various observables. The general outcome of this comparative analysis is that we can identify parameter regimes, where the approximate treatment provides quantitatively excellent results at least for heat currents. While Redfield + accounts to some extent for memory effects in the thermal reservoirs, it does so for sufficiently weak system-bath interaction. Consequently, the dynamics of the system operator in Eq. (20) σ x (t) reflects the bare dynamics, i.e.,  However, for stronger WM-reservoir coupling (thermal contact) and/or long correlation time of the reservoir, this bare dynamics is influenced by higher order correlations between WM and reservoirs (see below). This includes higher order quanta exchange between the periodically driven TLS and the reservoir oscillator modes with frequencies around Ω s and Ω f , respectively. Figure 3 displays the dynamics of the ground state population and the heat current I h (t) in the periodic steady state with period τ s . Although Redfield + predicts a somewhat higher population than HEOM, the heat currents obtained by the two methods are in excellent agreement for the chosen parameters of relatively strong and fast driving (see Fig. 4). Figure 5a,b compares heat currents for a strictly Markovian treatment and the approximate non-Markovian approach Redfield + with exact data from HEOM. A resonance-like pattern is apparent when ω s is tuned away from the regime of very slow (adiabatic) driving to the fast driving regime. The location of these resonances is captured by the Markovian treatment which, however, completely fails to predict correct resonance heights in contrast to the accurate agreement between Redfield + and HEOM. Remarkably, broad maxima in the heat currents occurring at higher driving frequencies are completely absent in the Markovian treatment, a clear evidence that reservoir feedback and non-Markovianity play a dominant role in this domain as we will discuss in detail below. This is also true when the reservoir temperatures are interchanged such that the high frequency bath becomes the cold (hot) reservoir and the low-frequency reservoir the hot one. Then, as seen Fig. 6, the agreement between Redfield + and HEOM is less good. We conclude that within the chosen parameter domain non-Markovian effects are prominent at all driving frequencies, but particularly at faster driving and lower temperatures for the fast bath.
Even when the parameter set of Fig. 5 is adopted for the reservoirs, but the driving amplitude for the WM is substantially increased, see Fig. 7, the performance of Redfield + remains acceptable with minor deviations from HEOM only in the resonance range at slower driving frequencies again with the overall tendency to yield smaller absolute values for heat currents.
Based on the above analysis (and more systematic results that are not shown here), we conclude that for the setting studied here the Redfield + is quantitatively correct in the range of driving frequencies 0 ≤ ω s ω 0 and for driving amplitudes λ |Ω f + ∆ f /2 − ω 0 |, |Ω s − ∆ f /2 − ω 0 | (so that ω(t) does not exceed the full bandwidths of the reservoirs during one period) as long as Ω f /T is on the order of 1 with T being the temperature of the high frequency reservoir (i.e. the reservoir memory time is not too strong).

Sideband resonances
In order to understand the behavior of the heat currents when the driving frequency ω s is varied, we start with the range of slow to moderate driving, where, as we will show, weak coupling master equation predict at least qualitatively resonant behavior. This is no longer true for faster driving as we will discuss in the next Section.
In the regime of weak coupling he interaction between the WM and the reservoirs is governed by the real-valued transition rates which directly determine the heat fluxes (cf. Appendix). Consequently, the reservoirs are only probed at the central (ω (0) = ω 0 ) and sideband (ω (k) , |k| = 0) resonances, where only the sidebands contribute to the heat current and the contribution of each sideband is weighted by the Bessel function J k (λ/ω s ), see Eq. (4). The underlying approximation requires that ω s ω 0 , Ω α . This treatment leads to the prediction that for the distributions of Eq. (8) around Ω α , we expect a resonance-like pattern for the heat currents around ω (k) ≈ Ω α and thus around the driving frequencies In Fig. 5 this Markovian prediction for the heat current is shown together with the prediction from Redfield + and HEOM. For the parameters chosen (λ = |Ω α −ω 0 | = 1), all treatments yield pronounced resonances around ω (1) s = 1 (single-quantum exchange) and ω (1) s = 1/2 (two-quanta exchange). However, the accuracy of the Markovian description is rather poor: the precise location of the resonance is shifted from the Markovian resonance condition and the peak heights differ substantially from the exact ones. The Markovian approximation completely fails in the limit ω s → 0. In a more accurate description, higher order quanta-exchange resonances are blurred by the steep decrease of the heat currents towards ω s → 0. Apparently, the resonances are broadened by the finite bandwidths of the reservoirs of order ∆ α /2k, see Figs. 5-7. Memory effects of the reservoir response become even more prominent when the temperature gradient is reversed as in Fig. 6, where the high frequency reservoir has T c = 0. This induces much stronger memory effects (non-Markovian behavior) on time scales of order τ s so that (i) resonances occur slightly away from ω (k) s and (ii) higher order system-reservoir correlations cannot be neglected. They broaden the resonances and increase their magnitudes, thus, demonstrating that 'deep' quantum effects enhance the heat transfer. With increasing driving amplitude, see Fig. 7, ω(t) covers the full bandwidths of the reservoirs, i.e. max t {ω(t)} ≈ Ω f + ∆ f /2, min t {ω(t)} ≈ Ω s − ∆ s /2, so that resonances overlap.
With respect to the heat power, we observe in Figs. 5 and 6 the expected behavior of peaks at the sideband resonances according to the respective temperature gradients. This is shown in more detail in Fig. 8, where the heat power turns from being positive to being negative when the temperature of the low frequency reservoir increases. Note that in all cases, finite heat power also appears outside the range, where sideband resonances exist (ω s > ω (1) s ), i.e. outside the Markovian domain. We will discuss this latter range in the following.

Power boost by non-Markovianity for fast driving
We also observe broad extrema in the heat currents for modulation frequencies around ω s = 2. To explain their nature, we consider the limit, where the two reservoirs collapse to single oscillator modes with frequencies Ω α . Accordingly, one has with the corresponding Heisenberg equations of motion One can iterate the above equations to obtain for the reservoir oscillator modes with the indexᾱ = f, s for α = s, f . To order c 2 α , this equation describes linearly and parametrically driven harmonic systems while higher order couplings induce nonlinearities.
The time dependent heat current follows according to (26) from I α (t) ∝ c α σ y q α t . Since to leading order q α α = 0, the c α -dependent terms in Eq. (34 c) for σ y (t) are relevant which implies I α ∝ c 2 α q α (t)q α (s)σ z (s) . To second order in the couplings, σ z (t) in I α carries the bare frequencies ω (k) and one regains the resonance condition Eq. (32) for the time averaged heat current. Beyond this approximation, the coupled dynamics in Eqs. (34) and (35)  being integer) and their combinations, thus giving rise to beating between WM and the reservoirs as well as between the reservoirs. Hence, in the time averaged heat current we expect in the range ω s > ω (1) s leading order resonances behavior at the frequencies ω s ≈ Ω α , 2Ω α , at ω s = (ω 0 + Ω α )/2 and ω s ≈ (Ω f ∓ Ω s )/2.
The above picture can be conveniently verified by considering narrow spectral distributions of the form which describe weakly damped oscillator modes with effective damping rates (spectral widths, cf. Fig. 1) ∆ α Ω α , see Fig. 9(a). Corresponding results for the heat current are shown in Fig. 9 Fig. 10 where extrema are in complete agreement with the above predictions.
Of particular interest are reservoir-reservoir correlations that we expect to emerge in the domain of strong non-equilibrium, i.e. for ω s > ω (1) s . As Fig. 11 reveals, in the latter range these correlations are continuously built up with increasing modulation frequency approaching a constant level. In contrast, under slow and moderate modulation, these correlations only appear near resonances. Note that X s X f −correlations do not exist in the Born-Markov approximation (which is a perturbative treatment up to order c 2 α in the secular approximation). These correlations reflect higher order systembath contributions, at least of order c 4 α , where they match the size of higher order contributions of auto-correlations (cf. the normalization in Eq. (27)).
Hence, from the above analysis we conclude that in addition to sideband resonances predicted from a Markovian treatment for slow to moderate driving, non-Markovian feedback effects dominate for faster modulation and lead to finite heat currents which would be absent otherwise. This effect consitutes non-Markovian power boost. This is further demonstrated in Figs. 12 and 13, which depict the dependence of power and efficiency on the modulation amplitude and the coupling to the reservoir. While the maximum power around ω s ≈ 2 is sensitive to the driving strength, the efficiency is much less affected and even decreases for stronger modulation. A similar picture emerges for growing coupling: Stronger coupling provides more heat power but does not enhance the efficiency around the power maximum. Beyond the maximum of the power the engine turns into a heat dissipator with collapsing efficiency. Most remarkably, even without spectral overlap of the response frequency sidebands with the reservoir spectral density, at the modulation constant λ = 0.5, for which the modulated ω(t) always remains within the bandgap between the reservoirs, we find substantial heat power with high efficiency. The reason is that response broadening due to the quantum time-energy uncertainty relation gives rise to the required spectral overlap in the non-Markovian anti-Zeno regime [16].

Power boost through modulation of the thermal contact
So far we have considered a continuous coupling (denoted by scheme I) of the WM to its reservoirs. In order to further explore the machine performance in the non-Markovian domain, we extend the analysis to protocols which modulate periodically the WMreservoir thermal contact (coupling strength). For this purpose, we choose two different protocols. Both start with an initialization step during which the engine approaches a steady state density from a factorized initial state. After this step thermal coupling is periodically switched on and off according to schemes: (II) Abrupt decoupling/coupling, as described in [16], can be simulated by setting all ADOs in HEOM to either zero or non-zero. (III) Spectral decoupling/coupling by intermittently moving ω 0 to a very high frequency above both reservoir bands and back to the reservoir gap, while still maintaning the modulation ω(t) in Eq. cycle is described by a time dependencē with θ(t) denoting the Heaviside step function and ω 1 ω 0 , Ω α , λ being the offresonance transition frequency. The time t 0 is chosen such that t m /t 0 1 so that, effectively, the WM is able to exchange heat with the reservoirs only during time spans of duration 2t m while for times spans of duration (n + 2)t m it is decoupled. The cycle time is thus τ C = (2n + 6)t m with an integer n such that (n + 2)t m is on the order of the memory time of the reservoirs τ R .
The resulting time dependent heat currents are depicted in Fig. 14, lower panel. Note that the strongly oscillating part during and after the switch from ω 0 → ω 1 averages basically to zero so that the net heat is exchanged only when the transition frequency is modulated around ω 0 . The power performance is shown in Fig. 15 for both schemes II and III compared to a continuous process (scheme I, no active decoupling). We recall that a Markovian treatment predicts zero power in this case.
Both segmental schemes II and III lead to an enhanced power in the regime around ω s = 2, where heat currents are now averaged over a full cycle. The abrupt decoupling scheme II turns the engine into a dissipator and then again into a heat engine for higher frequencies, while in the continuous decoupling/coupling scheme III it operates as heat engine throughout the range of modulation frequencies. While scheme II is more of an idealized thought experiment, scheme III can be implemented experimentally. These power boosts constitute a quantum advantage in agreement with the anti-Zeno effect [16] that stems from the time-energy uncertainty relation under fast modulation of the system-reservoir interaction. Here the anti-Zeno effect arises in the strong-coupling regime and is evaluated non-perturbatively.

Conclusion
This paper has studied a minimal model for a quantum heat engine in a non-trivial setting: By periodically modulating the transition frequency of a TLS acting as WM, this WM is in alternating spectral overlap with a cold or a hot reservoir with localized spectra that are separated (disjoint) by a bandgap. Such a continuously operating heat engine [8,54,55,57] has advantages over conventional stroke engines (for example the four-stroke Otto engine) since it does not require abrupt on-off switching of the coupling between WM and the reservoirs. As a result, the energetic and entropic cost of the switching is avoided. The present design extends the previous proposal [16] to a much broader range of operational conditions, from weak to strong coupling of the WM and the reservoirs and weak to strong or slow to fast driving/ modulation of their coupling strength. Predictions have been obtained under any of the foregoing regimes without restrictions, based on the formally exact HEOM simulation technique and its comparison with perturbative treatments. While Markovian treatments predict vanishing heat flow at faster driving, the perturbative Redfield + can provide quantitatively correct predictions of the quantum heat engine performance for all driving frequencies. For slow and moderate driving the WM exchanges energy with thermal reservoirs via multi-sideband resonances.
The essential feature of the reservoir spectral functions considered here is the presence of a bandgap with spectrally abrupt edges, because it gives rise to strongcoupling effects [59], and allows for disjoint hot and cold reservoir spectra, as required for efficient heat engines that operate continuously [8,16,23,54]. As noted in the Introduction, such bandgap reservoirs can be realized in photonic crystals [59,60] their phononic analogs [44][45][46][47][48][49][50][51][52][53]55], and are envisaged also in superconducting circuits. Although the exact spectral shape of the reservoirs outside bandgaps is of no qualitative importance, we note that reservoir spectra and the resulting energy (Lamb) shifts can be engineered using the principles discussed in refs. [42,43,55].
The highlight finding has been the quantum anti-Zeno advantage of the thermal machine for both continuous and segmental modulation in the deep non-Markovian regime. In the former case, strong bath feedback emerges due to non-equilibrium processes, while in the latter case, strong memory effects govern the quantum dynamics. The latter modulation protocol can be implemented in actual experimental settings.

Author contributions
M. X. performed numerical simulations. All authors have been involved in model setting, results analysis, discussion of scientific results and in the writing of the manuscript.

Data availability
The data that support the figures within this article are available from the corresponding author upon reasonable request.

Appendix: Heat flux in Born-Markov approximation
The population dynamics of a driven two level system interacting with a bandgap reservoir as considered in the main text is governed by [16] and P 1 (t) = 1 − P 0 (t). Asymptotically forṖ i (t) = 0, this equation can be solved as with transition rates Γ 0/1 = k Γ Accordingly, the heat currents in steady state are given by with ω (k) as in Eq. (4) and the population ratio