From quantum heat engines to laser cooling: Floquet theory beyond the Born-Markov approximation

We combine the formalisms of Floquet theory and full counting statistics with a Markovian embedding strategy to access the dynamics and thermodynamics of a periodically driven thermal machine beyond the conventional Born-Markov approximation. The working medium is a two-level system and we drive the tunneling as well as the coupling to one bath with the same period. We identify four different operating regimes of our machine which include a heat engine and a refrigerator. As the coupling strength with one bath is increased, the refrigerator regime disappears, the heat engine regime narrows and their efficiency and coefficient of performance decrease. Furthermore, our model can reproduce the setup of laser cooling of trapped ions in a specific parameter limit.


Introduction
In the recent years, notable progress has been made towards the experimental realization of small-scale thermal machines [1][2][3][4][5]. Specific examples in the quantum regime include a quantum absorption refrigerator with trapped ions [6] and quantum heat engines using an ensemble of nitrogen-vacancy centres [7]. The theoretical study of these quantum thermal machines has in general been restricted to the weak coupling and Markovian regime [8][9][10]. This approach follows the logic of traditional thermodynamics, which is a weak coupling theory and relies on the distinction of systems interacting with each other due to the negligible contribution of their interface compared to the bulk properties [11]. Nevertheless, this assumption becomes increasingly questionable at the nanoscale where, as the volume of systems becomes very small, boundaries cannot be clearly distinguished. In the present work we will present a formalism which allows to formulate a consistent thermodynamic framework for periodically driven systems arXiv:1712.07032v2 [quant-ph] 15 Jun 2018 coupled to multiple heat reservoirs -even when the coupling is strong, driven and induces non-Markovian behaviour.
Time-periodic systems have been successfully studied by means of Floquet theory [12][13][14] and a broad diversity of interesting phenomena such as coherent destruction of tunnelling [15,16], dynamical localisation [17] and creation of new phases of matter [18][19][20] have been discovered. When connected to a heat reservoir, the steady-state dynamics of time-periodic quantum systems has been mainly described using the Floquet-Markov approach [14,[21][22][23][24][25][26]. It consists in deriving a weak coupling master equation [27] in the Floquet basis of the driven system. Under the secular approximation, conditions for the emergence of a Floquet-Gibbs state have been studied recently in [28], a full stochastic thermodynamic analysis is given in [29] and a thermodynamic analysis of laser cooling experiments by collisional redistribution [30,31] can be found in [26].
To access the dynamics or study thermal transport in periodically driven systems beyond the Born-Markov approximation, more sophisticated methods must be used such as stochastic Liouville-von Neumann equations [32], perturbative high-frequency expansions [33] or influence functional integral methods [34]. Their thermodynamic interpretation, however, is not always clear. Recent studies beyond the Born-Markov approximation with a consistent thermodynamic interpretation include non-equilibrium Green's functions [35,36], redefined system-reservoir partitions (collective coordinate mappings) [37][38][39][40][41], a quantum absorption refrigerator [42], quenched thermodynamic protocols [43,44] and the derivation of an exact expression for the entropy production of a finite arbitrary system in contact with one or several thermal reservoirs [45], but none of which were applied to periodically driven machines so far. Exceptions are [46] and [47]. In [46] a polaron transformation and Floquet theory were used to include arbitrary strong coupling effects but only as long as the rotating wave approximation for the system Hamiltonian and the Markovian approximation for the reservoir hold. In [47] a periodically driven three-level heat engine was studied using the numerical method of hierarchy of equations of motion, an approach that is based in a decomposition of the bath correlation function rather than an explicit representation of the bath. For this reason, access to more general, thermodynamically relevant thermal transport properties requires a tailored solution [48].
However, none of the aforementioned methods was successfully applied to the combinations of problems we want to tackle: a driven system coupled to multiple heat reservoirs with a possibly strong, driven and non-Markovian coupling. To achieve this we combine three methods: a collective coordinate (CC) mapping [49,50], Floquet theory for open systems [14,21] and full counting statistics [51]. This novel and unique combination provides access to steady state thermodynamics lifting the restriction of common assumptions such as a very fast driving, Markovian and weakly coupled reservoirs, and the secular approximation. First, in the CC mapping, we identify a collective degree of freedom in the reservoir, sometimes referred to as reaction coordinate [52] that is responsible for the strong coupling and non-Markovian effects. Second, using Floquet theory, we derive a master equation for the original system plus CC and apply full counting statistics methods to obtain the change in energy of the reservoirs unambiguously. This strategy allows us to perform a consistent thermodynamic analysis of periodically driven thermal machines. Furthermore, it also allows us to accurately treat periodic time dependencies in the interaction between system and reservoirs, capturing their thermodynamic effects, which are inaccessible with standard methods. It is worth mentioning that even though the master equation that will be used to study the original system plus CC is Markovian, once the degrees of freedom of the CC are traced out, the dynamics of the system is non-Markovian.
This work is organized as follows: We start by briefly presenting the CC mapping (Sec. 2.1), explaining the Floquet master equation used including the technique of counting fields (Sec. 2.2) and defining important thermodynamic quantities (Sec. 2.3). In Sec. 3 we begin by specifying the model for study. Results are then presented for the weak coupling and non-Markovian regime (Sec. 3.1) with the identification of the refrigerator and heat engine regimes. In Sec. 3.2 we explore the regime of finite coupling and establish the existence of two additional modes of operation, whereas in Sec. 3.3 we look at the performance of our thermal machine. Finally, we study the limit in which we reproduce the setup of laser cooling (Sec. 3.4) and follow with conclusions (Sec. 4).

Method
Our formalism is in general applicable to situations modeled by a system-bath Hamiltonian of the form where H S (t) is the Hamiltonian of the system, H with frequency ω kν , mass-weighted positions x kν and momenta p kν of the bath oscillators and S ν (t) a system operator which mediates the coupling to the bath with strength c kν . We consider the case of a periodic time dependence in the working medium H S (t) = H S (t + T ) and coupling operator S ν (t) = S ν (t + T ), with period T = 2π/ω L . The periodicity will allow us to study the problem by using Floquet theory.
To fully characterize the model, the spectral density of the heat reservoirs, defined by must be parametrized. The spectral density is a positive function for ω > 0 and fulfills J ν (ω) → 0 for ω → 0 and ω → ∞. A structure-less spectral density (e.g. linear form) Figure 1. Sketch of the model before and after the CC mapping. A driven system S coupled to a hot bath B (h) and a cold bath B (c) is mapped to a system coupled to the same hot bath and the CC, which is coupled to a residual bath B (c) . The oscillating arrows indicate the influence of the driving. The CC and the residual reservoir B (c) are related by a unitary transformation to the original reservoir B (c) .
usually allows for a Markovian treatment of the reservoirs due to the fast decay of its associated correlation functions, while a more structured spectral density (e.g. strongly peaked around a frequency ω res ) demands a more elaborate treatment. In the following, we will consider the case of a hot (ν = h) and cold (ν = c) reservoir where only the cold spectral density has a structured form: The parameters d c and d h describe the overall coupling strength to each respective reservoir. The spectral density of the cold reservoir is peaked around resonance frequency ω res and its structure can be tuned by the parameter γ, where the smaller γ the more strongly peaked the spectral density is. The parameter ω 0 is a reference frequency of the working medium. We will also consider only the coupling operator of the cold bath S c (t) to be time dependent. Figure 1 (top panel) shows a sketch of the model.

Collective Coordinate Mapping
To be able to capture non-Markovian and strong coupling effects of our model we follow the procedure discussed in [37-41, 53, 54] and introduce a CC, also known as reaction coordinate, from the reservoir as part of the system. This is done by a unitary transformation applied to the bath degrees of freedom, which we will apply to the cold bath only because of its structured form. In the following we will refer to this method as the collective coordinate mapping.
The CC is defined as with λ 0 an unspecified parameter so far. After the mapping (see figure 1), we obtain with a redefined "supersystem" Hamiltonian and residual cold reservoir B described by Here, X 1 and P 1 are the position and momentum operators of the CC and X k and P k position and momentum operators of the residual cold bath B . The CC has natural frequency ω CC and the spectral density of B is defined as J c (ω) ≡ π 2 k>1 C 2 k Ω k δ(ω − Ω k ). In the continuum limit, all new parameters and J c (ω) can be expressed in terms of the original spectral density J c (ω) (see Appendix A or [49]). For our choice of spectral density, equation (4), we have J c (ω) = γω, ω CC = ω res , λ 0 = d c .
(9) Figure 1 shows a sketch of the original and mapped model with the newly defined supersystem S . In terms of creation and annihilation operators the new working medium supersystem Hamiltonian can be written as One of the key features of the mapping is that an increase in the interaction strength between working medium and cold bath only increases the interaction between working medium and CC. The coupling strength between CC and residual bath is unaffected (see equation 9 or Appendix A), allowing to treat the supersystem with standard weak coupling master equations when γ can be considered small. The CC mapping has been used to study an Otto cycle stroke type engine [38], stochastic thermodynamics based on coarse-graining [41], continuously coupled but not periodically driven engines [37], a fermionic autonomous Maxwell demon [39] and a fermionic electronic Maxwell demon [40], all in the strong coupling and non-Markovian regime. It has offered an accurate method for the study of open quantum systems apart from the thermodynamic applications [53][54][55] and is closely related to the method of time evolving density matrix using orthogonal polynomials algorithm (TEDOPA) [50,[56][57][58].
It was, however, not applied to the study of periodically driven systems so far. We point out that even though the CC mapping was applied here to the specific case of a spectral density with a single peak (see equation (4)), spectral densities with many peaks or general non-Ohmic forms may be addressed with the mapping if the use of multiple collective coordinates is considered [59].

Floquet Master Equation with Counting Field
We want to obtain a master equation for the time-periodic supersystem S that permits us to study its thermodynamic properties based on a rigorous framework. For this purpose we introduce a counting field χ ν for each reservoir as in [29,51] (see Appendix B) and apply Floquet theory for open systems (see Appendix C and Appendix D).
Formally, the counting field χ ν is introduced by defining the modified density matrix with total (supersystem plus reservoirs) initial density matrix ρ tot (0) and modified evolution operator is the evolution operator corresponding to the Hamiltonian of equation (6). The introduction of a counting field will allow us to obtain the change in energy of each reservoir from the operator ρ(χ ν , t) ≡ Tr B {ρ tot (χ ν , t)} (see Appendix B for details). To compute the reduced density operator dressed with counting field ρ(χ ν , t), a generalized master equation can be derived by performing the Born and Markov approximations to equation (11). In the Schrödinger picture it has the form ∂ t ρ(χ ν , t) = L(χ ν , t)ρ(χ ν , t), where L(χ ν , t) is a time periodic superoperator that has the same frequency as the driving. The master equation for the reduced density operator of the supersystem is then obtained by taking χ ν = 0. In general, ρ(t) is not a time periodic function, but we are only interested in the dynamics in the long-time limit (t → ∞). Here, we expect ρ(t) to acquire the same periodicity as the driving, such that the master equation in the long-time limit may be written as where ρ q and L q are Fourier components defined by ρ(t) = q e iω L qt ρ q and L(t) = q e iω L qt L q , respectively. The change in energy of the reservoir can then be obtained from (see Appendix B) where the prime indicates a derivative with respect to χ ν . Although this formulation is associated with the statistics of two consecutive measurements [51], the effect of the first measurement only shows up in high order moments and can be accounted for by means of a generalized fluctuation theorem [48].
In the standard weak coupling approach, it is usually assumed that the system's level broadening is much smaller than its level spacing such that non-resonant terms may be neglected. This is known as the secular approximation [27], whose advantage is that the master equation generator is of Lindblad form. It also has the advantage of simplifying the study numerically since at steady state one only has to deal with a vector of populations and not the full density matrix. However, the secular approximation is problematic here because, due to the periodic time dependence of H S (t), we resort to Floquet theory for the solution, where the dynamics of the supersystem is determined by its quasienergies (chosen to lie in the same Brillouin zone) and respective Floquet modes (see Appendix C). Due to the high (infinite) dimension of the CC, there is a large amount of near-degenerancies in the chosen Brillouin zone, thus making the secular approximation in general not suitable for our analysis [60].

Steady state Thermodynamics
Due to the periodicity of the working medium, expectation values in the long-time limit will asymptotically oscillate as a function of time with period T . It is natural then to consider their average over one period, i.e. A = T 0 dt A/T . The first law (energy balance) and second law (positivity of entropy production rate) in the long-time limit stipulatė −β hQh − β cQc ≥ 0 .
Due to the weak coupling between supersystem and reservoirs, we identify the change in energy of the reservoir, calculated using the counting field χ ν , with the heat floẇ Q ν [37], defined to be positive if it flows into the system. For the hot bath and cold bath we define respectivelẏ The rate of workẆ performed on or exerted by the system is then fixed by the first law (14). Negative values indicate that work is extracted. Positivity of the entropy production rate, equation (15), can be shown from the definition of entropy production [45] and the existence of a periodic steady state in the long-time limit. Under the assumption of weak coupling between supersystem and reservoirs, our use of master equation (12) showed no violation of the second law at steady state for all presented calculations.
The relation between the change in energy of the residual cold reservoir B and the change in energy of the original cold reservoir B is given by the change in energy of the CC. To see this, note that in the CC mapping the reservoir Hamiltonian is mapped to three different terms (see equations (6), (7) and Appendix A) CCB . As the coupling between CC and residual reservoir is weak, we can approximate the expectation value of the residual reservoir by In the long-time limit, if we average over one period, we have B . Note that the CC mapping allows for a thermodynamically consistent study of a time-dependent interaction S ν (t) between the system and the cold reservoir. This is due to the fact that, after the mapping, the time-dependent coupling operator only appears in the supersystem Hamiltonian. The interaction with the residual reservoir has no time dependence (see Sec. 2.1 or Appendix A).

Results
All previously presented ideas apply to the case of an arbitrary working medium described with Hamiltonian H S (t) and arbitrary coupling operators S ν (t). We now focus specifically on the case where the working medium is a driven two level system Here, σ i denotes a Pauli matrix, ω 0 > 0 is the energy splitting of the two-level system in the absence of driving, g the driving strength and ω L is the frequency of the driving responsible for work extraction and injection. The system coupling operators are where we have included a time dependence on the cold reservoir coupling operator. The mapped working medium consist then of a driven two level system coupled with the CC via a time dependent interaction. This bears resemblance to the setup of laser cooling of trapped ions [61], where an ion interacts with an electric field. The ion is approximated by a two-level system and its motional degree of freedom is quantized in a harmonic potential. Assuming the ion is well confined, one ends up with an interaction Hamiltonian between the ion and the harmonic potential that is time-dependent in a similar fashion as the coupling operator S c (t) (see Appendix F). Cooling in this context refers to the preparation of the state of the harmonic mode in a low occupation number. It is optimized around the resonance condition where the frequency of the laser and of the harmonic potential add up to the energy splitting of the two-level system. Based on the similarities of this setup and our mapped model we explore this resonance condition which for our model reads Importantly, the setup of laser cooling follows from an exact limit of our model by considering a weak amplitude driving g ω 0 and taking the limit γ → 0. In this case, there is no coupling to the residual reservoir, we have a two level system (ion) coupled to the CC (vibrational mode of the ion) and one bath (see inset in figure 5). This limit will be considered later in Sec. 3.4. For now we keep γ finite so that the supersystem is coupled to both reservoirs.

Weak coupling
First, we will consider weak coupling between working medium and cold bath. In the mapped model this translates to a weak interaction λ 0 = d c between the two-level system and CC. It is important to note that even when the coupling is weak, the structured shape of the spectral density can have a strong effect on the behavior of the system invalidating a conventional master equation approach for the two-level system (also see Appendix E for benchmarks related to that question). Figure 2 (bottom panel) shows the heat flow to the two reservoirs and the work flow as a function of the driving frequency while the resonance condition is maintained. Two main regions can be identified in figure 2 by looking at the conditionsẆ < 0 anḋ Q c > 0. When one of them is fulfilled the model either behaves as a heat engine or as a refrigerator. The first region is illustrated by the red (left) shaded area. Here, heat flows from the hot reservoir into the systemQ h > 0, work is extractedẆ < 0 and heat flows from the system into the cold reservoirQ c < 0 making our model act as a heat engine. The blue shaded area shows the region where the model acts as a refrigerator. Here, … … ω CC Figure 3. Mapped working medium acting as a refrigerator (dashed arrows) or as a heat engine (continuous arrows) at the resonance condition (20). |e and |g refer to the excited and ground state of the two-level system. |n refers to the Fock state n of the CC.
work is applied to the systemẆ > 0, heat flows from the system to the hot reservoiṙ Q h < 0 and from the cold reservoir to the systemQ c > 0. We also note the existence of a small finite white region separating the transition between refrigerator and heat engine not seen in a previous study [22], based on the standard Floquet-Markov approach, for a similar but simpler model.
To physically explain the heat engine and refrigerator regime, one can use figure 3. It shows the energy levels of the mapped working medium [see equation (10)] for very weak coupling and driving amplitude. In this parameter regime the cold reservoir mainly induces transitions between states of the CC, leaving the state of the two-level system unchanged. This is illustrated for example by the transitions |g, n ↔ |g, n − 1 (blue arrows). On the other hand, the hot bath preferentially induces transitions that leave the Fock state of the CC fixed and flip the state of the two-level system as in |e, n − 1 ↔ |g, n − 1 (red arrows). Under the chosen resonance condition, the driving part addresses transitions that alter both the state of the two-level system and the CC, and they are responsible for work extraction or injection. Such transitions are illustrated by e.g. |e, n−1 ↔ |g, n (green arrows) . Fig 3 only shows the most dominant transitions, additional transitions are induced by the reservoirs or the driving but they are suppressed due to the weak coupling or the resonance conditions and become less likely.
To further understand why our model behaves either as a refrigerator or a heat engine, let us consider the effective inverse temperature of the CC β CC , implicitly defined by the equation where a † a is the occupation number of the CC in the long-time limit averaged over a period. Whenever β CC is bigger than β c heat flows from the cold bath into supersystem S and increases the Fock state of the CC, therefore cooling the cold reservoir. This cooling is achieved by the injection of work simultaneously exciting the two-level system and lowering the state of the CC followed by heat flowing from the supersystem into the hot reservoir via a decay of the two-level system. This process is illustrated by the dashed arrows in figure 3. In the case that our model acts as a heat engine, heat flows from the hot reservoir to the two-level system and work is extracted by simultaneously exciting the CC and the two-level system decaying. Also, since now β CC < β c , heat no longer flows from the cold reservoir into the system. This process is indicated by the continuous arrows in figure 3. The top panel of figure 2 shows a plot of β CC as a function of the frequency. We see that the effective temperature of the CC (dashed line) and the temperature of the cold bath (continuous line) agree precisely at the boundary of the refrigerator region. We stress the fact that the discussion so far only applies for weak driving strength and weak coupling between working medium and cold bath, since otherwise the energy level diagram of figure 3 does not apply. The implicit definition of the effective temperature in equation (21) is just a particular parametrisation. We have analysed the fluctuations around the average occupation number of the collective coordinate and found good agreement with those typical for a thermal state (see Appendix G). In addition, we have checked explicitly that in this limit, the populations of the CC are close to a thermal distribution.

Finite coupling strength
As mentioned earlier, the transition between refrigerator and heat engine is not immediate, there exists a gap between the two regimes. Although already present in figure 2, we further explore this gap in the phase diagram of figure 4 (left) where four different regions can be identified depending on the coupling strength between working medium and cold reservoir. Regions I and IV indicate the regions previously introduced, where the model behaves as a heat engine or as a refrigerator, respectively. In both of the two middle regions (II and III) work is being done on the system and heat flows from the system into the cold reservoir. Regions II and III are distinguished by the direction of the heat flowing between the hot reservoir and the supersystem. For region III, heat flows from the system into the hot reservoir and for region II the opposite way. In region II work is applied to the system to assist the flow of heat from the hot reservoir to the cold reservoir and in region III work is converted into heat flowing into both reservoirs. The existence of these regions should not be ignored since their area grows as the coupling is increased. The refrigerator regime actually completely disappears. The immediate transition from a refrigerator to a heat engine as the one obtained in [22] only occurs in the ideal case of vanishing coupling strength.

Performance
To quantify the performance of our model we introduce the efficiency of the heat engine η and the coefficient of performance (COP) of the refrigerator κ. Both have ideal Carnot bounds given in terms of the bath temperatures. They are defined as: Figure 4 (right) shows the efficiency and COP of our model for different coupling strengths. The shaded area only applies for the weakest coupling case d c = 0.001ω 2 0 . A clear decrease in performance can be seen as the coupling increases for both, heat engine and refrigerator. This is consistent with performance results obtained for the quantum Otto cycle [38] and a continuously coupled but not driven three-level heat engine [37]. We also see that, as the coupling decreases, the maximum in efficiency η for the heat engine approaches the border between the white and red region, where the power goes to zero (see figure 2).
The exploration of stronger coupling regimes than the ones in figure 4 (left) requires a high level of truncation for the CC for convergence, making numerical simulations very demanding. This is also the case for a blue detuned resonance condition where ω L > ω 0 , since the CC is heated increasing its occupation number, and no proper truncation of the Fock levels may be defined. In all presented results the dimension of the CC was considered finite and truncated at a level where convergent results were obtained.

Single Reservoir: Laser cooling limit
We now consider the case of vanishing coupling between supersystem and residual bath B : γ → 0. This means that the bath is exclusively represented by the CC, and therefore there is no need to treat it perturbatively. The supersystem is then coupled only to a The goal of such experiments is the preparation of the harmonic mode in its ground state, and the success is usually measured in terms of a low occupation number of the vibrational mode n . The inset of figure 5 shows a sketch of this set up. Figure 5 shows the occupation number of the CC as a function of the detuning ∆ = ω L − ω 0 for different driving amplitudes. We see that the optimal preparation (maximum cooling) is obtained for small driving amplitudes and around the resonance condition studied previously ∆ = −ω CC . The continuous line is obtained from the analytical results in [61] where they showed that at steady state the occupation number is given by n = A + A − −A + , with A ± = Γf (∆ ± ν) and f (∆) = (Γ/2) 2 (Γ/2) 2 +∆ 2 , ν the frequency of the oscillator and Γ the decay rate of the two-level system (ion). In figure 5 we take ν = ω CC and Γ = d h . It is important to state that the analytical results in [61] are obtained under approximations that we do not perform here such as the rotating wave approximation and the adiabatic elimination of the internal degrees of freedom of the ion. We have also only considered the case of low amplitude driving such that scattering effects, as the recoil of the ion in the emission of a photon, can be neglected.

Conclusions
We have provided a novel framework to study periodically driven thermal machines beyond the weak coupling and Markovian approximations. The CC mapping was used to overcome the loss of separability that appears for systems as their interaction strength with the bath increases, Floquet theory allowed us to accurately treat periodic time dependencies of working media and, with the use of full counting statistics, we were able to consistently perform a steady state thermodynamic analysis.
We specifically considered as the working medium a driven two-level system with time-dependent interaction to one of its reservoirs and dealt with the case where only this heat reservoir had a structured spectral density. We focus on the limit of a strongly peaked spectral density since that also guaranteed a weak coupling to the residual bath such that a master equation could be used. In the weak coupling and non-Markovian scenario, we confirmed the existence of two main regimes where the thermal machine works either as a refrigerator or a heat engine and gave an intuitive physical explanation of the observed difference and transition between these two regimes in terms of the occupation number of the CC and its effective temperature. However, we also additionally identified a small gap between both operational regimes.
Upon increasing the interaction strength between working medium and the cold reservoir, we observed an enlargement of the gap between the refrigerator and heat engine regime. We identified four different operation regimes for our model and witnessed the eventual disappearance of the refrigerator regime as the coupling increased. In terms of efficiency, both the refrigerator and the heat engine showed a monotonic decrease of performance as a function of the coupling strength. Similar performance results for non periodically driven systems were observed in [37,38] and higher efficiency in the weak coupling but non-Markovian regime was observed in [37]. Nonetheless, as it was the case here, these results were focused on particular models, and the general question of whether strong coupling and non-Markovian effects are (un)favorable for the performance of quantum thermal machines is still open. Strong coupling corrections to the second law due to the interaction strength were derived in [44], also showing a lower performance, but restricted to thermodynamic protocols where the system is adiabatically driven and never simultaneously coupled to more than one thermal bath.
We also considered the case where the inclusion of the CC is enough to capture all the effects of one of the two baths, leaving the supersystem coupled only to one bath. By doing this we recovered the setup of laser cooling of trapped ions. The predicted occupation of the harmonic mode (CC) was shown to be in good agreement with the analytical results of [61], confirming the adequacy of the proposed method.
Finally, we emphasize that the method presented here applies in general to any periodically driven thermal machine coupled linearly with a bosonic bath, so that more complex models may be addressed in the future.

Appendix A. Collective Coordinate Mapping
In the following it will be assumed that the system (supersystem) is in contact only with one bath, but additional reservoirs can be incorporated additively. The relevant quantities of the mapped Hamiltonian are given in terms of the original spectral density [49] as with J 0 (ω) = J c (ω) and frequency of CC ω CC = λ 2 k Ω k δ(ω − Ω k ) which can be related to the original spectral density by W 0 (ω + iε) and P indicates the principal value. Note that J 0 (ω) is extended to negative values of ω via J 0 (−ω) = −J 0 (ω). From these expressions it is easy to see that increasing the coupling between system and reservoir by J 0 (ω) → αJ 0 (ω) (with α = 1) has no effect on the coupling between supersystem and residual reservoir since the multiplicative factor α cancels out in J 1 (ω). It is important to note that the CC mapping can actually be applied more than once. If one wishes to apply the CC mapping recursively, it is then necessary that the behaviour of the spectral density J n (ω) as ω → ∞, guarantees convergence of expressions for the residual reservoir with spectral density J n+1 (ω). The index n = 0 indicates the original spectral density of the bath. This can be guaranteed from the start by introducing a hard cutoff in the original spectral density J 0 (ω) at a sufficient high frequency ω R . Here we apply the CC mapping only once. For our particular model [see equation (4)], an analytical expression can be obtained for J c (ω). For 4ω 2 res > γ 2 it follows that [52] J c (ω) = γω. = ω res . In terms of creation and annihilation operators defined by (7) and (8) can be written as where terms proportional to the identity and of second order in C k have been disregarded since we consider the case where the coupling to the residual reservoir is weak. Now the total Hamiltonian can be written as H = H S (t)+H CC +H SCC +H Comparing with equation (1) we have H (c) CCB .

Appendix B. Counting Field
In order to study the thermodynamics of our system we follow the full counting statistics formalism reviewed in [51]. The total Hamiltonian, including system and reservoir is H = H S (t) + H B + H SB and we assume an initial factorizing state of the form ρ tot (0) = ρ(0) ⊗ ρ B , with ρ B ∼ e −βH B . Note that we have considered the interaction Hamiltonian H SB to be time-independent. This term has a periodic time dependence in our model but after the CC mapping this time dependence only appears in the supersystem Hamiltonian and not in the interaction to the residual bath. What is presented in this section also applies for the case of replacing S by S and B by B . Let us define the modified density matrix with total (supersystem plus reservoirs) initial density matrix ρ tot (0) and modified evolution operator U (χ, t) = e −i χH B /2 U (t) e i χH B /2 , where U (t) is the evolution operator corresponding to Hamiltonian H. The variable χ is usually referred to as counting field. The evolution of operator ρ tot (χ, t) is given by Taking the trace over the reservoir degrees of freedom we define Note that the total density matrix and the reduced density matrix of our system are both recovered by setting χ = 0. The moment generating function [51]  with ∆E = E t − E 0 . It allows us to obtain the statistics of the energy transferred between system and reservoir by simple differentiation (B.5) Note that for a steady-state thermodynamic analysis we are only interested in the change in time of the energy of the reservoir given by equation (13).

Appendix C. Floquet Theory and Extended Space
Floquet's theorem establishes that, for a time-periodic Hamiltonian H S (t) = k e ik ω L t H k , with period T = 2π ω L , a solution to Schrödinger's equation is given by |ψ r (t) = e −iεrt |r(t) , where ε r are called quasienergies and |r(t) Floquet modes (states). The Floquet modes are time periodic and form a complete basis. To find them one solves the eigenvalue problem (H S (t) − i∂ t ) |r(t) = ε r |r(t) . (C.1) The periodicity of the Floquet modes allows us to map the eigenvalue problem (C.1) to a time independent one in an extended Hilbert space, sometimes referred to as Sambe space [13]. This is done by introducing an infinite dimensional space with integer quantum numbers. Its basis is given by In this space we define the operators F k and F z by their action on a state from B T L .
These operators will help us write equation (C.1) in a simpler time independent form. We also have a basis for our quantum system living in Hilbert space H S . Its basis is denoted by B S = {|φ 1 , |φ 2 , |φ 3 . . . |φ K }. Now, we can build the basis of the extended Hilbert space, which is the space in which our eigenvalue problem will become time independent.
Vectors in the extended Hilbert space will be denoted by a double ket notation |r . Its corresponding state at time t in H S will be denoted by |r(t) . Conversely, a periodic state |u(t) = |u(t + T ) is denoted in the extended Hilbert space by |u . If an operator has a periodic time dependence O(t) = k O k e iω L t , its form in extended space is defined as Furthermore, the scalar product in the extended space is defined by With the previous definitions it is now possible to write the operator Q = H(t) − i∂ t and eigenvalue problem (C.1) in extended space as where H k are the Fourier components of the Hamiltonian H(t). The reason why the eigenstates and quasienergies have been denoted with an additional index m in the extended space is because they actually carry redundant information. This becomes apparent when one considers a solution of (C.1) and defines ε rm = ε r + mω L and |r m (t) = |r(t) e imω L t such that |ψ r (t) = e −iεrt |r(t) = e −iεrmt |r m (t) . This means that the state |ψ r (t) can actually be constructed from both solutions |r m or |r m . To completely characterize the problem one has to choose the Floquet modes in the extended space whose quasienergies lie in the same Brillouin zone. From now on we will denote Floquet modes in the extended space just by |r , assuming that all lie in the same Brillouin zone.
With the Floquet modes at hand, one can find the evolution of any operator. Let us consider an arbitrary operator S.
where S ω,n = T 0 dt T r(t)|S e −inω L t |r (t) |r(0) r (0)| and ε r − ε r = ω. Depending on the form of |r(t) and S, calculating the integral in square brackets might not be trivial. In extended space this is easily calculated just by With decomposition (C.8) it is now straight forward to obtain a master equation for a driven open quantum system. For a more complete study and review of Floquet theory and driven systems the reader is referred to [14,62].

Appendix D. Floquet Master Equation
Starting from equation (B.2), going to the interaction picture and performing the standard Born and Markov approximations [27] we obtain with U 0 (t) the evolution operator associated to Hamiltonian H 0 (t) = H S (t) + H B . We take the interaction Hamiltonian to have the form H SB = S ⊗ B = S ⊗ k c k x k and define the correlation function C(χ, t) ≡ B (χ, t)B = Tr B B (χ, t)Bρ B . Using the fact that (D. 2) The correlation functions can actually be written in terms of the spectral density , with the Bose-Einstein distribution N (ω) = 1 e βω −1 . The periodicity of H S (t) allows us to decompose the system operators in the interaction picture as in equation (C.8). Performing the integrals over s and ω with the help of ∞ 0 ds e iωs = πδ(ω) + i P 1 ω and disregarding the principal value P term, one ends up with an equation in the Schrödinger picture of the form − J(∆ ω,n ) [1 + N (∆ ω,n )] Sρ(χ, t)S ω,n (t)e i∆ω,nχ − J(∆ ω,n )N (∆ ω,n )S ω,n (t)ρ(χ, t)Se −i∆ω,nχ + J(∆ ω,n ) [1 + N (∆ ω,n )] ρ(χ, t)S ω,n (t)S}, with ∆ ω,n = ω + nω L and S ω,n (t) = T 0 dt T r(t)|S e −inω L t |r (t) |r(t) r (t)| such that ω = ε r − ε r . Note that due to the periodicity of the Floquet modes |r(t) , the superoperator L(χ, t) also has the same periodicity. The heat flow is obtained by (see Sec. 2.3 and Appendix B) (D.5) Assuming that in the long-time limit the density matrix ρ(t) is time periodic with the same period as the Floquet modes we obtain equation (12). In the extended space this equation has the form with ρ a vector containing all Fourier components of ρ(t).

Appendix E. Time-Independent Interaction
We consider the case where the working medium is coupled to the cold bath and hot bath with the same system coupling operator S c = S h = σ x / √ 2ω 0 . Note that we have removed the time dependence in the coupling operator that appeared in our original model. This is done in order to compare results of the heat flows calculated with and without the CC mapping. The overall coupling strength between system and reservoir is controlled by the parameter d c in (4). Parameter γ, on the other hand, can tune how strongly peaked the spectral density is. A smaller value of γ means a more strongly peaked spectral density. Figure E1 (a) shows the heat flow of the hot reservoir as a function of the coupling strength for different values of γ. Continuous lines are the result of a standard Markovian theory, where the Born, Markov and secular approximations have been performed and dashed lines show the result of the non secular master equation of the supersystem obtained after the CC mapping. We can clearly see that the structure of the spectral density can have a strong impact on steady-state quantities, since the smaller the value of γ the bigger the difference between the standard Markovian theory and the results performing the CC mapping. The reason for this relies in the separation of timescales typically performed during the Markovian approximation. In this approximation, it is assumed that the correlation functions of the reservoir decay at a much faster rate than the typical timescales of the system. As the spectral density gets more peaked, correlation functions tend to decay at a slower rate, thus making the separation of timescales not valid. The heat flow of the cold reservoir and the work flow also follow a similar behavior as a function of the coupling strength for different values of γ. Figure E1 (b) shows the work flow as a function of γ for small couplings d c . This time we have included the results performing the CC mapping with the secular approximation (dots) for completeness. As explained in Sec. 2.2, the study of a periodically driven system with Floquet theory can make the secular approximation not suitable for analysis [60]. The secular approximation requires the systems level broadening to be much smaller than its level spacing. After the CC mapping is performed the dimension of the supersystem becomes infinite. Since for a correct treatment of the driven problem all quasienergies (see Appendix C) are mapped to the same Brillouin zone, the assumption regarding the system's level spacing looses its validity. Figure E1 (b) shows that the results of performing the CC mapping with the secular approximation might agree with the results of only performing the CC mapping in some regimes, nevertheless to fully identify this regimes independently of the model, a further study is needed. As in figure E1(a) we also see that bigger differences between the standard Markoivan theory and the CC mapping formalism are obtained for strongly peaked spectral densities. Similar results to the ones in figure E1 (a) were also obtained in [37] for the case of a non-driven heat engine.
We stress the point that even though higher values of γ show a better agreement with the standard Markovian theory after the CC mapping is performed, the coupling between supersystem and residual reservoir must remain weak in order to obtain a master equation. Since in our model the mapped spectral density is proportional to γ (see Appendix A), one needs to be careful because increasing the value of γ might make the master equation for the supersystem not valid.

Appendix F. Laser Cooling
We consider an ion interacting with a laser field [61,63]. The interaction is given by H I = d · E, where d is the dipole moment and E the electric field. We consider only two electronic states of the ion and define the electric field to be along the x-direction, H I = d S · E( r, t) = d S · E cos(kx − ωt) , (F.1) with S the spin of the ion and x its position. We also assume that the motional degree of freedom of the ion is quantized in a harmonic potential such that x = x 0 (a + a † ) and define η = kx 0 , so that If the ion is well confined we can assume that η 1. Defining Ω ≡ dE/2 we obtain H I = Ω σ x cos(ωt) − Ω η σ x (a + a † ) sin(ωt) , (F. 3) which is analogous to the time dependent part of the mapped working medium in equation (10).

Appendix G. Fluctuations
Here we analyze the fluctuations in the occupation number of the CC, defined as (∆n) 2 ≡ n 2 − n 2 , with n = a † a . For a thermal state (th) it is easy to show that these fulfil (∆n) 2 th = n th (1 + n th ) ≥ n th .  Figure G1. Left: Occupation number of the CC and its fluctuations as a function of the driving frequency ω L while keeping the resonance condition (20). Parameters used are as in figure 2. Figure G1 shows a good agreement between the variance of the occupation number of the CC and that of a thermal state, further supporting the choice of parametrisation for the effective temperature in equation (21).