{'fraction_non_alphanumeric': 0.046113662973971, 'fraction_numerical': 0.013029511543426702, 'mean_word_length': 4.4065898393118, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cooling condition for multilevel quantum absorption refrigerators

Models for quantum absorption refrigerators serve as test beds for exploring concepts and developing methods in quantum thermodynamics. Here, we depart from the minimal, ideal design and consider a generic multilevel model for a quantum absorption refrigerator, which potentially suffers from lossy processes. Based on a full-counting statistics approach, we derive a formal cooling condition for the refrigerator, which can be feasibly evaluated analytically and numerically. We exemplify our approach on a three-level model for a quantum absorption refrigerator that suffers from different forms of non-ideality (heat leakage, competition between different cooling pathways), and examine the cooling current with different designs. This study assists in identifying the cooling window of imperfect thermal machines.


I. INTRODUCTION
Recent years have seen an explosion in research aiming to tie the fields of quantum mechanics and classical thermodynamics, with a focus on quantum thermal machines [1][2][3] . Of particular interest is the quantum absorption refrigerator (QAR) 4-12 , a machine, analogous to the classical counterpart, which employs a heat source rather than an external power input to achieve refrigeration of a target component. Here, the "working fluid" that converts energy to refrigerate is quantum mechanical in nature. QARs are particularly useful in nanotechnological applications because they utilize waste heat and operate autonomously without driving. Potential applications include quantum state preparation, quantum computing 13 , and biological function in proteins 14 .
Different mechanisms are responsible for irreversibility in QARs, impacting performance bounds. This includes heat leaks: the parasitic coupling of heat baths to the cooling and driving transitions, internal dissipation, which corresponds to the competition between different cooling pathways, and delocalized dissipation. The effect of such lossy mechanisms on the performance of QARs was discussed in specific model systems, see e.g. Refs. 8,36 . A graph theory analysis of multistate QARs (with their dynamics described by classical rate equations) was re-cently presented in Ref. 10 . This treatment allows decomposition of the cooling current in terms of the circuits of the machine, bringing in strategies to enhance the cooling performance. Most notably, it was demonstrated in Ref. 10 that the performance of incoherent multilevel QAR was smaller than or equal to the performance of the best-performing circuit component. A graph theory treatment expresses the cooling current as a decomposition of different cycles. Nevertheless, a framework for efficiently resolving the cooling window, cooling current and associated noise of a QAR device directly from the generator of the dynamics is highly desirable. Particularly, given efforts to realize QARs, a viable analysis of the performance of imperfect devices suffering lossy processes is crucial.
In this work, we show that a cooling condition for QARs can be defined quite generally and analytically for a multilevel system coupled to multiple heat baths, where one of the baths (labeled C) is refrigerated by the setup. The flexibility of our analytical results allows us to demonstrate that a three-level QAR 4,5,7 can operate even in the presence of heat leaks and internal dissipation to all baths in the setup. More broadly, our formalism allows us to efficiently calculate, analytically and numerically, heat exchange in devices with multiplexed couplings to thermal reservoirs using linear algebra operations applied directly on the Liouvillian of the dynamics. These results are achieved using an open quantum system full-counting statistics approach.
The paper is organized as follows. In Section II, we introduce the non-ideal multilevel refrigerator and derive expressions for its cooling current and the cooling condition using a truncated cumulant generating function. Theoretical calculations are exemplified with simulations of a three-level QAR in Sec. III. Further examples and details are delegated to Appendixes A and B. We conclude and discuss future directions in Section IV.

A. Setup and equations of motion
Our setup includes an N -level, nondegenerate system, which is coupled to multiple heat baths (enumerated by µ) maintained at different temperatures. These baths are responsible for inducing transitions (excitations and relaxations) between levels, N ≥ 2. We specifically identify the coldest heat bath (denoted by C); our goal is to extract heat from this reservoir, assisted by other baths, and dump the heat into other heat reservoirs. Note that only thermal energy (heat) is exchanged here between system and baths. A three-level example is displayed in Figure (1).
We are interested in the steady state behavior of this nonequilibrium system. To this end, we make several standard assumptions: weak system-bath coupling, factorized system-bath initial state, Markovian dynamics, and the secular approximation (decoupling population and coherences). The latter approximation is valid when energy levels of the system are nondegenerate and the temperature is high. Under these assumptions, the dynamics of the N -level system can be organized as a Markovian quantum master equation (MQME), an equation of motion for the system's population, p 42 , Here, L is a Liouvillian matrix of size N by N , which is built from bath-induced rate constants that lead to changes in the system's state populations. Under the weak-coupling approximation, in additive interaction models L is given by a sum of the Liouvillians (transitions rates) for each bath, µ L µ .
FIG. 1: A scheme of the model refrigerator. The threelevel system acts as a working material. Transitions between levels are induced by the heat baths. Unlike the optimal design 6,9 , several reservoirs may simultaneously couple to the same transition.
Given this out-of-equilibrium, N -level, multi-bath model, our objective is to calculate the heat current extracted from the cold bath. We approach this challenge using a full-counting statistics (FCS) formalism [43][44][45] .
This tool provides the cumulant generating function (CGF) for heat exchange, thus fully characterizing energy exchange in steady state. The FCS treatment is appealing for several reasons: (i) It allows validation of the fluctuation symmetry, thus ensuring the thermodynamical consistency of approximate techniques. (ii) An FCS treatment automatically hands over the averaged current and high order cumulants. (iii) Practically, extending the Markovian population dynamics, Eq. (1), to report on the FCS of heat exchange is quite straightforward, as was exemplified in numerous studies 43,46,47 .
To derive the characteristic function for heat exchange, we follow the two-time measurement protocol 43 and introduce a counting parameter χ which tracks thermal energy entering/leaving the cold reservoir. According to our sign convention, heat exchange (heat current) is positive when it flows towards the central system. In principle, one could keep track of energy transfer at each boundary with the system by introducing a counting parameter χ µ for each reservoir µ 46 . However, to simplify our analysis we limit it to a single counting parameter.
In the MQME approach, the FCS is obtained by writing down the counting field-dependent quantum master equation,ṗ(t, χ) = L(χ) p(t, χ), with counting-field dependent Liouvillian and population. The microscopic derivation of this equation was detailed e.g. in Ref. 46 , and exemplified in Appendix A and B for two-and threelevel models, respectively.
The characteristic function for heat exchange is given by Z(t, χ) = I|p(t, χ) , where I| = (1, 1, 1, ...) T is the identity vector. The CGF for steady state heat exchange is obtained in the long time limit as Differentiating this function k times with respect to (iχ) brings the kth cumulant of the process. Specifically, the first cumulant is the energy current between the system and the reservoir where counting is performed, The second cumulant is the noise of that current. We obtain the CGF by diagonalizing L(χ) and selecting the eigenvalue corresponding to the long-time limit. This is the eigenvalue with the smallest magnitude for its real part, λ min. real magnitude = G(χ). Formally, we solve the eigenvalue problem by writing down the characteristic polynomial for the matrix L(χ) in terms of its eigenvalues, λ 1,2,...,N , The coefficients a j (χ) are functions of the counting parameter. In what follows, we distinguish between a j (χ) and a j (0), which are the coefficients of the characteristic polynomials of L(χ) and L(0) = L, respectively. We now list several important properties of these elements.
The coefficients of the characteristic polynomial can be expressed in terms of the eigenvalues: a j (χ) of λ N −j is given by the sum of all products of j eigenvalues 48 , Specifically, For the physical problems that we consider, the trace, a 1 (χ) does not depend on the counting field, see Appendices A and B. As for the other coefficients, from the general rule we get e.g. that for a 3 × 3 matrix, a 2 (χ) = λ 1 λ 2 +λ 1 λ 3 +λ 2 λ 3 . The last property in Eq. (6) stems from the fact that under our dynamics, the system reaches a unique steady state, thus L(0) must acquire a zero eigenvalue, which nullifies the determinant. From the detailed balance condition it can be shown that the eigenvalues of L(0) must be real since the eigenvalue problem can be mapped to a real symmetric matrix. Furthermore, positive values for the eigenvalues are unphysical since that would lead to exponentially growing probabilities. 42,49 . Therefore, the eigenvalues of the rate matrix L(0) are real, negative numbers. With that in mind, based on Eq. (5), we conclude that the coefficients a 1,2,...N (0) must be positive. In particular, we later use the fact that a N −1 (0) > 0.

B. Current and cooling condition
For an N × N problem, L(χ) has N eigenvalues, with G(χ) as the eigenvalue with the smallest-magnitude real part. As discussed in Ref. 11 , the eigenvalue problem can be greatly simplified if we focus on the heat current only, Eq. (3): Because G(χ = 0) = 0, Eq. (4) can be truncated to first order in λ -and would still solve exactly for the first cumulant. More generally, Eq. (4) can be truncated to the order of the desired cumulant.
Truncating the characteristic polynomial (4) to first order we get To obtain the heat current we differentiate Eq. (7) with respect to (iχ). Using a N (0) = 0, we arrive at We recall that a N (χ) = (−1) N det [L(χ)]. Applying Jacobi's formula, which connects the derivative of the determinant to a trace relation, we get an analytic expression for the heat current, from the reservoir where counting is performed towards the quantum system, Here, adj [L] is the adjugate matrix of the Liouvillian; the transpose of the cofactor matrix of L. If counting is performed on the cold bath, Eq. (9) provides the cooling current J C . However, the equation can be readily applied to calculate other currents by assigning the counting parameter to the respective heat source. We now write down a general condition on the cooling window. Since a N −1 (0) > 0, cooling is achieved if Equations (9)-(10) are the main results of this paper. They offer the following benefits over standard calculations of the heat current 50 : (i) Equation (9) allows us to feasibly calculate heat currents analytically for systems described by the generic evolution, Eq. (1), without a prior knowledge of the steady state population. Similarly, the cooling condition (10) can be evaluated analytically to easily determine if a quantum system is suitable to act as a refrigerator.
(ii) In both Equations (9) and (10) only the diagonal elements of the matrix product are required to perform the trace operation. The matrix ∂L(χ)/∂(iχ)| χ=0 is sparse in typical models, thus drastically reducing the complexity of the calculation. Importantly, this matrix can be readily-intuitively constructed, and it does not require knowledge of the principles of the FCS formalism: As we illustrate in the Appendices, when calculating J ν the terms surviving in ∂L(χ)/∂(iχ)| χ=0 are energy transfer rates to/from the reservoir ν and the system.
(iii) Eq. (9) can be used to study heat exchange at each contact, thus providing the cooling coefficient of performance, that is the ratio of cooling current to input power extracted from the work reservoir.
(iv) Our derivation for the cooling current does not discriminate between different sources of irreversibility: whether these are due to the competition between different cooling pathways ("internal dissipation") or due to heat leaks with multiple reservoirs coupled to the driving (work) transition.
(v) Our derivation does not rely on the additivity of the total Liouvillian with the different baths. Therefore, expressions (9)-(10) can be applied to more general setups, such as the non-additive system-bath models discussed in Ref. 46 .
Appendix A exemplifies the evaluation of the heat current in a two-level two-bath setup. In Appendix B, we construct three-level three-bath models and illustrate the calculation of the cooling current and cooling condition in ideal and non-ideal QARs.

C. Noise power
It is possible to obtain the noise power of a specific cur- by following a similar procedure to the calculation of the heat current. To achieve that, we use a characteristic polynomial truncated to second order in λ, a N −2 (χ)λ 2 + a N −1 (χ)λ + a N (χ) = 0, which yields the CGF, The resulting noise expression [after a second derivative in (iχ)] is cumbersome. Therefore, we only consider classes of models for which both a N −1 (χ) and a N −2 (χ) do not depend on the counting parameter. This is the case, for example, for a three level system if we further assume that the cold bath is coupled to a specific transition that no other bath is coupled to 51 . Under this approximation, the noise can be received via For a given device setup, the noise could be minimized with respect to tunable parameters to determine the fundamental limit of noise reduction. While this expression seems cumbersome, the extreme sparsity of the differentiated Liouvillian matrices drastically simplifies calculations.

III. CASE STUDY: THREE-LEVEL QAR WITH COMPETING CYCLES AND LEAKS
A QAR is a continuous-cycle heat machine, which operates without external driving. It extracts energy from a cold (C) reservoir assisted by heat supplied from a socalled work reservoir (W ), and it dumps the heat into a hot reservoir (H). The simplest version of such a device, the three-level quantum absorption refrigerator 4 , has been explored in detail in the weak system-bath coupling limit 5,7 .
We exemplify the analytic results derived in Sec. (II) on a general, non-ideal three-level QAR. The setup includes three levels, |1 , |2 , and |3 . Transitions between states |i and |j are induced by three heat baths, with β C > β H > β W ; β µ is the inverse temperature of the µ bath. The Hamiltonian of the setup is given by, Here,â † µ,k (â µ,k ) is a bosonic creation (annihilation) operator of mode k in the µ bath.B µ ij = k u µ,k ij â † µ,k +â µ,k is a bath operator that couples to the i ↔ j transition with energy u µ ij , which is assumed to be a real number.
The Hamiltonian (13) is a general version of the ideal 3-level quantum absorption refrigerator, in which only the u C 12 , u W 23 and u H 13 transition elements are nonzero; this scenario is 'ideal' because tuning the system's parameters (energy levels) can achieve the Carnot cooling limit, see Appendix B and e.g. Refs. 9,11 .
In what follows, we study the cooling performance of non-ideal setups. Specifically, in our work, the term 'heat leaks' refers to having several baths coupled to the same transition. This leads to a competition between the different baths, to either excite or relax population. Internal dissipation arrises when cooling can be achieved by more than one cycle 37 .
We can feasibly predict the parameters for the cooling window by writing down the Liouvillian L(χ) = L C (χ) + L W + L H . For harmonic baths and bilinear system-bath couplings the Liouvillians are made from the rate constants, Here, n µ (ω) = [e βµω − 1] −1 is the Bose-Einstein occupation function. Γ µ ij (ω) = 2π k (u µ,k ij ) 2 δ(ω − ω µ,k ) is the spectral density function of the baths, E i,j = E i − E j . Specifically, we employ an ohmic model, Γ µ ij (ω) = γ µ ij ωe −|ω|/ωc with γ as a dimensionless coupling parameter and ω c the cutoff frequency, assumed to be high. Since we study the problem in the weak-coupling limit, the full functional form of the spectral function carries no impact on the cooling current; only the value in resonance with the system's transitions is used in calculations.
We examine different configurations of the QAR model in the four panels of Fig. (2). In each case, we indicate schematically the allowed transitions as activated by the different baths using colored arrows: blue for the C bath, red for H, and green for W . The relative width of the arrows indicate the relative strength of that bath-induced transition. The red curve encloses the region where cooling is predicted according to Eq. (10), and the contours show the positive cooling current calculated numerically. For visibility, negative current (heat flow directed towards the cold bath) is not presented (gray region). As can be seen, the cooling condition (red curve) perfectly predicts the cooling window for each of the QAR setups.
Our analysis demonstrates that the current at the cold contact can be generally written as Here, F i,j cyc identifies a circuit that can realize a QAR with the cold bath coupled to the i ↔ j transition and the work and hot baths completing a cycle, e.g. i − → i, and the reversed process. 'Cycle' contributions must involve the three baths. The second term, F i,j leak,µ , is always negative. It comprises heat leaks from the µ bath to the cold bath when both reservoirs are coupled to the i ↔ j transition. Thus, leak terms involve two reservoirs, either H or W , and C. This classification (see Appendix B) agrees with 10 . We now use Eq. (15) to explain the results of Fig. 2.
Model A. The ideal QAR setup is constructed by taking nonzero values for γ C 12 , γ W 23 , γ H 13 ; all other couplings are null. As we show in panel a, this setup achieves the widest cooling window and the most substantial cooling current. We note that there is an optimal level spacing E 2,1 for maximizing the cooling current, which depends on the temperature of the cold bath. This setup is termed 'ideal' since it can achieve the Carnot bound. Increasing β H (decreasing T H ) is beneficial for cooling: When the temperature gradient between H and C is reduced, heat flow from H to C is minimized. In Appendix B we use Eq. (10) and derive the cooling condition, F 2,1 cyc > 0, or explicitly It agrees with previous studies, see e.g. 9 . Next, the ideal configuration is compromised by including additional couplings, resulting in heat leaks. Model B. Here we distinguish between 'dominant' couplings, which reproduce the ideal setup, γ = γ C 12 = γ W 23 = γ H 13 , and secondary-weaker couplings, which lead to heat leakage and competition between cycles,γ = γ C,W 13 = γ C,H

12
, withγ ≪ γ. Results are presented in panel b. We find that the cooling function disappears for small E 2,1 , but it survives in the intermediate E 2,1 ∼ E 3,2 regime. We can rationalize our observations as follows (Appendix B): While heat directly leaks to the cold bath from both the hot and work reservoirs, we may extract energy from the cold bath using the different cycles. We note that F 3,1 cyc is always negative thus we are left with two cycles, F 2,1 cyc and F 3,2 cyc . The cooling condition for cycle F 2,1 cyc is In contrast, cooling via the cycle F 3,2 cyc is achieved if we satisfy The two cycles have conflicting requirements on the energy gap E 2,1 , thus cooling is not achievable for either small or large E 2,1 , when one of the cycles delivers heat into the cold bath. Furthermore, heat leaks and the noncooling F 3,1 cyc cycle suppress the cooling current for both small and large E 2,1 spacing.
Model C. In this model (panel c) the hot bath directly competes with the cold bath on the cooling performance. As we show in Appendix B, the cooling condition of the setup can be compactly written as with F 2,1 leak,H ∝ e −βCE2,1 − e −βH E2,1 , which is negative. Therefore, the extracted current is smaller than in Model A, and the cooling condition is more difficult to satisfy, resulting in a more limited window of operation. Specifically, beyond a certain value for E 2,1 , the cooling function is lost since the heat leak term takes over.
Model D. This setup (panel d) is similar to Model C, but with the work bath directly competing with the cold bath. A cooling condition, analogous to Eq. (19), can be written. However, given that β W ≪ β C , the cooling performance is significantly reduced, and it survives only for very small energy splitting, E 2,1 .
The secular approximation breaks down when two quasidegenerate levels are coupled to a third level, with the same bath responsible for transitions between the ground and excited states. In Fig. 2, the energy level E 2 is allowed to move in the full range between E 1 and E 3 , thus reaching degeneracy in the extreme limits. However, our analysis and conclusions do not pertain to these special limits. Specifically, in panel b the region of no-cooling at small E 2,1 is controlled byγ, and not by coherent effects. In panel d, our main observation is the suppression of the central cooling region (when E 2 is set about halfway in the E 3,1 gap). This suppression is controlled by a leakage process from the work bath.
In Fig. 3 we display examples based on Fig. 2 and highlight the different magnitudes of the cooling currents. Because T H is close to T C , the cooling behavior is moderately-negatively impacted in models B and C relative to model A. In contrast, heat leaks from the work to the cold reservoir (model D) dramatically lessen both the cooling window and the cooling current.  Fig. 2 for β H = 0.9.

IV. SUMMARY
We derived a cooling condition for a multilevel quantum absorption refrigerator described at the level of a Markovian quantum master equation. The obtained cooling condition is given in terms of the generator of the dynamics, and it provides analytical insight to the machine operation. It allows us to determine the effect of heat leaks and competing cycles in the setup, and feasibly identify the parameter space that achieves the desired function. We further obtained a closed-form expression for the currents (cooling or input power) in the model, and described the calculation of the current noise. We expect this work to assist in the quest to realize QARs; a realistic device evidently suffers from leakage and internal dissipation processes. While we exemplified our results on a three-level QAR, the formalism holds for a general N -level system, with the heat current combining cycles and leaks as in Eq. (15).
Steady-state coherences arise in V and Λ-type level structure when nonequilibrium baths are coupled to the energy degenerate transitions 18 . The impact of such coherences on the cooling current can be studied numerically with a non-secular MQME [18][19][20]23 . A recent study showed that under some assumptions, quantum coherent devices can be mapped onto a "classical emulator" that reproduces the same thermodynamic performance in the long time limit 52 , where our formalism could be applied. Achieving a compact, feasible cooling condition in cases that cannot be described by a classical-incoherent emulator is left for future work. We show that Eqs. (9) and (12) reproduce known results for the heat current and its noise in the nonequilibrium spin boson model (N = 2), displayed in Fig. (4). The Hamiltonian of the model iŝ . The counting field dependent MQME readsṗ(t, χ) = L(χ)p(t, χ) with the Liouvillian rate matrix 46 , In this equation, rate constants describing heat exchange with the cold bath are decorated by the counting fields.
The 'bare' rate constants are k µ 1→2 = ∞ −∞ dt e −iω0t B µ (t)B µ (0) , and the detailed balance condition enforces k µ 1→2 = e −βµω0 k µ 2→1 . To use equation (9) we calculate the derivative, We note the natural form of this matrix, allowing us to create it intuitively: It includes heat transfer rates to/from the cold bath with the appropriate sign convention. Using Eq. (9) we can now determine the current from the C bath, The matrix filled with the cofactor terms C ij is the adjugate matrix, adj[L(0)]. a 1 (0) is the coefficient of λ in the characteristic polynomial, Eq. (4). Simplifying, and marking only the relevant terms, while using the fact that a 1 (0) = Tr [L(0)], we get We need only solve for C 12 and C 21 from the adjugate matrix. Trivially, Using the harmonic bath model and bilinear system-bath coupling the rate constants are k µ 2→1 = Γ µ (ω 0 )[n µ (ω 0 ) + 1] and k µ 1→2 = Γ µ (ω 0 )n µ (ω 0 ), therefore Here, Γ µ (ω 0 ) = 2π k (u µ,k ) 2 δ(ω 0 − ω µ,k ) is the spectral density of the µ bath evaluated at ω 0 and n µ (ω 0 ) = (e βµω0 − 1) −1 is the Bose-Einstein distribution function. Equation (A4) agrees with previous studies [53][54][55] . However, the present evaluation, which is based on Eq. (A2) is almost trivial compared to the standard approach, which requires one to solve the equation of motion for the population in steady state, then substitute the population in the expression for the current.
Beyond the current, we follow the procedure for calculating the current noise using Eq. (12). In this case we need to examine two matrices ∂ 2 L(χ)/∂(iχ) 2 and ∂adj [L(χ)] /∂(iχ). The former is easy to calculate, As for the latter, again, we only need to study two elements, Altogether, we get which reduces to the known result 55 Here, Γ µ and n µ are evaluated at the frequency of the spin, ω 0 .
FIG. 4: Scheme of the nonequilibrium two-level system. By introducing a counting parameter on the cold bath, we calculate the heat current extracted from this bath to the system, which is negative in this two-terminal problem.

APPENDIX B: IDEAL AND NON-IDEAL THREE-LEVEL QAR
We describe here the calculation of the cooling window in the different 3-level models presented in the four panels of Fig. 2.

A. Ideal QAR: Model A
We exemplify our procedure to obtain a cooling condition [Eq. (10)] on the canonical three-level quantum absorption refrigerator with 'ideal' coupling operator. The Hamiltonian of the model readŝ whereB µ ij = k u µ,k ij (â † µ,k + a µ,k ) is an operator of the µ heat bath. The equations of motion for the populations are given byṗ(t) = L(0) p(t) with the population vector p(t) = (p 1 (t), p 2 (t), p 3 (t)) T and the Liouvillian The rate constants are The detailed balance relation dictates the rate constants of reversed processes, e.g., k C 2→1 = Γ C (E 2,1 ) [n C (E 2,1 ) + 1].
The cumulant generating function of the model can be derived by following a rigorous procedure 43,46 . Here, we employ a classical, intuitive derivation of the FCS; it agrees with the rigorous method under the weakcoupling, Markovian and secular approximations 46 . We begin by defining P t (j, nE 2,1 ) as the probability that by the (long) time t, the thermal energy nE 2,1 had been absorbed by the system (transferred from the cold bath), and the system occupies the state j = 1, 2, 3. Note that, for simplicity, we count energy exchanged with the cold bath only. Here, n is an integer since under the weak coupling approximation heat is transferred into and out of the system in discrete quanta, in resonance with system's level spacing. We readily write down an equation of motion for P t (j, nE 2,1 ) 43,55,56 , The equation is Fourier-transformed with the counting field χ, and we define the counting field dependent population Note: the sign convention in the Fourier transform dictates the sign convention for the heat current, with positive current corresponding to heat absorbed by the system. Using Eq. (B4), we write down the set of differential equations,ṗ with the rate matrix (Liouvillian) The characteristic function is given by Z(χ, t) = I|p(χ, t) with the identity vector I|. We now employ Eq. (10) to determine the cooling condition for the setup. We calculate the derivatives of the elements in L(χ), Again, we remark on the intuitive form of this matrix, pointing to the cooling processes. It can be constructed without prior knowledge of the FCS method. The cooling condition is The matrix filled with the cofactors C ij is the adjugate matrix of the Liouvillian, L(0). This cooling condition condenses into If we interpret the cofactors as effective population, C 21 ↔ p 1 and C 12 ↔ p 2 , we rationalize the cooling condition as (nonequilibrium) deviations from the detailed balance relation. We now evaluate the cofactors, and substitute them in Eq. (B9). Altogether we get This inequality can be simplified to arrive at the wellknown cooling condition The formalism provides the cooling current, and following similar steps, we get the heat absorbed from the work reservoir as well, thus the cooling coefficient of performance, η ≡ J C / J W = E 2,1 /E 3,2 . Combining this with Eq. (B12) we conclude that η ≤ η c with the Carnot cooling efficiency η c = (β H − β W )/(β C − β H ). Since one can tune parameters to reach the maximal cooling performance, we refer to this model as 'ideal'. This concludes the derivation of the cooling window for Model A presented in Fig. 2(a).

B. Non-ideal QAR: Model B
We analyze here Model B, as described and presented in Fig. 2(b). The Hamiltonian is given bŷ In this model, each transition in the system is coupled to the three reservoirs, albeit we distinguish between 'dominant' couplings as in the ideal QAR model, and secondary (weaker) couplings. Following the standard steps as performed for the ideal model, we arrive at an equation of motion for the counting field-dependent population, p(χ, t) = L(χ)p(χ, t) with the Liouvillian Recall that counting is performed only at the cold bath. The derivative is given as which generalizes (B8), encompassing additional coldbath mediated heat exchange processes. Using Eq. (9), we obtain a condition on cooling, The three contributions arise due to the coupling of the cold bath to the three transitions, |1 ↔ |2 , |1 ↔ |3 , and |2 ↔ |3 . We begin by analyzing the second term in Eq. (B16), and we write it down in terms of the rate constants, Since β C > β H > β W , this contribution is always negative hindering cooling through the 1 ↔ 3 transition. Next, we consider the first row, (C 21 e −βCE2,1 − C 12 ) in Eq. (B16) and write it explicitly in terms of the rate constants. We find that it includes heat leaks that are always negative, as well as one term that could be positivedepending on the energy structure, Similarly, the only potentially-positive combination in the last row of Eq. (B16), (C 32 e −βCE3,2 − C 23 ), is F 3,1 cyc is given by the first contribution in Eq. (B17), and it is negative.
We conclude that there are two possible cooling cycles in the three level model: (i) F 2,1 cyc : The cold bath excites the system allowing the transition |1 → |2 . The work bath pumps heat and induces the transition |2 → |3 . Heat is released to the H bath with the system relaxing to |1 . (ii) F 3,2 cyc : The cold bath excites the system, |2 → |3 . The work bath induces the transition |1 → |2 . Heat is released to the H bath with the system relaxing to |1 . These two circuits compete ("internal dissipation") and it is impossible to maximize the cooling performance simultaneously and achieve maximal (Carnot) cooling: Trying to optimize cycle (i) leads to suppression of the cooling current through cycle (ii), and vice versa 10,37 . Particularly, conditions (B22) are organized into the cooling conditions (17)- (18), illustrating the conflicting requirements that the different cycles impose on the device.
C. Non-ideal QAR: Models C and D We analyze here Models C and D, as described in Fig.  2(c)-(d). Recalling, in Model C (D) the hot (work) bath leaks energy directly to the cold bath; the cold bath couples at the |1 → |2 transition. The last term describes the 'leaky' (l) coupling: l = H (W ) in Model C (D). Following similar definitions as for the ideal model, we arrive at the equation of motion for the counting fielddependent populations,ṗ(χ, t) = L(χ)p(χ, t) with the Liouvillian The derivative matrix ∂L/∂(iχ)| χ=0 is the same as in the ideal case, resulting in a formally identical cooling condition, yet with different cofactors C 21 and C 12 , We plug these expressions back into the cooling condition (B25) and get e −βC E2,1−βW E3,2 − e −βH E3,1 The first square-bracket term corresponds to the cooling condition of the ideal model A. The second square bracket describes the heat flow from a hot bath (H or W ) to the cold bath. Since the second term is always negative, we conclude that (as expected) the cooling condition is more limited in Model C or D, compared to the ideal design, Model A.