Qubit absorption refrigerator at strong coupling

We demonstrate that a quantum absorption refrigerator can be realized from the smallest quantum system, a qubit, by coupling it in a non-additive (strong) manner to three heat baths. This function is un-attainable for the qubit model under the weak system-bath coupling limit, when the dissipation is additive. In an optimal design, the reservoirs are engineered and characterized by a single frequency component. We obtain then closed expressions for the cooling window and refrigeration efficiency, as well as bounds for the maximal cooling efficiency and the efficiency at maximal power. Our results agree with macroscopic designs and with three-level models for quantum absorption refrigerators, which are based on the weak system-bath coupling assumption. Beyond the optimal limit, we show with analytic calculations and numerical simulations that the cooling efficiency varies in a non-universal manner with model parameters. Our work demonstrates that strongly-coupled quantum machines can exhibit function that is un-attainable under the weak system-bath coupling assumption.

A three-level quantum system, with each transition thermalized separately with respect to an independent bath, hot (h), cold (c) and work (w). This model, which can operate as an absorption refrigerator [4,30], is not examined in the present work. (b) A two-level system with an additive interaction model, Eq. (1). This system cannot act as a QAR under the weak-coupling approximation. (c) A two-level system with an inseparable interaction model to the baths, Eq. (4). The concerted action of the baths is represented by the arrows, acting together to e.g., excite the qubit. This model is examined in the present work. When optimized, it performs as a QAR under conditions analogous to the three-level QAR of panel (a).

II. MODEL
A common design of an autonomous QAR consists of a three-level quantum system and three independent thermal reservoirs [4]. Each transition between a pair of levels is weakly coupled to only one of the three heat baths, c, h and w, where T w > T h > T c , see Fig. 1a. In the steady state limit, the (ultra-hot) work bath provides energy to the system. This allows the extraction of energy from the cold bath, to be dumped into the hot reservoir. The opposite heating process, from the hot bath to the cold, takes place as well, but it can be minimized by manipulating the frequencies of the system. The three-level QAR was discussed in details in several recent studies, see e.g. Refs. [4,30,31]. It is designed to perform optimally under the weak coupling approximation, when each bath individually couples to a different transition. Quantum coherent and strong coupling effects are expected to reduce the cooling performance of a three-level QAR.
In this paper we focus on a QAR made of a two-level system coupled to three independent thermal baths. When the baths couple to the qubit in an additive manner (Fig. 1b) we prove next that it is impossible to cool down the cold bath when the system-bath coupling is weak. By allowing for cooperative system-bath interaction between the qubit and the reservoirs (Fig. 1c), we are able to receive a cooling condition, as well as derive bounds for the maximal efficiency of the QAR and its maximal power efficiency.
A. Un-attainability of cooling for an additive dissipation model The additive model comprises a two-level system (spin, qubit) and three independent thermal reservoirs ν = c, h, w, T w > T h > T c , β ν = 1/T ν with k B = 1. The generic Hamiltonian is written aŝ Here,σ are the Pauli spin matrices.Ĥ B,ν is the Hamiltonian of the ν-th reservoir. It includes, for example, a collection of harmonic oscillators of frequencies ω j,ν ,Ĥ B,ν = j ω j,νb † j,νb j,ν with b † (b) as bosonic creation (annihilation) operators. The bath operatorÂ ν is assumed to be hermitian. It couples the bath ν to the spin, where e.g.Â ν = j λ j,ν b † j,ν +b j,ν with coupling strength λ j,ν . Assuming a factorized-product initial state, Â ν = 0 with the average performed over the initial-canonical state of the bath, weak system-bath coupling and Markovian dynamics, we obtain a second order perturbative, Markovian quantum master equation [32]. This standard Born-Markov scheme results in the stationary populations of the excited and ground state, respectively, with k d,u = ν k (ν) d,u . The decay (d) and excitation (u) rate constants k (ν) d,u , induced by the ν-th bath, depend on the details of the model. The detailed balance relation dictates local thermal equilibrium, k The energy  current, defined positive when flowing towards the qubit, can be similarly derived from the Born-Markov approximation, and it is given [33,34,36,37]. Substituting the steady state population (2) we obtain Using the detailed-balance relation and the fact that (e −β h,w ω0 − e −βcω0 ) > 0, we conclude that J c < 0 irrespective of the details of the model. Equation (3) reveals that under the additive model at weak coupling, every two reservoirs exchange energy independently. The prefactor in the denominator, k d + k u , which includes contributions from the three baths, only renormalizes the current. Since every two baths separately communicate, thermal energy always flows towards the colder bath, and a chiller performance is un-attainable. It should be pointed out that time-dependent, driven or stochastic models can realize refrigeration based on a qubit as a working medium even at weak coupling, see e.g. Refs. [38][39][40][41][42][43][44][45]. These type of driven machines are beyond the scope of our work.

B. Non-additive (strong) coupling model
It is evident that to realize a QAR with a qubit as the working substance, we must go beyond the model Hamiltonian (1), or the weak-coupling approximation. Our starting point is a revised Hamiltonian with a built-in strongcoupling characteristic, a non-additive system-bath interaction operator, Here,B ν are bath operators, assumed to be hermitian, and γ is an energy parameter characterizing the interaction energy. The non-additivity of our model is assumed to arise from a more fundamental Hamiltonian with strong interactions between the quantum system and individual reservoirs [33,34,36], see Sec. V. Non-additive models such as (4) can be also accomplished by engineering many-body Hamiltonians based on e.g. resonant conditions and selection rules.
We emphasize that our model (4) differs in a fundamental way from the QAR model analyzed theoretically in e.g. Refs. [4,46,47] and realized experimentally in a recent study [48]. In Refs. [4,[46][47][48], the working medium includes three degrees of freedom such as three harmonic oscillators, which interact via a three-body interaction term. Each oscillator is independently coupled to its own thermal bath, taken into account by introducing additive Lindblad dissipators into the time evolution equation. In contrast, in our model (4) the quantum system is as simple as it can be, a qubit. Nonlinearity is encoded into the model by assuming a non-additive interaction Hamiltonian with the three baths. This inseparability prevents us from arriving at standard perturbative quantum master equations with additive dissipators (standard multi-terminal Lindblad or Redfield).
Back to Eq. (4), we study the system's dynamics assuming a fully factorized initial state by using the Born-Markov approximation with the perturbative parameter γ. While this is analogous to a weak coupling treatment, we emphasize again that the model is defined with an inherent strong-coupling feature, the non-additivity of the interaction.

Two-bath model
Equations of motion for the spin polarization, as well as the energy current, were derived in Refs. [36,49] for the model (4) with two baths (hot and cold). The derivation relies on the assumption B ν = 0, which could be satisfied exactly or under conditions such as strong coupling or high temperature [36]. Further, this assumption can be relaxed by re-defining the model Hamiltonian to add and subtract the thermal average of the interaction Hamiltonian, rediagonalizing then the system's Hamiltonian and proceeding with the perturbative treatment [50,51]. The population dynamics [33,34,36] satisfiesṗ with rate constants Here, M ν (t) = B ν (t)B ν (0) is the two-time correlation function with the average performed with respect to the canonical (initial) state of the ν thermal bath. In Fourier space In what follows, we refer to this function as the "Fourier bath-correlation function" (FBCF). This function is real valued and positive. In our work, the FBCF has a physical dimension of inverse energy ( = 1). Formally similar to P (E) theory [35], the detailed-balance condition is satisfied for the individual components, Mν (ω) Mν (−ω) = e βν ω , but we do not have such a relation for the convoluted rate constant M (±ω 0 ). Within the same treatment, the thermal energy current, flowing from the cold bath to the system, is given by a rather intuitive expression [36,49], Here, we have absorbed a factor C 2 = 1 2π γ 2 2 in the definition of the current. The heat current exhibits cooperative energy transfer processes: The first term describes contributions to the current due to the decay of the qubit, with ω energy exchanged with the cold bath and the rest ω 0 − ω absorbed or released by the hot bath. A similar reasoning applies to the second term in Eq. (7), which describes the excitation of the qubit. For simplicity, in what follows we eliminate the spin gap and take ω 0 = 0. We then immediately conclude that p e = p g = 1/2, thus the heat current (7) simplifies to Eq. (7) can be derived from a full counting statistics perspective, and the resulting cumulant generating function satisfies the steady state heat exchange fluctuation theorem [36,[49][50][51]. Therefore, our framework satisfies the second law of thermodynamics, which is a direct consequence of the exchange fluctuation relation. From this, it is easy to prove that if β c > β h , J c < 0, meaning that the heat current flows towards the cold bath.

Three-bath case
Back to the QAR model Hamiltonian (4), a full-counting statistics analysis allows us to describe energy exchange with three thermal reservoirs [49], in a complete analogy to the two-bath case described in Sec. II B 1. This formalism directly provides both population dynamics and the dynamics of the so-called cumulant generating function, handing over all current cumulants. The population dynamics follows Eq. (5), with the FBCF now however given by The energy current, from the cold bath towards the qubit, is given by Again, we have absorbed a prefactor C 3 = 1 (2π) 2 γ 2 2 into the definition of the current. This intuitive expression, which can be also suggested phenomenologically [52], describes coordinated three-bath energy exchange processes, with an overall conservation of energy. An amount of energy ω 1 is delivered to the cold bath or absorbed from it, while the other reservoirs assist by providing or absorbing the rest of the energy, so as to complete a decay (first integral) or an excitation (second line) process. When ω 0 = 0, we receive While we write here an expression for J c only, the full-counting statistics approach [49] readily hands over analogous expressions for J w and J h . Particularly, J h is received from Eq. (10) by interchanging the 'c' and 'h' indices, and J w then follows from energy conservation.

A. Cooling window and efficiency
Can we realize a QAR based on a qubit, using the model (4)? Our objective is to derive a cooling condition from Eq. (11), i.e., find out whether we can engineer the system and the baths to achieve refrigeration, J c ≥ 0.
The FBCF M ν (ω) is related to the spectral density of the ν thermal bath, see Sec. V. Let us assume that the baths are engineered such that these functions are characterized by the frequency θ ν > 0, satisfying a resonant assumption We now analyze the performance of our system as a QAR in an idealized limit, then under more practical settings.
The resonant assumption will be assumed throughout, though it is not a necessary condition for refrigeration in non-ideal designs as we demonstrate through simulations and discuss in our conclusions.

Ideal design
In an optimal design, the FBCF M ν (ω) is restricted to be nonzero within a narrow spectral window. In the most extreme case, we filter all frequency components besides the central mode, Here, ǫ ν > 0 is a dimensionless parameter. We evaluate the integrals in Eq. (11), and find that there are only two contributions to the cooling current, when ω 1 ≤ 0 and ω 2 ≥ 0, and the other way around. We then receive the cooling condition, J c ≥ 0, The first term describes the removal of heat from the cold bath-assisted by the work reservoir-and its release to the hot bath. The second term accounts for the reverse process, with energy absorbed from the hot bath and released into both the cold and work reservoirs. Equation (14) can be also organized as follows, which precisely corresponds to the cooling condition as obtained for a three-level QAR-analyzed with a Markovian master equation with additive dissipators [4,30].

Non-ideal design
Let us now consider a more physical-and less optimal model. We engineer the FBCFs as follows (ω > 0), For negative frequencies, the detailed balance relation multiplies each function by a thermal factor. Here, H(x) is the Heaviside step function, δ(x) is the Dirac Delta function and ǫ ν is a dimensionless parameter. As we show below, these parameters, which characterize the reservoirs' spectral functions, neither affect the cooling window nor the efficiency. We assume the resonant assumption (12) to hold, and that the width parameter δ is small, δ ≪ θ ν .
We calculate the cooling current based on Eq. (11) by breaking it into four contributions. Let us first inspect the ω 1 ≥ 0 and ω 2 ≥ 0 term, Since θ c + θ h − δ > θ w + δ, the integral collapses to null. A similar argument brings a zero contribution from the negative branch, J (−,−) c = 0, which includes ω 1 ≤ 0 and ω 2 ≤ 0. We proceed and evaluate the contribution to the cooling current from We had utilized the resonant assumption (12) in the last line. A similar analysis gives Putting together Eqs. (18) and (19), we organize the cooling condition, J c ≥ 0, as We Taylor-expand this result to the first nontrivial order in δ, which is the third order, and receive Using Re-arranging this expression, we get a cooling condition In the limit δ → 0, we retrieve Eq. (14), which corresponds to the three-level QAR analyzed with a Markovian master equation with additive dissipators [4,30]. It is important to recognize that the δ dependence in Eq. (23) is non-universal, and it depends on our choice (16). For example, using a different setting, with M c (ω) and M h (ω) as step functions of width 2δ, but M w (ω) a Dirac delta function, we derive an alternative cooling condition, The role of the width parameter δ is non-trivial, and it can reduce or increase the cooling window. Arriving at equation (15), and receiving its generalizations to finite width, Eqs. (23) and (24), are central results of our work. The cooling window depends on δ in a non-universal manner. In the optimal limit, δ → 0, we recover the three-level weak-coupling condition, which is bounded by the Carnot limit of a macroscopic absorption refrigerator, as we show next.
The efficiency of a refrigerator is defined by the coefficient of performance (COP), the ratio between the heat current removed from the cold bath J c and the input heat from the work bath J w , η ≡ Jc Jw . For convenience, in what follows, we sometimes refer to the cooling COP as the "efficiency" of the refrigerator, in the sense that this measure characterizes the competence of the machine, but remind the reader that it can e.g. assume values larger than one.
It is convenient to write down an equation for J h , analogous to Eq. (11), and evaluate it with the model (16). The heat current from the work bath is given by J w = −J c − J h . The currents from the work and cold baths are We expand the currents to third order in δ and receive the cooling COP The width parameter δ could affect the efficiency in a non-monotonic way. However, to the lowest (linear) order in δ we gather The cooling condition (15) can be used to derive a bound on efficiency. In the ideal limit, the cooling window is defined from We inverse this expression and reach which can be also expressed as finally receiving θc θw ≤ Tw −T h T h −Tc

Tc
Tw . Comparing this expression to Eq. (27), we gain a bound on the cooling COP (in the δ → 0 limit), which is nothing but the Carnot bound. Specifically, when T w ≫ T c , T h , we find that θc θw < T c /(T h − T c ), which is the Carnot bound for cooling machines. We emphasize that Eq. (31) was received in previous studies, see e.g. Ref. [30], yet restricted to quantum systems that evolve with additive dissipators. The COP of non-ideal models with a finite width M ν (ω), can be also calculated and simulated, as we do in Sec. IV. It is important to note that the COP of our model is always limited by the Carnot bound (31) as was recently proven in Ref. [12] based on the heat exchange fluctuation theorem.
The analysis presented in this subsection can be extended to describe a non-degenerate (biased) qubit ω 0 = 0 as we show in Appendix A. We find that for small ω 0 we can adjust the functions M ν (ω) so as to recover the ideal cooling window (14) and the corresponding efficiency bounds.
Another natural question concerns the potential to realize a QAR with continuous functions M ν (ω), rather than the hard-cutoff model (16). In Appendix B we calculate the cooling current (11) using Gaussian functions for M ν (ω), and prove that the system cannot act as a QAR. It would be interesting to find whether there is a general proof for the un-attainability of cooling once the functions M ν (ω) fully overlap.

B. Efficiency at maximum power
The discussion above concerns a bound for the maximal COP-obtained for vanishing cooling power. For practical purposes, however, one needs to operate machines at non-vanishing output power. A more practical figure-of-merit is the maximal power efficiency (MPE) η * , meaning, that we calculate the cooling COP when maximizing the cooling current. An intriguing question, which recently received much attention, is whether η * can approach η c -albeit at finite cooling power. For heat engines, bounds for the MPE were extensively investigated [53][54][55][56][57][58][59]. The MPE of refrigerators was examined in e.g., Ref. [30], demonstrating that it is constrained by the spectral properties of the thermal reservoirs. Nevertheless, Ref. [30] was concerned with additive (weakly-coupled) system-bath models.
To derive an expression for the efficiency at maximal cooling power, we go back to Eq. (25), and analyze it to the lowest order in δ, This function is linear in θ c for small values, but it decays exponentially with θ c beyond that. We thus search for θ c which maximizes J c (θ c ) by solving ∂Jc ∂θc = 0. We find that ln [ Since ln(1 + x) ≤ x, we turn the equality (33) into an inequality, or The MPE is thus bounded by Since β w ≤ β h , we receive the bound Interestingly, the efficiency at maximum cooling power is upper-bounded by half of the Carnot bound for cooling. The second order term in the expansion is − 1 4 η 2 c , which is different than the value obtained for low dissipation engines [57].
We repeat this exercise using functions M ν (ω) of a finite width δ. For example, we start from equation (25), expand it to third order in δ, and solve ∂Jc ∂θc = 0 to obtain a relation for θ c which maximizes the cooling current. Using again the fact that e x ≥ 1 + x, we receive The left hand side in this inequality corresponds to the COP, even at large δ, as long as we work in the "linear response limit" of β h → β w → 0, see Eq. (26). The right hand side then collapses to η c /2. We therefore conclude that in the linear response limit the maximum power efficiency is upper bounded by η c /2, independent of the details of the model.

IV. SIMULATIONS
In this section, we examine the performance of the qubit-QAR beyond the ideal (or close to ideal) limit as explored in Sec. III. We perform simulation directly by using Eqs. (10) and the analogous expression for J w . For the FBCF we use the following model (ω ≥ 0), Here, δ ν is not necessarily small, n ν (ω) = [exp(β ν ω) − 1] −1 is the Bose-Einstein distribution function, H(x) is the Heaviside step function, and p is a dimensionless parameter. The same form holds for negative frequencies, but with a detailed balance factor such that Mν (ω) Mν (−ω) = e βν ω . The model is motivated from physical grounds in Sec. V. In simulations we use p = 1, thus one should add a prefactor with units of an inverse energy to (39). This prefactor is taken as unity here. We found that results did not qualitatively change when using other powers p. We assume the resonance assumption (12) to hold, though this is not a necessary condition for realizing a qubit-QAR once we operate the system beyond the ideal limit (large δ ν ).
As a concrete example, we assume that the cold bath has a finite-fixed width of δ c = 0.2, while the other two windows are tuned with δ w = δ and δ h = 2δ. For a schematic representation, see Fig. 2. In Fig. 3, we display the cooling current and the cooling COP for the model of Fig. 2. The current J c /C 3 is dimensionless, with C 3 ∝ γ 2 , γ as the system-bath interaction energy, see text below Eq. (10). To receive the current in physical units of e.g. J/sec, one should further divide it by . We find that the system can act as a refrigerator within a certain domain: According to the ideal case, Eq. (15), cooling takes place when β h ≥ θc θ h (β c − β w ) + β w . In our parameters, this reduces to β h ≥ 0.9/3 + 0.1 = 0.4. This estimate qualitatively agrees with simulations at small δ.
The width parameter is important as well, and we find that refrigeration takes place as long as δ < ∼ 1.6. Note that when δ ≥ 1.8, the functions M c (ω) and M w (ω) overlap leading to energy leakage directly from the work bath to the cold environment. Nevertheless, it is interesting to note that the hot and work reservoirs already touch when δ ≥ 2/3, yet their overlap does not prohibit the cooling process. Overall, temperature, θ ν , and the width parameters δ ν intermingle in the expression for cooling. As a result, the cooling window depends on δ in a non-trivial manner. As an additional comment, in Fig. 3, we consistently receive η < θ c /θ w = 0.5. However, when using a narrower window for the cold bath, e.g., δ c ∼ 0.1, we find that in fact the cooling COP can exceed 1/2, yet obviously it still lies below the Carnot bound.
In Fig. 4, we display slices of the cooling current as a function of the width δ and temperature β h , taken from the contour plot, Fig. 3. For small value of δ the cooling current grows linearly. However, at a certain point (δ = 1.3), when the cold and work reservoirs almost touch, the cooling current begins to drop with δ, eventually missing the cooling operation altogether. The behavior of the cooling current with β h is monotonic.
The role of a nonzero gap ω 0 = 0 on the operation of the qubit-QAR is presented in Fig. 5, using again parameters as in Fig. 2. In panel (a), we set δ = 1. We find that within the cooling window, the performance of the QAR is intact for small ω 0 , but it deteriorates and eventually disappears once ω 0 is comparable to differences of the central frequencies, θ ν . Outside the cooling window, for large δ, in panel (b) we reveal a non-trivial non-monotonic behavior of J c with ω 0 , with a large value of ω 0 manifesting a cooling function that is missing in the degenerate case. Overall, as long as ω 0 ≪ δ, it plays an insignificant role in the refrigeration behavior. Beyond that, it can introduce cooling, enhance or reduce performance, depending on the particular choice of parameters.

V. PHYSICAL MODEL
Before concluding, we discuss physical models corresponding to the model Hamiltonian (4) and the functions M ν (ω). We assume that the baths include collections of harmonic oscillators,Ĥ B,ν = j ω j,νb † j,νb j,ν , and that bath operators, which are coupled to the system, are of the form The two-time correlation functions are thus given by with n ν (ω) = [e βν ω − 1] −1 as the Bose-Einstein distribution function. In the frequency domain, We now define the spectral density function, which is related to the density of states of the bath, and find that, with g ν (ω) analytically continued to the complete real axis as an odd function. According to Eq. (44), the function M ν (ω) is linearly related to the spectral density function g ν (ω) of the respective bath. Thus, if we can engineer the density of states of the bath, or filter its frequencies, we can realize a structured model of the form (16).
Let us now also comment on the relation of the Hamiltonian (4) to the nonequilibrium (multi-bath) spin-boson Hamiltonian. We begin with the additive model, and transform it to the displaced bath-oscillators basis using the small polaron transformation [60],Ĥ p =Û †ĤÛ , Here,σ ± = 1 2 (σ x ± iσ y ) are the auxiliary Pauli matrices,Ω = νΩ ν ,Ω ν = 2i j λj,ν ωj (b † j,ν −b j,ν ). It can be shown that under the so-called noninteracting blip approximation [61][62][63], the population dynamics follows Eq. (5) with the time correlation function The average is performed with respect to the initial, product thermal state of the baths. The function Q(t) = ν Q ν (t) is complex, with real and imaginary components, The multi-bath spin-boson example clearly demonstrates that non-additivity of the interaction Hamiltonian embodies strong coupling effects. Using Eq. (46) as a starting point, a second order time-evolution scheme with respect to ∆ provides an equation of motion for the population dynamics in the form (5), which is beyond second order in the system-bath interaction energy λ [36]. Nevertheless, to realize a QAR we need to design the function M ν (t), which is the central object in the expressions for the population and energy current. For the spin-boson model, this function, [see Eqs. (47)- (48)], depends in a non-linear manner on the reservoirs' density of states, thus it is not immediately clear how to physically-rationally engineer it to receive the form (16).

VI. SUMMARY AND OUTLOOK
What is the smallest possible fridge? Using a quantum master equation in Lindblad form with additive dissipators, it was argued in Ref. [64] that a three-level system (qutrit) is the smallest refrigerator, if each transition is thermalized independently. Here, we demonstrate that in fact a qubit could serve as smallest refrigerator-if the three thermal reservoirs c, h and w couple in a non-additive manner to the qubit, to transfer energy in a cooperative manner. Our formalism is consistent with classical thermodynamics, as it complies with the heat exchange fluctuation theorem.
In the ideal limit we receive refrigeration if (i) the reservoirs are engineered, and are allowed to exchange energy within a highly restricted spectral window, (ii) a resonant condition θ c + θ w = θ h is satisfied for the characteristic frequencies of the three reservoirs. In this special situation, we derived closed expressions for the cooling window, the efficiency of the refrigerator (characterized by its COP), the maximal efficiency (bounded by the Carnot limit) and the maximum-power efficiency. We found that a qubit-QAR with a non-additive dissipation form can embody function that is prohibited under the weak system-bath coupling assumption. We had further studied, analytically and numerically, the operation of the system beyond the ideal limit and found that it can operate as a QAR for a fair range of parameters.
Throughout this study, we had assumed a resonance condition (12) for the characteristic frequencies of the three reservoirs, so as to optimize performance. Nevertheless, it is imperative to realize that one could construct a qubit-QAR without structuring the work reservoir. Going back to Eq. (11), we evaluate the cooling current when the cold and hot FBCF take the form of a Dirac delta function as in Eq. (13), with θ h > θ c > 0, but the work reservoir's FBCF is simply a constant, M w (ω) = ǫ w for ω ≥ 0. We then obtain a modified cooling window (compare to Eq. (14)), The new, last two terms correspond to the extraction of energy from the cold bath, assisted by the hot bath, to be dumped into the work reservoir, and the reversed process. We can reorganize the cooling condition as where A ≡ e −βw(θc+θ h ) e β h θ h − e −βcθc . After some manipulations, we arrive at Since A > 0, it is clear that the cooling window is reduced here relative to the case with a structured work bath, (15). This new situation is significant: We do not place here stringent conditions on the work reservoir, which could in fact be featureless and still support cooling. Furthermore, this scenario, where the work reservoir is allowed to provide or absorb both low and high frequencies, is obviously un-accounted for in the traditional three-level weak-coupling setting where each bath excites a particular transition in a resonant manner (see Fig. 1)(a). A strongly-coupled qubit-QAR thus offers new regimes of operation, missing in multi-level, weakly coupled designs. We examined an additive dissipation model and proved that it cannot support a QAR performance based on a qubit. In contrast, a non-additive model, with the three reservoirs acting in a concerted manner, can achieve refrigeration. It is of an interest to examine an in-between model, which could be more practical. For example, the cold and work reservoirs could couple strongly to build up a non-additive dissipator, but the hot bath would weakly-separately couple to the system, bringing in an additional, local dissipation term. This scenario could be treated using the reaction coordinate method [23] or the polaron-transformed master equation [50,51,65,66], which can smoothly interpolate between additive and non-additive dissipation problems.
The description of quantum systems that are strongly coupled to multiple thermal reservoirs poses a significant theoretical challenge. Markovian master equations of Lindblad form with additive (local) dissipators provide a consistent thermodynamical description of observables [3]. However, it is not yet established how to formulate quantum thermodynamics beyond the weak coupling limit. Our work here exemplifies a strong coupling framework which is thermodynamically consistent. It is of interest to compare our approach and predictions to other recent studies on strongly-coupled energy conversion devices [20][21][22][23][25][26][27], and further develop numerically exact techniques that could target such problems. Understanding the role of both strong coupling and non-Markovianity of the reservoirs on the operation of driven and autonomous thermal machines remains a challenge for future work. financial support by the DFG (SCHA 1646/3-1, GRK 1558). Discussions with Javier Cerrillo are gratefully acknowledged.

Appendix A: Cooling condition for a non-degenerate Qubit-QAR
We investigate here the cooling performance of a qubit of a finite energy gap ω 0 . For simplicity, we assume the following functions (ω > 0). The negative-frequency branch is decorated by detailed-balance (thermal) factors, Here, ǫ c,h,w are dimensionless parameters. Since these parameters do not influence the cooling window and the cooling efficiency, we ignore them below. We maintain the resonant condition, θ c + θ w = θ h . We further assume that ω 0 is much smaller than θ ν , ν = c, h, w.
One should note that the steady-state populations of the qubit are no-longer equal. We will begin by evaluating them from Eq. (5). The three-bath convoluted decay rate is We break the integral into four contributions. For ω 1 ≥ 0, ω 2 ≥ 0, since we assume that ω 0 is sufficiently small, such that θ c + θ h − 2ω 0 > θ w . Similarly, when ω 1 ≤ 0, ω 2 ≤ 0, We receive a finite contribution when ω 1 ≥ 0 and ω 2 ≤ 0, as well as for ω 1 ≤ 0, and ω 2 > 0, We now turn to the expression for the cooling current, Eq. (10). It is easy to calculate the integrals in this equation since there is only one extra ω 1 term, which is trivial to account for under (A1). We evaluate the two integrals, and In total, the cooling current is The cooling window, J c ≥ 0, precisely follows Eq. (14)-obtained for the case ω 0 = 0.