Qubit absorption refrigerator at strong coupling

An autonomous absorption refrigerator transfers thermal energy from a cold (c) bath to a hot (h) bath without input power, by utilizing heat from an additional heat bath, a so-called work (w) reservoir. Classical, large-scale absorption refrigerators were realized in the 19th century [1], playing an important role in the development of the theory of irreversible thermodynamics. Proposals for quantum, nanoscale analogs of such machines aspire to establish the theory of thermodynamics from quantum principles [2–5].

Quantum thermodynamical machines differ from their classical counterparts in two central aspects. First, their performance relies on quantum phenomena such as the discreteness of the energy spectrum of the working medium and quantum statistics. Moreover, nontrivial quantum effects such as quantum coherence in the system, [6–11] or in the bath [12–15], quantum correlations [16, 17], non-locality, measurement [18, 19], and quantum driving and control [20, 21], may offer new principles for thermal machines. Beyond quantum resources, a second, fundamental aspect of nanoscale heat machines is that they may operate beyond the weak system-bath coupling limit [22–30]. Classical-macroscopic thermodynamics is a weak-coupling theory; the impact of the system-bath interface is small relative to the bulk behavior. In contrast, small systems can strongly couple to their surroundings, in the sense that the interaction energy between the system and the bath becomes comparable to frequencies of the isolated system.

The goal of the present paper is to demonstrate that strongly-coupled system-bath quantum machines can exhibit function that is un-attainable under the weak coupling assumption. We do so by analyzing additive and non-additive system-bath interaction models, i.e. where the generator may or may not be additively decomposed into individual generators from the connected reservoirs. In the additive case, the system (working medium) separately-independently exchanges energy with the hot, cold and work reservoirs. In the non-additive model, the reservoirs inseparably interact with the system, thus acting in a concerted-cooperative manner.

Specifically, we show that a two-level system (TLS) cannot operate as an autonomous quantum absorption refrigerator (QAR) under the weak system-bath coupling approximation with additive dissipators. However, the same system does function as a QAR once it is allowed to couple to its surrounding reservoirs in a non-additive manner—representing strong coupling. Moreover, the qubit QAR can be optimized to perform at the maximal Carnot efficiency, and its performance is compatible with previous designs using three-level or three-qubit models, which were constrained to operate under Lindblad dynamics with additive dissipators [31–34]. The smallest possible QAR described here relies on quantum principles and strong system-bath coupling effects. These unique aspects are inherent to nanoscale devices.

This work is organized as follows: we first introduce our model, showing that a QAR mode is impossible in the additive case (section 2.1), and afterwards present the non-additive model (section 2.2), for which we present the basic definitions of energy currents for two and three reservoirs. Next, we present analytical results for the non-additive model, first on the cooling window and efficiency for specific spectral densities in section 3.1 and then on the cooling efficiency at maximum power in section 3.2. We confirm these results by numerical simulations in section 4. Finally, we show explicitly how such a non-additive dissipator may arise in the strong-coupling limit in section 5.

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}\gt {T}_{h}\gt {T}_{c}$, see figure 1(a). 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. [4, 32, 33]. 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.

Zoom In Zoom Out Reset image size

In this paper we focus on a QAR made of a TLS coupled to three independent thermal baths. When the baths couple to the qubit in an additive manner (figure 1(b)) 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 (figure 1(c)), we are able to obtain a cooling condition, as well as derive bounds for the maximal efficiency of the QAR and its maximal power efficiency (MPE).

2.1. Un-attainability of cooling for an additive dissipation model

The additive model comprises a TLS (spin, qubit) and three independent thermal reservoirs $\nu =c,h,w$, ${T}_{w}\gt {T}_{h}\gt {T}_{c}$, ${\beta }_{\nu }=1/{T}_{\nu }$ with k_B = 1. The generic Hamiltonian is written as

$$\begin{eqnarray}&&\hat{H}=\displaystyle \frac{{\omega }_{0}}{2}{\hat{\sigma }}_{z}+\displaystyle \sum _{\nu }{\hat{H}}_{B,\nu }+\displaystyle \frac{{\hat{\sigma }}_{x}}{2}\otimes ({\hat{A}}_{c}+{\hat{A}}_{h}+{\hat{A}}_{w}).\end{eqnarray} \tag{ 1 }$$

Here, $\hat{\sigma }$ are the Pauli spin matrices and ${\hat{H}}_{B,\nu }$ is the Hamiltonian of the νth reservoir. It includes, for example, a collection of harmonic oscillators of frequencies ${\omega }_{j,\nu }$, ${\hat{H}}_{B,\nu }={\sum }_{j}{\omega }_{j,\nu }{\hat{b}}_{j,\nu }^{\dagger }{\hat{b}}_{j,\nu }$ with ${b}^{\dagger }$ ($\hat{b}$) as bosonic creation (annihilation) operators. The bath operator ${\hat{A}}_{\nu }$ is assumed to be hermitian. It couples the $\nu th$ bath to the spin, where e.g. ${\hat{A}}_{\nu }={\sum }_{j}{\lambda }_{j,\nu }({\hat{b}}_{j,\nu }^{\dagger }+{\hat{b}}_{j,\nu })$ with coupling strength ${\lambda }_{j,\nu }$.

Assuming a factorized-product initial state, $\langle {\hat{A}}_{\nu }\rangle =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 [35], with the dissipator given by the sum of three (c, h, w) independent contributions. This standard Born–Markov scheme results in the stationary populations of the excited and ground state, respectively,

$$\begin{eqnarray}&&{p}_{e}=\displaystyle \frac{{k}_{u}}{{k}_{d}+{k}_{u}},\,\,{p}_{g}=\displaystyle \frac{{k}_{d}}{{k}_{d}+{k}_{u}},\end{eqnarray} \tag{ 2 }$$

with ${k}_{d,u}={\sum }_{\nu }{k}_{d,u}^{(\nu )}$. The decay (d) and excitation (u) rate constants ${k}_{d,u}^{(\nu )}$, induced by the νth bath, depend on the details of the model. The detailed balance relation dictates local thermal equilibrium, ${k}_{u}^{(\nu )}/{k}_{d}^{(\nu )}={{\rm{e}}}^{-{\beta }_{\nu }{\omega }_{0}}$. The energy current, defined positive when flowing towards the qubit, can be similarly derived from the Born–Markov approximation, and it is given by ${J}_{c}=-{\omega }_{0}({k}_{d}^{(c)}{p}_{e}-{k}_{u}^{(c)}{p}_{g})$ [36–39]. Substituting the steady state population (2) we obtain

$$\begin{eqnarray}&&{J}_{c}=-\displaystyle \frac{{\omega }_{0}}{{k}_{d}+{k}_{u}}\left[{k}_{d}^{(c)}{k}_{d}^{(h)}\left(\displaystyle \frac{{k}_{u}^{(h)}}{{k}_{d}^{(h)}}-\displaystyle \frac{{k}_{u}^{(c)}}{{k}_{d}^{(c)}}\right)+{k}_{d}^{(c)}{k}_{d}^{(w)}\left(\displaystyle \frac{{k}_{u}^{(w)}}{{k}_{d}^{(w)}}-\displaystyle \frac{{k}_{u}^{(c)}}{{k}_{d}^{(c)}}\right)\right].\end{eqnarray} \tag{ 3 }$$

Using the detailed-balance relation and the fact that $({{\rm{e}}}^{-{\beta }_{h,w}{\omega }_{0}}-{{\rm{e}}}^{-{\beta }_{c}{\omega }_{0}})\gt 0$, we conclude that ${J}_{c}\lt 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. [40–47]. These type of driven machines are beyond the scope of our work.

Quite generally, the dissipator of the weak coupling, Markovian master equation derived for the model Hamiltonian (1) is additive in the different reservoirs. The breakdown of additivity, e.g. beyond the weak coupling limit, is a topic of recent interest [30, 48, 49].

2.2. 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 strong-coupling characteristic, a non-additive system-bath interaction operator,

$$\begin{eqnarray}&&\hat{H}=\displaystyle \frac{{\omega }_{0}}{2}{\hat{\sigma }}_{z}+\displaystyle \sum _{\nu }{\hat{H}}_{B,\nu }+\gamma \displaystyle \frac{{\hat{\sigma }}_{x}}{2}\otimes ({\hat{B}}_{c}\otimes {\hat{B}}_{h}\otimes {\hat{B}}_{w}).\end{eqnarray} \tag{ 4 }$$

Here, ${\hat{B}}_{\nu }$ 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 [36–38], see section 5. 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. [4, 50, 51] and realized experimentally in a recent study [52]. In [4, 50–52], 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 equation (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.

2.2.1. Two-bath model

Equations of motion for the spin polarization, as well as the energy current, were derived in [38, 53] for the model (4) with two baths (hot and cold). The derivation relies on the assumption $\langle {\hat{B}}_{\nu }\rangle =0$, which could be satisfied exactly or under conditions such as strong coupling or high temperature [38]. Further, this assumption can be relaxed by re-defining the model Hamiltonian to add and subtract the thermal average of the interaction Hamiltonian, re-diagonalizing then the system's Hamiltonian and proceeding with the perturbative treatment [54, 55]. The population dynamics [36–38] satisfies

$$\begin{eqnarray}&&{\dot{p}}_{e}=-M({\omega }_{0}){p}_{e}(t)+M(-{\omega }_{0}){p}_{g}(t),\end{eqnarray} \tag{ 5 }$$

with the normalization condition ${p}_{g}(t)+{p}_{e}(t)=1$. The rate constants are

$$\begin{eqnarray}M({\omega }_{0}) & = & {\left(\displaystyle \frac{\gamma }{2}\right)}^{2}{\int }_{-\infty }^{\infty }{{\rm{e}}}^{{\rm{i}}{\omega }_{0}t}{M}_{h}(t){M}_{c}(t){\rm{d}}{t}\\ & = & \displaystyle \frac{1}{2\pi }{\left(\displaystyle \frac{\gamma }{2}\right)}^{2}{\displaystyle \int }_{-\infty }^{\infty }{M}_{h}({\omega }_{0}-\omega ){M}_{c}(\omega ){\rm{d}}\omega .\end{eqnarray} \tag{ 6 }$$

Here, ${M}_{\nu }(t)=\langle {\hat{B}}_{\nu }(t){\hat{B}}_{\nu }(0)\rangle$ is the two-time correlation function with ${B}_{\nu }(t)={{\rm{e}}}^{{{\rm{i}}{H}}_{B,\nu }t}{B}_{\nu }{{\rm{e}}}^{-{{\rm{i}}{H}}_{B,\nu }t}$ denoting the interaction picture and the average performed with respect to the canonical (initial) state of the $\nu th$ thermal bath. In Fourier space we introduce ${M}_{\nu }(\omega )={\int }_{-\infty }^{\infty }{{\rm{e}}}^{{\rm{i}}\omega t}{M}_{\nu }(t){\rm{d}}{t}$. 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 (${\hslash }=1$). Formally similar to P(E) theory [56], the detailed-balance condition is satisfied for the individual components, $\tfrac{{M}_{\nu }(\omega )}{{M}_{\nu }(-\omega )}={{\rm{e}}}^{{\beta }_{\nu }\omega }$, but we do not have such a relation for the convoluted rate constant $M(\pm {\omega }_{0})$.

Since the reservoirs act in a cooperative manner, one cannot identify a dissipation process of the system to an individual bath. Certainly, the energy current should reflect the coordinated multi-bath nature of energy transfer processes in the system. A formal expression for the thermal energy current in this non-additive model can be derived based on a full counting statistics approach. In [38], we derived the cumulant generating function of the model (4) by energy resolving the Markovian quantum master equation (5). Alternatively, in [53–55] the cumulant generating function was derived based on a two-time measurement protocol, by devising a counting field-dependent reduced density operator. The result of either procedures is a rather intuitive expression for the thermal energy current entering the system from the cold reservoir [38, 53],

$$\begin{eqnarray}&&{J}_{c}=-{p}_{{e}}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}\omega \omega {M}_{c}(\omega ){M}_{h}({\omega }_{0}-\omega )+{p}_{g}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}\omega \omega {M}_{c}(-\omega ){M}_{h}(\omega -{\omega }_{0}).\end{eqnarray} \tag{ 7 }$$

Here, we have absorbed a factor ${C}_{2}=\tfrac{1}{2\pi }{\left(\tfrac{\gamma }{2}\right)}^{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 ${\omega }_{0}-\omega$ absorbed or released by the hot bath. A similar reasoning applies to the second term in equation (7), which describes the excitation of the qubit. For simplicity, in what follows we eliminate the spin gap and take ${\omega }_{0}=0$. We then immediately conclude that ${p}_{e}={p}_{g}=1/2$, thus the heat current (7) simplifies to

$$\begin{eqnarray}{J}_{c} & = & -\displaystyle \frac{1}{2}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}\omega \omega {M}_{c}(\omega ){M}_{h}(-\omega )+\displaystyle \frac{1}{2}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}\omega \omega {M}_{c}(-\omega ){M}_{h}(\omega )\\ & = & \displaystyle \frac{1}{2}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}\omega \omega {M}_{c}(\omega ){M}_{h}(\omega )[{{\rm{e}}}^{-{\beta }_{c}\omega }-{{\rm{e}}}^{-{\beta }_{h}\omega }]\\ & = & {\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}\omega \omega {M}_{c}(-\omega ){M}_{h}(\omega ).\end{eqnarray} \tag{ 8 }$$

These three expressions offer three different useful representations of the heat current. We obtain the second representation from the first line by using the detailed balance relation for ${M}_{c,h}(-\omega )$. The third expression is obtained from the first one after changing the integration variable in the first integral. The cumulant generating function of the model satisfies the steady state heat exchange fluctuation theorem [38, 53–55]. Therefore, our framework upholds the second law of thermodynamics, which is a direct consequence of the exchange fluctuation relation. From this, it is easy to prove that if ${\beta }_{c}\gt {\beta }_{h}$, ${J}_{c}\lt 0$, meaning that the heat current flows towards the cold bath.

2.2.2. 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 [53], in a complete analogy to the two-bath case described in section 2.2.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 equation (5), with the FBCF now however given by

$$\begin{eqnarray}{M}({{\omega }}_{0}) & = & {\left(\displaystyle \frac{{\gamma }}{2}\right)}^{2}{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}{{\rm{e}}}^{{\rm{i}}{{\omega }}_{0}{t}}{{M}}_{{c}}({t}){{M}}_{{h}}({t}){{M}}_{{w}}({t}){\rm{d}}{t}\\ & = & \displaystyle \frac{1}{{(2{\pi })}^{2}}{\left(\displaystyle \frac{{\gamma }}{2}\right)}^{2}{\displaystyle \int }_{-{\rm{\infty }}}^{{\rm{\infty }}}{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}{{M}}_{{c}}({{\omega }}_{1}){{M}}_{{h}}({{\omega }}_{2}){{M}}_{{w}}({{\omega }}_{0}-{{\omega }}_{1}-{{\omega }}_{2}){\rm{d}}{{\omega }}_{1}{\rm{d}}{{\omega }}_{2}.\end{eqnarray} \tag{ 9 }$$

The energy current, from the cold bath towards the qubit, is given by

$$\begin{eqnarray}{J}_{c} & = & -{p}_{e}{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{\omega }_{1}{M}_{c}({\omega }_{1}){M}_{h}({\omega }_{2}){M}_{w}({\omega }_{0}-{\omega }_{1}-{\omega }_{2})\\ & & +{p}_{g}{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{\omega }_{1}{M}_{c}(-{\omega }_{1}){M}_{h}(-{\omega }_{2}){M}_{w}({\omega }_{1}+{\omega }_{2}-{\omega }_{0}).\end{eqnarray} \tag{ 10 }$$

Again, we have absorbed a prefactor ${C}_{3}=\tfrac{1}{{(2\pi )}^{2}}{\left(\tfrac{\gamma }{2}\right)}^{2}$ into the definition of the current. This intuitive expression, which can be also suggested phenomenologically [57], describes coordinated three-bath energy exchange processes, with an overall conservation of energy. An amount of energy ${\omega }_{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. One can further simplify the general model Hamiltonian and consider the degenerate model ${\omega }_{0}=0$. In this case, we obtain from equation (10),

$$\begin{eqnarray}{J}_{c} & = & -\displaystyle \frac{1}{2}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\omega }_{1}{M}_{c}({\omega }_{1}){\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){M}_{w}({\omega }_{1}+{\omega }_{2})[{{\rm{e}}}^{-{\beta }_{w}({\omega }_{1}+{\omega }_{2})}-{{\rm{e}}}^{-{\beta }_{c}{\omega }_{1}}{{\rm{e}}}^{-{\beta }_{h}{\omega }_{2}}]\\ & = & -{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{\omega }_{1}{M}_{c}({\omega }_{1}){M}_{h}({\omega }_{2}){M}_{w}(-{\omega }_{1}-{\omega }_{2}).\end{eqnarray} \tag{ 11 }$$

The first representation was derived from equation (10) by using the bath-specific detailed balance relations. The second line was obtained from equation (10) by changing integration variables in the second integral.

In what follows, we first analyze the degenerate model with ${\omega }_{0}=0$. In appendix A, we extend our study to cover the non-degenerate model, and show that it can be optimally operated, in an equivalent way to the ${\omega }_{0}=0$ case. Simulations at finite ${\omega }_{0}$ exhibit rich features; the energy gap ${\omega }_{0}$ offers a non-monotonic control over the cooling current.

While we write here an expression for J_c only, a full-counting statistics approach [53] readily hands over analogous expressions for J_w and J_h. Particularly, J_h is obtained from equation (10) by interchanging the 'c' and 'h' indices, and J_w then follows from energy conservation.

3.1. 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 equation (11), i.e., find out whether we can engineer the system and the baths to achieve refrigeration, ${J}_{c}\geqslant 0$.

The FBCF ${M}_{\nu }(\omega )$ is related to the spectral density of the ν thermal bath, see section 5. Let us assume that the baths are engineered such that these functions are characterized by the frequency ${\theta }_{\nu }\gt 0$, satisfying a resonant assumption

$$\begin{eqnarray}&&{\theta }_{c}+{\theta }_{w}={\theta }_{h}.\end{eqnarray} \tag{ 12 }$$

We now analyze the performance of our system as a QAR in an idealized limit, then under less restricting 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.

3.1.1. Ideal design

In an optimal design, the FBCF ${M}_{\nu }(\omega )$ is restricted to be nonzero within a narrow spectral window. In the most extreme case, we filter for each reservoir all frequency components besides the central mode⁵

$$\begin{eqnarray}&&{M}_{\nu }(\omega )={\epsilon }_{\nu }[\delta (\omega -{\theta }_{\nu })+\delta (\omega +{\theta }_{\nu }){{\rm{e}}}^{+{\beta }_{\nu }\omega }],\end{eqnarray} \tag{ 13 }$$

where ${\epsilon }_{\nu }\gt 0$ is a dimensionless parameter. The Dirac delta function here represents a sharp, normalized function. For example, it could be viewed as a limit of a Gaussian function with a small width parameter. We evaluate the integrals in equation (11), and find that there are only two contributions to the cooling current, when ${\omega }_{1}\leqslant 0$ and ${\omega }_{2}\geqslant 0$, and the other way around. We then gather the cooling condition, ${J}_{c}\geqslant 0$,

$$\begin{eqnarray}&&{\theta }_{c}({{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}-{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}})\geqslant 0.\end{eqnarray} \tag{ 14 }$$

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,

$$\begin{eqnarray}&&\left(\displaystyle \frac{{T}_{w}-{T}_{h}}{{T}_{w}-{T}_{c}}\right)\displaystyle \frac{{T}_{c}}{{T}_{h}}\geqslant \displaystyle \frac{{\theta }_{c}}{{\theta }_{h}},\end{eqnarray} \tag{ 15 }$$

which precisely corresponds to the cooling condition as obtained for a three-level or three-qubit QAR—analyzed with a Markovian master equation with additive dissipators [4, 32, 34]. While our cooling condition is given in terms of baths characteristic frequencies, in weak-coupling multi-level designs the analogous cooling condition is expressed in terms of the subsystem energies.

This result, along with an expression for the efficiency of the refrigerator (equation (27) below), are significant conclusions of our work: to realize a weak coupling three-qubit or three-level QAR, a resonant condition on the subsystem energies is enforced. As well, one needs to engineer ('filter') the interaction within the subsystem, to permit only particular transitions. Our analysis here demonstrated that an equivalent optimal performance for a QAR is achieved in the collective interaction model of the baths. In this inseparable case, the spectral functions of the reservoirs are filtered. Moreover, these spectral functions are enforced to satisfy a resonant condition so as to only admit processes in which an amount of energy of ${\theta }_{c}$ + ${\theta }_{w}$, from/to the cold and work baths is exchanged with quanta of ${\theta }_{h}$ to/from the hot reservoir.

3.1.2. Non-ideal design

Let us now consider a potentially more feasible—and less optimal model. We engineer the FBCFs as follows ($\omega \gt 0$),

$$\begin{eqnarray}{M}_{c}(\omega ) & = & {\epsilon }_{c}\delta (\omega -{\theta }_{c}),\\ {M}_{h}(\omega ) & = & {\epsilon }_{h}[H(\omega -{\theta }_{h}+\delta )-H(\omega -{\theta }_{h}-\delta )],\\ {M}_{w}(\omega ) & = & {\epsilon }_{w}[H(\omega -{\theta }_{w}+\delta )-H(\omega -{\theta }_{w}-\delta )].\end{eqnarray} \tag{ 16 }$$

For negative frequencies, the detailed balance relation multiplies each function by a thermal factor. Here, H(x) is the Heaviside step function, $\delta (x)$ is the Dirac Delta function, δ (without an argument) is the width parameter. The prefactors ${\epsilon }_{\nu }$ determine the magnitude of the FBCF⁶ . As we show below, the parameters ${\theta }_{\nu }$ and δ critically determine the performance of the QAR. In contrast, ${\epsilon }_{\nu }$, which characterizes the 'strength' of the reservoirs' spectral functions, affects neither the cooling window nor the efficiency. We assume the resonant assumption (12) to hold, and that the width parameter δ is small, $\delta \ll {\theta }_{\nu }$.

We calculate the cooling current based on equation (11) by breaking it into four contributions. Let us first inspect the ${\omega }_{1}\geqslant 0$ and ${\omega }_{2}\geqslant 0$ term,

$$\begin{eqnarray}{{J}}_{{c}}^{(+,+)} & = & -{\displaystyle \int }_{0}^{{\rm{\infty }}}{\rm{d}}{{\omega }}_{1}{{\omega }}_{1}{{M}}_{{c}}({{\omega }}_{1}){\displaystyle \int }_{0}^{{\rm{\infty }}}{\rm{d}}{{\omega }}_{2}{{M}}_{{h}}({{\omega }}_{2}){{M}}_{{w}}({{\omega }}_{1}+{{\omega }}_{2}){{\rm{e}}}^{-{{\beta }}_{{w}}({{\omega }}_{1}+{{\omega }}_{2})}\\ & = & -{{\theta }}_{{c}}{{\epsilon }}_{{c}}{{\varepsilon }}_{{h}}{\displaystyle \int }_{{{\theta }}_{{h}}-{\delta }}^{{{\theta }}_{{h}}+{\delta }}{\rm{d}}{{\omega }}_{2}{{M}}_{{w}}({{\theta }}_{{c}}+{{\omega }}_{2}){{\rm{e}}}^{-{{\beta }}_{{w}}({{\theta }}_{{c}}+{{\omega }}_{2})}.\end{eqnarray} \tag{ 17 }$$

Since ${\theta }_{c}+{\theta }_{h}-\delta \gt {\theta }_{w}+\delta$, the integral collapses to null. A similar argument brings a zero contribution from the negative branch, ${J}_{c}^{(-,-)}=0$, which includes ${\omega }_{1}\leqslant 0$ and ${\omega }_{2}\leqslant 0$. We proceed and evaluate the contribution to the cooling current from ${\omega }_{1}\leqslant 0$ but ${\omega }_{2}\geqslant 0$. Recall that ${\theta }_{h}-\delta \gt {\theta }_{c}$,

$$\begin{eqnarray}{J}_{c}^{(-,+)} & = & {\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{1}{\omega }_{1}{M}_{c}({\omega }_{1}){{\rm{e}}}^{-{\beta }_{c}{\omega }_{1}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){M}_{w}(-{\omega }_{1}+{\omega }_{2}){{\rm{e}}}^{-{\beta }_{w}(-{\omega }_{1}+{\omega }_{2})}\\ & = & {\epsilon }_{c}{\theta }_{c}{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){M}_{w}({\omega }_{2}-{\theta }_{c}){{\rm{e}}}^{-{\beta }_{w}({\omega }_{2}-{\theta }_{c})}\\ & = & {\epsilon }_{c}{\epsilon }_{h}{\theta }_{c}{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{\displaystyle \int }_{{\theta }_{h}-\delta }^{{\theta }_{h}+\delta }{\rm{d}}{\omega }_{2}{M}_{w}({\omega }_{2}-{\theta }_{c}){{\rm{e}}}^{-{\beta }_{w}({\omega }_{2}-{\theta }_{c})}\\ & = & \displaystyle \frac{2{\epsilon }_{c}{\epsilon }_{h}{\epsilon }_{w}}{{\beta }_{w}}{\theta }_{c}{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}\sinh ({\beta }_{w}\delta ).\end{eqnarray} \tag{ 18 }$$

We had utilized the resonant assumption (12) in the last line. A similar analysis gives

$$\begin{eqnarray}&&{J}_{c}^{(+,-)}=-\displaystyle \frac{2{\epsilon }_{c}{\epsilon }_{h}{\epsilon }_{w}}{{\beta }_{h}}{\theta }_{c}{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}\sinh ({\beta }_{h}\delta ).\end{eqnarray} \tag{ 19 }$$

Putting together equations (18) and (19), we organize the cooling condition, ${J}_{c}\geqslant 0$, as

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\beta }_{w}}{\theta }_{c}{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}\sinh ({\beta }_{w}\delta )-\displaystyle \frac{1}{{\beta }_{h}}{\theta }_{c}{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}\sinh ({\beta }_{h}\delta )\geqslant 0.\end{eqnarray} \tag{ 20 }$$

We Taylor-expand this result to the first nontrivial order in δ, which is the third order, and get

$$\begin{eqnarray}&&\left(2\delta +\displaystyle \frac{1}{3}{\beta }_{w}^{2}{\delta }^{3}\right){{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}-\left(2\delta +\displaystyle \frac{1}{3}{\beta }_{h}^{2}{\delta }^{3}\right){{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}\geqslant 0.\end{eqnarray} \tag{ 21 }$$

Using ${\theta }_{w}={\theta }_{h}-{\theta }_{c}$, we find that

$$\begin{eqnarray}&&({\beta }_{h}-{\beta }_{w}){\theta }_{h}-({\beta }_{c}-{\beta }_{w}){\theta }_{c}\geqslant \mathrm{ln}\displaystyle \frac{2\delta +\tfrac{1}{3}{\beta }_{h}^{2}{\delta }^{3}}{2\delta +\tfrac{1}{3}{\beta }_{w}^{2}{\delta }^{3}}.\end{eqnarray} \tag{ 22 }$$

Re-arranging this expression, we get a cooling condition

$$\begin{eqnarray}&&\left(\displaystyle \frac{{T}_{w}-{T}_{h}}{{T}_{w}-{T}_{c}}\right)\displaystyle \frac{{T}_{c}}{{T}_{h}}\geqslant \displaystyle \frac{{\theta }_{c}}{{\theta }_{h}}-\displaystyle \frac{{T}_{c}{T}_{w}}{{\theta }_{h}({T}_{w}-{T}_{c})}\mathrm{ln}\left[\displaystyle \frac{2+\tfrac{{{\beta }_{w}}^{2}{\delta }^{2}}{3}}{2+\tfrac{{{\beta }_{h}}^{2}{\delta }^{2}}{3}}\right].\end{eqnarray} \tag{ 23 }$$

In the limit $\delta \to 0$, we retrieve equation (14), which corresponds to the three-level QAR analyzed with a Markovian master equation with additive dissipators [4, 32]. It is important to recognize that the δ dependence in equation (23) is non-universal, and it depends on our choice (16). For example, using a different setting, with ${M}_{c}(\omega )$ and ${M}_{h}(\omega )$ as step functions of width $2\delta$, but ${M}_{w}(\omega )$ a Dirac delta function, we derive an alternative cooling condition,

$$\begin{eqnarray}&&\left(\displaystyle \frac{{T}_{w}-{T}_{h}}{{T}_{w}-{T}_{c}}\right)\displaystyle \frac{{T}_{c}}{{T}_{h}}\geqslant \displaystyle \frac{{\theta }_{c}}{{\theta }_{h}}-\displaystyle \frac{{T}_{c}{T}_{w}}{{\theta }_{h}({T}_{w}-{T}_{c})}\mathrm{ln}\left[\displaystyle \frac{2{\theta }_{c}+\tfrac{1}{3}{\theta }_{c}{{\beta }_{c}}^{2}{\delta }^{2}-\tfrac{2}{3}{\beta }_{c}{\delta }^{2}}{2{\theta }_{c}+\tfrac{1}{3}{\theta }_{c}{{\beta }_{h}}^{2}{\delta }^{2}-\tfrac{2}{3}{\beta }_{h}{\delta }^{2}}\right].\end{eqnarray} \tag{ 24 }$$

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, equations (23) and (24), are central results of our work. The cooling window depends on δ in a non-universal manner. In the optimal limit, $\delta \to 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, $\eta \equiv \tfrac{{J}_{c}}{{J}_{w}}$. 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 equation (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

$$\begin{eqnarray}{J}_{w} & = & \displaystyle \frac{{\epsilon }_{c}{\epsilon }_{h}{\epsilon }_{w}}{{\beta }_{h}}{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}\left[\left({\theta }_{w}+\delta +\displaystyle \frac{1}{{\beta }_{h}}\right){{\rm{e}}}^{-{\beta }_{h}\delta }-\left({\theta }_{w}-\delta +\displaystyle \frac{1}{{\beta }_{h}}\right){{\rm{e}}}^{{\beta }_{h}\delta }\right]\\ & & -\displaystyle \frac{{\epsilon }_{c}{\epsilon }_{h}{\epsilon }_{w}}{{\beta }_{w}}{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}\left[\left({\theta }_{w}+\delta +\displaystyle \frac{1}{{\beta }_{w}}\right){{\rm{e}}}^{-{\beta }_{w}\delta }-\left({\theta }_{w}-\delta +\displaystyle \frac{1}{{\beta }_{w}}\right){{\rm{e}}}^{{\beta }_{w}\delta }\right],\\ {J}_{c} & = & {\theta }_{c}{\epsilon }_{c}{\epsilon }_{h}{\epsilon }_{w}\left\{{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}[{{\rm{e}}}^{-{\beta }_{w}\delta }-{{\rm{e}}}^{{\beta }_{w}\delta }]\displaystyle \frac{1}{-{\beta }_{w}}+\displaystyle \frac{1}{{\beta }_{h}}{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}[{{\rm{e}}}^{-{\beta }_{h}\delta }-{{\rm{e}}}^{{\beta }_{h}\delta }]\right\}.\end{eqnarray} \tag{ 25 }$$

We expand the currents to third order in δ and accomplish the cooling COP

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{\theta }_{c}{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}[2\delta +\tfrac{1}{3}{{\beta }_{w}}^{2}{\delta }^{3}]-{\theta }_{c}{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}\left[2\delta +\tfrac{1}{3}{{\beta }_{h}}^{2}{\delta }^{3}\right]}{{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}\left[2{\theta }_{w}\delta -\tfrac{2}{3}{\beta }_{w}{\delta }^{3}+\tfrac{1}{3}{\theta }_{w}{{\beta }_{w}}^{2}{\delta }^{3}\right]-{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}\left[2{\theta }_{w}\delta -\tfrac{2}{3}{\beta }_{h}{\delta }^{3}+\tfrac{1}{3}{\theta }_{w}{{\beta }_{h}}^{2}{\delta }^{3}\right]}.\end{eqnarray} \tag{ 26 }$$

The width parameter δ could affect the efficiency in a non-monotonic way. However, to the lowest (linear) order in δ we gather

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{\theta }_{c}}{{\theta }_{w}}.\end{eqnarray} \tag{ 27 }$$

This result agrees with the efficiency derived for the three-level heat pump in an early, seminal work by Geusic et al [58]. The cooling condition (15) can be used to derive a bound on efficiency. In the ideal limit, the cooling window is defined from

$$\begin{eqnarray}&&\left(\displaystyle \frac{{T}_{w}-{T}_{h}}{{T}_{w}-{T}_{c}}\right)\displaystyle \frac{{T}_{c}}{{T}_{h}}\geqslant \displaystyle \frac{{\theta }_{c}}{{\theta }_{h}}=\displaystyle \frac{{\theta }_{c}}{{\theta }_{c}+{\theta }_{w}}.\end{eqnarray} \tag{ 28 }$$

We inverse this expression and reach

$$\begin{eqnarray}&&\left(\displaystyle \frac{{T}_{w}-{T}_{c}}{{T}_{w}-{T}_{h}}\right)\displaystyle \frac{{T}_{h}}{{T}_{c}}\leqslant 1+\displaystyle \frac{{\theta }_{w}}{{\theta }_{c}},\end{eqnarray} \tag{ 29 }$$

which can be also expressed as

$$\begin{eqnarray}&&\left(\displaystyle \frac{{T}_{h}-{T}_{c}}{{T}_{w}-{T}_{h}}\right)\displaystyle \frac{{T}_{w}}{{T}_{c}}\leqslant \displaystyle \frac{{\theta }_{w}}{{\theta }_{c}},\end{eqnarray} \tag{ 30 }$$

finally receiving $\tfrac{{\theta }_{c}}{{\theta }_{w}}\leqslant \left(\tfrac{{T}_{w}-{T}_{h}}{{T}_{h}-{T}_{c}}\right)\tfrac{{T}_{c}}{{T}_{w}}$. Comparing this expression to equation (27), we gain a bound on the cooling COP (in the $\delta \to 0$ limit),

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{\theta }_{c}}{{\theta }_{w}}\leqslant \left(\displaystyle \frac{{T}_{w}-{T}_{h}}{{T}_{h}-{T}_{c}}\right)\displaystyle \frac{{T}_{c}}{{T}_{w}}\equiv {\eta }_{c},\end{eqnarray} \tag{ 31 }$$

which is nothing but the Carnot bound. Specifically, when ${T}_{w}\gg {T}_{c},{T}_{h}$, we find that $\tfrac{{\theta }_{c}}{{\theta }_{w}}\leqslant {T}_{c}/({T}_{h}-{T}_{c})$, which is the Carnot bound for cooling machines. We emphasize that equation (31) was developed in previous studies, see e.g. [32], yet restricted to quantum systems that evolve with additive dissipators. The COP of non-ideal models with a finite width ${M}_{\nu }(\omega )$, can be also calculated and simulated, as we do in section 4. It is important to note that the COP of our model is always limited by the Carnot bound (31) as was recently proven in [14] based on the heat exchange fluctuation theorem.

The analysis presented in this subsection can be extended to describe a non-degenerate (biased) qubit ${\omega }_{0}\ne 0$ as we show in appendix A. We find that for small ${\omega }_{0}$ we can adjust the functions ${M}_{\nu }(\omega )$ 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 soft-cutoff functions ${M}_{\nu }(\omega )$, rather than the hard-cutoff model (16). In appendix B we calculate the cooling current (11) using specific, continuous, unimodal Gaussian ${M}_{\nu }(\omega )$ functions that are physically motivated from the eminent spin-boson model. In this model, the width of each FBCF is determined by the temperature of the associated bath. We prove that this model cannot act as a QAR. We further suggest a more complex, bimodal Gaussian function for the FBCF, in which the width parameter is determined independently of the temperature. We find that a QAR performance can be achieved in this case—with soft cutoffs for the spectral functions—for a fair range of parameters. For example, we realize a QAR when the function ${M}_{w}(\omega )$ quickly decays at high frequencies, beyond ${\theta }_{h}+{\theta }_{c}$. This requirement limits processes in which the work reservoir supplies energy to both the hot and cold baths.

3.2. 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 MPE ${\eta }^{* }$, meaning, that we calculate the cooling COP when maximizing the cooling current. An intriguing question, which recently received much attention, is whether ${\eta }^{* }$ can approach ${\eta }_{c}$— albeit at finite cooling power. For heat engines, bounds for the MPE were extensively investigated [59–65]. The MPE of refrigerators was examined in e.g. [32], demonstrating that it is constrained by the spectral properties of the thermal reservoirs. Nevertheless, [32] was concerned with additive (weakly-coupled) system-bath models.

To derive an expression for the efficiency at maximal cooling power, we go back to equation (25), and analyze it to the lowest order in δ,

$$\begin{eqnarray}{J}_{c} & = & 2\delta \times {\epsilon }_{c}{\epsilon }_{h}{\epsilon }_{w}{\theta }_{c}({{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}-{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}})\\ & = & 2\delta \times {\epsilon }_{c}{\epsilon }_{h}{\epsilon }_{w}{\theta }_{c}({{\rm{e}}}^{({\beta }_{w}-{\beta }_{c}){\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{h}}-{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}).\end{eqnarray} \tag{ 32 }$$

This function is linear in ${\theta }_{c}$ for small values, but it decays exponentially with ${\theta }_{c}$ beyond that. We thus search for ${\theta }_{c}$ which maximizes ${J}_{c}({\theta }_{c})$ by solving $\tfrac{\partial {J}_{c}}{\partial {\theta }_{c}}=0$. We find that

$$\begin{eqnarray}&&\mathrm{ln}[1+{\theta }_{c}({\beta }_{w}-{\beta }_{c})]+({\beta }_{w}-{\beta }_{c}){\theta }_{c}-{\beta }_{w}{\theta }_{h}=-{\beta }_{h}{\theta }_{h}.\end{eqnarray} \tag{ 33 }$$

Since $\mathrm{ln}(1+x)\leqslant x$, we turn the equality (33) into an inequality,

$$\begin{eqnarray}&&({\beta }_{w}-{\beta }_{h}){\theta }_{h}\leqslant 2({\beta }_{w}-{\beta }_{c}){\theta }_{c},\end{eqnarray} \tag{ 34 }$$

or

$$\begin{eqnarray}&&{\theta }_{c}(2{\beta }_{c}-{\beta }_{w}-{\beta }_{h})\leqslant {\theta }_{w}({\beta }_{h}-{\beta }_{w}).\end{eqnarray} \tag{ 35 }$$

The MPE is thus bounded by

$$\begin{eqnarray}&&{\eta }^{* }=\displaystyle \frac{{\theta }_{c}}{{\theta }_{w}}\leqslant \displaystyle \frac{{\beta }_{h}-{\beta }_{w}}{2{\beta }_{c}-{\beta }_{w}-{\beta }_{h}}=\displaystyle \frac{{\eta }_{c}}{{\eta }_{c}+2}.\end{eqnarray} \tag{ 36 }$$

Since ${\beta }_{w}\leqslant {\beta }_{h}$, we get the bound

$$\begin{eqnarray}&&{\eta }^{* }\leqslant \displaystyle \frac{1}{2}\displaystyle \frac{{\beta }_{h}-{\beta }_{w}}{{\beta }_{c}-{\beta }_{h}}=\displaystyle \frac{1}{2}{\eta }_{c}.\end{eqnarray} \tag{ 37 }$$

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 $-\tfrac{1}{4}{\eta }_{c}^{2}$, which is different than the value obtained for low dissipation engines [63].

We repeat this exercise using functions ${M}_{\nu }(\omega )$ of a finite width δ. For example, we start from equation (25), expand it to third order in δ, and solve $\tfrac{\partial {J}_{c}}{\partial {\theta }_{c}}=0$ to obtain a relation for ${\theta }_{c}$ which maximizes the cooling current. Using again the fact that ${{\rm{e}}}^{x}\geqslant 1+x$, we find that

$$\begin{eqnarray}&&\displaystyle \frac{{\theta }_{c}}{{\theta }_{w}}\leqslant \displaystyle \frac{1}{2}{\eta }_{c}+\displaystyle \frac{1}{{\theta }_{w}}\displaystyle \frac{1}{2{\beta }_{c}-{\beta }_{w}-{\beta }_{h}}\mathrm{ln}\left[\displaystyle \frac{2+\tfrac{1}{3}{\beta }_{w}^{2}{\delta }^{2}}{2+\tfrac{1}{3}{\beta }_{h}^{2}{\delta }^{2}}\right].\end{eqnarray} \tag{ 38 }$$

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 ${\beta }_{h}\to {\beta }_{w}\to 0$, see equation (26). The right-hand side then collapses to ${\eta }_{c}/2$. We therefore conclude that in the linear response limit the maximum power efficiency is upper bounded by ${\eta }_{c}/2$, independent of the details of the model.

In this section, we examine the performance of the qubit-QAR beyond the ideal (or close to ideal) limit as explored in section 3. We perform simulation directly by using equations (10) and the analogous expression for J_w. For the FBCF we use the following model ($\omega \geqslant 0$),

$$\begin{eqnarray}&&{M}_{\nu }(\omega )={\omega }^{p}[{n}_{\nu }(\omega )+1][H(\omega -{\theta }_{\nu }+{\delta }_{\nu })-H(\omega -{\theta }_{\nu }-{\delta }_{\nu })].\end{eqnarray} \tag{ 39 }$$

Here, ${\delta }_{\nu }$ is not necessarily small, ${n}_{\nu }(\omega )={[\exp ({\beta }_{\nu }\omega )-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 $\tfrac{{M}_{\nu }(\omega )}{{M}_{\nu }(-\omega )}={{\rm{e}}}^{{\beta }_{\nu }\omega }$. The model is motivated from physical grounds in section 5. In simulations we use p = 1, thus one should add a prefactor to organize 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 ${\delta }_{\nu }$).

As a concrete example, we assume that the cold bath has a finite-fixed width of ${\delta }_{c}=0.2$. The other two windows are tuned with ${\delta }_{w}=\delta$ and ${\delta }_{h}=2\delta$. For a schematic representation, see figure 2. The motivation behind this choice has been to show that the system can act as a QAR even when the spectral windows of the reservoirs overlap, as long as the work bath does not acquire high frequency modes, and the cold bath is characterized by a limited range window. Nevertheless, our simulations demonstrate that the device is quite robust, and that it can operate as a QAR for varied choices for ${\delta }_{c,h,w}$.

Zoom In Zoom Out Reset image size

In figure 3, we display the cooling current and the cooling COP for the model of figure 2. The current ${J}_{c}/{C}_{3}$ is dimensionless, with ${C}_{3}\propto {\gamma }^{2}$, γ as the system-bath interaction energy, see text below equation (10). To reach the current in physical units of e.g. Joule s⁻¹, 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, equation (15), cooling takes place when ${\beta }_{h}\geqslant \tfrac{{\theta }_{c}}{{\theta }_{h}}({\beta }_{c}-{\beta }_{w})+{\beta }_{w}$. In our parameters, this reduces to ${\beta }_{h}\geqslant 0.9/3+0.1=0.4$. This estimate qualitatively agrees with simulations at small δ.

Zoom In Zoom Out Reset image size

The width parameter is important as well, and we find that refrigeration takes place as long as $\delta \lesssim 1.6$. Note that when $\delta \geqslant 1.8$, the functions ${M}_{c}(\omega )$ and ${M}_{w}(\omega )$ 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 $\delta \geqslant 2/3$, yet their overlap does not prohibit the cooling process. Overall, temperature, ${\theta }_{\nu }$, and the width parameters ${\delta }_{\nu }$ intermingle in the expression for cooling. As a result, the cooling window depends on δ in a non-trivial manner. As an additional comment, in figure 3, we consistently find that $\eta \lt {\theta }_{c}/{\theta }_{w}$, which is 1/2 in our parameters. In other words, the QAR performs below the efficiency obtained in equation (27)—valid to linear order in δ. However, when using a narrower window for the cold bath, e.g. ${\delta }_{c}\sim 0.1$, we find that in fact the cooling COP can exceed the value 1/2, yet obviously it still lies below the Carnot bound.

In figure 4, we display slices of the cooling current as a function of the width δ and temperature ${\beta }_{h}$, taken from the contour plot, figure 3. For small value of δ the cooling current grows linearly. However, at a certain point ($\delta =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 ${\beta }_{h}$ is monotonic.

Zoom In Zoom Out Reset image size

The role of a nonzero gap ${\omega }_{0}\ne 0$ on the operation of the qubit-QAR is presented in figure 5, using again parameters as in figure 2. In panel (a), we set $\delta =1$. We find that within the cooling window, the performance of the QAR is intact for small ${\omega }_{0}$, but it deteriorates and eventually disappears once ${\omega }_{0}$ is comparable to differences of the central frequencies, ${\theta }_{w}-{\theta }_{c}$. Outside the cooling window, for large δ, in panel (b) we reveal a non-trivial non-monotonic behavior of J_c with ${\omega }_{0}$, with a large value of ${\omega }_{0}$ manifesting a cooling function that is missing in the degenerate case. Overall, as long as ${\omega }_{0}\ll \delta$, 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.

Zoom In Zoom Out Reset image size

Before concluding, we discuss physical models corresponding to the model Hamiltonian (4) and the functions ${M}_{\nu }(\omega )$. We assume that the baths include collections of harmonic oscillators, ${\hat{H}}_{B,\nu }={\sum }_{j}{\omega }_{j,\nu }{\hat{b}}_{j,\nu }^{\dagger }{\hat{b}}_{j,\nu }$, and that bath operators, which are coupled to the system, are of the form

$$\begin{eqnarray}&&{\hat{B}}_{\nu }=\displaystyle \sum _{j}\displaystyle \frac{{\lambda }_{j,\nu }}{{\omega }_{j,\nu }}({\hat{b}}_{j,\nu }^{\dagger }+{\hat{b}}_{j,\nu }).\end{eqnarray} \tag{ 40 }$$

Returning to the Hamiltonian (4), in this model the displacements from equilibrium of the baths' oscillators are coupled to the transition operator ${\hat{\sigma }}_{x}$, thus the model generalizes the Caldeira–Leggett model [66]. The two-time correlation functions are thus given by

$$\begin{eqnarray}&&{{M}}_{{\nu }}({t})=\langle {\hat{{B}}}_{{\nu }}({t}){\hat{{B}}}_{{\nu }}(0)\rangle =\displaystyle \sum _{{j}}\displaystyle \frac{{{\lambda }}_{{j},{\nu }}^{2}}{{{\omega }}_{{j},{\nu }}^{2}}[{{n}}_{{\nu }}({{\omega }}_{{j},{\nu }}){{\rm{e}}}^{{\rm{i}}{{\omega }}_{{j},{\nu }}{t}}+({{n}}_{{\nu }}({{\omega }}_{{j},{\nu }})+1){{\rm{e}}}^{-{\rm{i}}{{\omega }}_{{j},{\nu }}{t}}],\end{eqnarray} \tag{ 41 }$$

with ${n}_{\nu }(\omega )={[{{\rm{e}}}^{{\beta }_{\nu }\omega }-1]}^{-1}$ as the Bose–Einstein distribution function. In the frequency domain,

$$\begin{eqnarray}&&{M}_{\nu }(\omega )=2\pi \displaystyle \sum _{j}\displaystyle \frac{{\lambda }_{j,\nu }^{2}}{{\omega }_{j,\nu }^{2}}[{n}_{\nu }({\omega }_{j,\nu })\delta (\omega +{\omega }_{j,\nu })+({n}_{\nu }({\omega }_{j,\nu })+1)\delta (\omega -{\omega }_{j,\nu })].\end{eqnarray} \tag{ 42 }$$

We now define the spectral density function, which is related to the density of states of the bath,

$$\begin{eqnarray}&&{g}_{\nu }(\omega )=2\pi \displaystyle \sum _{j}\displaystyle \frac{{\lambda }_{j,\nu }^{2}}{{\omega }_{j,\nu }^{2}}\delta (\omega -{\omega }_{j,\nu }),\end{eqnarray} \tag{ 43 }$$

and find that,

$$\begin{eqnarray}&&{M}_{\nu }(\omega )=[{n}_{\nu }(\omega )+1]{g}_{\nu }(\omega )\end{eqnarray} \tag{ 44 }$$

with ${g}_{\nu }(\omega )$ analytically continued to the complete real axis as an odd function. According to equation (44), the function ${M}_{\nu }(\omega )$ is linearly related to the spectral density function ${g}_{\nu }(\omega )$ 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,

$$\begin{eqnarray}&&\hat{H}=\displaystyle \frac{{\omega }_{0}}{2}{\hat{\sigma }}_{z}+\displaystyle \frac{{\rm{\Delta }}}{2}{\hat{\sigma }}_{x}+\displaystyle \sum _{\nu ,j}{\omega }_{j,\nu }{\hat{b}}_{j,\nu }^{\dagger }{\hat{b}}_{j,\nu }+{\hat{\sigma }}_{z}\displaystyle \sum _{\nu ,j}{\lambda }_{j,\nu }({\hat{b}}_{j,\nu }^{\dagger }+{\hat{b}}_{j,\nu }),\end{eqnarray} \tag{ 45 }$$

and transform it to the displaced bath-oscillators basis using the small polaron transformation [67], ${\hat{H}}_{p}={\hat{U}}^{\dagger }\hat{H}\hat{U}$, $\hat{U}={{\rm{e}}}^{{\rm{i}}{\hat{\sigma }}_{z}\hat{{\rm{\Omega }}}/2}$,

$$\begin{eqnarray}&&{\hat{H}}_{p}=\displaystyle \frac{{\omega }_{0}}{2}{\hat{\sigma }}_{z}+\displaystyle \frac{{\rm{\Delta }}}{2}({\hat{\sigma }}_{+}{{\rm{e}}}^{{\rm{i}}\hat{{\rm{\Omega }}}}+{\hat{\sigma }}_{-}{{\rm{e}}}^{-{\rm{i}}\hat{{\rm{\Omega }}}})+\displaystyle \sum _{\nu ,j}{\omega }_{j,\nu }{\hat{b}}_{j,\nu }^{\dagger }{\hat{b}}_{j,\nu }.\end{eqnarray} \tag{ 46 }$$

Here, ${\hat{\sigma }}_{\pm }=\tfrac{1}{2}({\hat{\sigma }}_{x}\pm {\rm{i}}{\hat{\sigma }}_{y})$ are the auxiliary Pauli matrices, $\hat{{\rm{\Omega }}}={\sum }_{\nu }{\hat{{\rm{\Omega }}}}_{\nu }$, ${\hat{{\rm{\Omega }}}}_{\nu }=2{\rm{i}}{\sum }_{j}\tfrac{{\lambda }_{j,\nu }}{{\omega }_{j}}({\hat{b}}_{j,\nu }^{\dagger }-{\hat{b}}_{j,\nu })$. It can be shown that under the so-called noninteracting blip approximation [68–70], the population dynamics follows equation (5) with the time correlation function

$$\begin{eqnarray}&&{M}_{\nu }(t)=\langle {{\rm{e}}}^{{\rm{i}}{\rm{\Omega }}(t)}{{\rm{e}}}^{-{\rm{i}}{\rm{\Omega }}(0)}\rangle ={{\rm{e}}}^{-{Q}_{\nu }(t)}.\end{eqnarray} \tag{ 47 }$$

The average is performed with respect to the initial, product thermal state of the baths. The function $Q(t)={\sum }_{\nu }{Q}_{\nu }(t)$ is complex, with real and imaginary components, ${Q}_{\nu }(t)={Q}_{\nu }^{{\prime} }(t)+{\rm{i}}{Q}{{\prime\prime} }_{\nu }(t)$,

$$\begin{eqnarray}{Q}_{\nu }^{{\prime\prime} }(t) & = & 2{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{g}_{\nu }(\omega )}{\pi }\sin (\omega t){\rm{d}}\omega ,\\ {Q}_{\nu }^{{\prime} }(t) & = & 2{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{g}_{\nu }(\omega )}{\pi }[1-\cos (\omega t)][1+2{n}_{\nu }(\omega )]{\rm{d}}\omega .\end{eqnarray} \tag{ 48 }$$

The multi-bath spin-boson example clearly demonstrates that non-additivity of the interaction Hamiltonian embodies strong coupling effects. Using equation (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 λ [38]. Nevertheless, to realize a QAR we need to design the function ${M}_{\nu }(t)$, which is the central object in the expressions for the population and energy current. For the spin-boson model, this function, (see equations (47) and (48)), depends in a nonlinear manner on the reservoirs' density of states, thus it is not immediately clear how to physically-rationally engineer it to obtain the form (16).

What is the smallest possible fridge? Using a quantum master equation in Lindblad form with additive dissipators, it was argued in [71, 72] 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 achieve refrigeration if (i) the reservoirs are engineered, and are allowed to exchange energy within a highly restricted spectral window, (ii) a resonant condition ${\theta }_{c}+{\theta }_{w}={\theta }_{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 showed that the optimal behavior of our model—with inseparable interaction of the baths—is precisely analogous to that derived for weakly coupled three-level or three-qubit models, in which the subsystem itself (energy structure and transitions) is engineered, rather than the spectral functions of the baths as in our case. 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 equation (11), we evaluate the cooling current when the cold and hot FBCF take the form of a Dirac delta function as in equation (13), with ${\theta }_{h}\gt {\theta }_{c}\gt 0$, but the work reservoir's FBCF is simply a constant, ${M}_{w}(\omega )={\epsilon }_{w}$ for $\omega \geqslant 0$. We then obtain a modified cooling window (compare to equation (14)),

$$\begin{eqnarray}&&{\theta }_{c}({{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}({\theta }_{h}-{\theta }_{c})}-{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}})+{\theta }_{c}({{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}-{{\rm{e}}}^{-{\beta }_{w}({\theta }_{c}+{\theta }_{h})})\geqslant 0.\end{eqnarray} \tag{ 49 }$$

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

$$\begin{eqnarray}&&{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}({\theta }_{h}-{\theta }_{c})}-{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}(1+A)\geqslant 0,\end{eqnarray} \tag{ 50 }$$

where $A\equiv {{\rm{e}}}^{-{\beta }_{w}({\theta }_{c}+{\theta }_{h})}{{\rm{e}}}^{{\beta }_{h}{\theta }_{h}}-{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}$. After some manipulations, we arrive at

$$\begin{eqnarray}&&\displaystyle \frac{{\theta }_{c}}{{\theta }_{h}}\leqslant \displaystyle \frac{{\beta }_{h}-{\beta }_{w}}{{\beta }_{c}-{\beta }_{w}}-\displaystyle \frac{1}{{\theta }_{h}({\beta }_{c}-{\beta }_{w})}\mathrm{ln}(1+A).\end{eqnarray} \tag{ 51 }$$

Since $A\gt 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 figure 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 [25] or the polaron-transformed master equation [54, 55, 73, 74], 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 [22–25, 27–29], 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.

DS acknowledges support from an NSERC Discovery Grant and the Canada Research Chair program. The work of AM was supported by the CQIQC summer fellowship at the University of Toronto. GS gratefully acknowledges financial support by the DFG (SCHA 1646/3-1, GRK 1558). Discussions with Javier Cerrillo are gratefully acknowledged.

We investigate here the cooling performance of a qubit of a finite energy gap ${\omega }_{0}$. For simplicity, we assume the following functions ($\omega \gt 0$). The negative-frequency branch is decorated by detailed-balance (thermal) factors,

$$\begin{eqnarray}{M}_{c}(\omega ) & = & {\epsilon }_{c}\delta (\omega -{\theta }_{c}),\\ {M}_{h}(\omega ) & = & {\epsilon }_{h}[\delta (\omega -{\theta }_{h}-{\omega }_{0})+\delta (\omega -{\theta }_{h}+{\omega }_{0})],\\ {M}_{w}(\omega ) & = & {\epsilon }_{w}\delta (\omega -{\theta }_{w}).\end{eqnarray} \tag{ A1 }$$

Here, ${\epsilon }_{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, ${\theta }_{c}+{\theta }_{w}={\theta }_{h}$. We further assume that ${\omega }_{0}$ is much smaller than ${\theta }_{\nu },\nu =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 equation (5). The three-bath convoluted decay rate is

$$\begin{eqnarray}&&M({\omega }_{0})={\int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{M}_{c}({\omega }_{1}){\int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){M}_{w}({\omega }_{0}-{\omega }_{1}-{\omega }_{2}).\end{eqnarray} \tag{ A2 }$$

We break the integral into four contributions. For ${\omega }_{1}\geqslant 0$, ${\omega }_{2}\geqslant 0$,

$$\begin{eqnarray}M({\omega }_{0}) & = & {\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{1}{M}_{c}({\omega }_{1}){\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){M}_{w}({\omega }_{0}-{\omega }_{1}-{\omega }_{2})\\ & = & 0,\end{eqnarray} \tag{ A3 }$$

since we assume that ${\omega }_{0}$ is sufficiently small, such that ${\theta }_{c}+{\theta }_{h}-2{\omega }_{0}\gt {\theta }_{w}$. Similarly, when ${\omega }_{1}\leqslant 0$, ${\omega }_{2}\leqslant 0$,

$$\begin{eqnarray}M({\omega }_{0}) & = & {\displaystyle \int }_{-\infty }^{0}{\rm{d}}{\omega }_{1}{M}_{c}({\omega }_{1}){\displaystyle \int }_{-\infty }^{0}{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){M}_{w}({\omega }_{0}-{\omega }_{1}-{\omega }_{2})\\ & = & 0.\end{eqnarray} \tag{ A4 }$$

We receive a finite contribution when ${\omega }_{1}\geqslant 0$ and ${\omega }_{2}\leqslant 0$,

$$\begin{eqnarray}M({\omega }_{0}) & = & {\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{1}{M}_{c}({\omega }_{1}){\displaystyle \int }_{-\infty }^{0}{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){M}_{w}({\omega }_{0}-{\omega }_{1}-{\omega }_{2})\\ & = & {\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{1}{M}_{c}({\omega }_{1}){\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){{\rm{e}}}^{-{\beta }_{h}{\omega }_{2}}{M}_{w}({\omega }_{0}-{\omega }_{1}+{\omega }_{2})\\ & = & {{\rm{e}}}^{-{\beta }_{h}({\theta }_{h}-{\omega }_{0})},\end{eqnarray} \tag{ A5 }$$

as well as for ${\omega }_{1}\leqslant 0$, and ${\omega }_{2}\gt 0$,

$$\begin{eqnarray}M({\omega }_{0}) & = & {\displaystyle \int }_{-\infty }^{0}{\rm{d}}{\omega }_{1}{M}_{c}({\omega }_{1}){\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){M}_{w}({\omega }_{0}-{\omega }_{1}-{\omega }_{2})\\ & = & {\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{1}{M}_{c}({\omega }_{1}){{\rm{e}}}^{-{\beta }_{c}{\omega }_{1}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){M}_{w}({\omega }_{0}+{\omega }_{1}-{\omega }_{2})\\ & = & {\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{1}{M}_{c}({\omega }_{1}){{\rm{e}}}^{-{\beta }_{c}{\omega }_{1}}{\displaystyle \int }_{0}^{\infty }{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){M}_{w}({\omega }_{2}-{\omega }_{1}-{\omega }_{0}){{\rm{e}}}^{-{\beta }_{w}({\omega }_{2}-{\omega }_{1}-{\omega }_{0})}\\ & = & {{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}.\end{eqnarray} \tag{ A6 }$$

Therefore, we find that $M({\omega }_{0})={{\rm{e}}}^{-{\beta }_{h}({\theta }_{h}-{\omega }_{0})}+{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}$. Similarly, $M(-{\omega }_{0})={{\rm{e}}}^{-{\beta }_{h}({\theta }_{h}+{\omega }_{0})}+{{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}$. In steady state, the populations follow ${p}_{g}=\tfrac{M({\omega }_{0})}{M({\omega }_{0})+M(-{\omega }_{0})}$ and ${p}_{e}=\tfrac{M(-{\omega }_{0})}{M({\omega }_{0})+M(-{\omega }_{0})}$. Our results reduce to ${p}_{g,e}=1/2$ when ${\omega }_{0}=0$.

We now turn to the expression for the cooling current, equation (10). It is easy to calculate the integrals in this equation since there is only one extra ${\omega }_{1}$ term, which is trivial to account for under (A1). We evaluate the two integrals,

$$\begin{eqnarray}{I}_{1} & \equiv & -{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{\omega }_{1}{M}_{c}({\omega }_{1}){M}_{h}({\omega }_{2}){M}_{w}({\omega }_{0}-{\omega }_{1}-{\omega }_{2})\\ & = & {\theta }_{c}({{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}-{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}{{\rm{e}}}^{{\beta }_{h}{\omega }_{0}}),\end{eqnarray} \tag{ A7 }$$

and

$$\begin{eqnarray}{I}_{2} & \equiv & {\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{\omega }_{1}{M}_{c}(-{\omega }_{1}){M}_{h}(-{\omega }_{2}){M}_{w}(-{\omega }_{0}+{\omega }_{1}+{\omega }_{2})\\ & = & {\theta }_{c}({{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}-{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}{{\rm{e}}}^{-{\beta }_{h}{\omega }_{0}}).\end{eqnarray} \tag{ A8 }$$

In total, the cooling current is

$$\begin{eqnarray}{J}_{c} & = & -{p}_{{\rm{e}}}{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{\omega }_{1}{M}_{c}({\omega }_{1}){M}_{h}({\omega }_{2}){M}_{w}({\omega }_{0}-{\omega }_{1}-{\omega }_{2})\\ & & +{p}_{g}{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{\omega }_{1}{M}_{c}(-{\omega }_{1}){M}_{h}(-{\omega }_{2}){M}_{w}(-{\omega }_{0}+{\omega }_{1}+{\omega }_{2})\\ & = & \displaystyle \frac{{\theta }_{c}}{M({\omega }_{0})+M(-{\omega }_{0})}({{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}+{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}{{\rm{e}}}^{-{\beta }_{h}{\omega }_{0}})({{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}-{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}{{\rm{e}}}^{{\beta }_{h}{\omega }_{0}})\\ & & +\displaystyle \frac{{\theta }_{c}}{M({\omega }_{0})+M(-{\omega }_{0})}({{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}+{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}{{\rm{e}}}^{{\beta }_{h}{\omega }_{0}})({{\rm{e}}}^{-{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-{\beta }_{w}{\theta }_{w}}-{{\rm{e}}}^{-{\beta }_{h}{\theta }_{h}}{{\rm{e}}}^{-{\beta }_{h}{\omega }_{0}})\\ & = & \displaystyle \frac{2{\theta }_{c}}{M({\omega }_{0})+M(-{\omega }_{0})}({{\rm{e}}}^{-2{\beta }_{c}{\theta }_{c}}{{\rm{e}}}^{-2{\beta }_{w}{\theta }_{w}}-{{\rm{e}}}^{-2{\beta }_{h}{\theta }_{h}}).\end{eqnarray} \tag{ A9 }$$

The cooling window, ${J}_{c}\geqslant 0$, precisely follows equation (14)—obtained for the case ${\omega }_{0}=0$.

As we showed in the main text, functions ${M}_{\nu }(\omega )$ of restricted range support the cooling function in a qubit-QAR when the interaction of the qubit with the separate baths is made non-additive. We now examine examples with continuous (soft cutoff) functions, without any filtering.

First, we assume a unimodal Gaussian form for the FBCF, ${M}_{\nu }(\omega )=\sqrt{\tfrac{\pi }{{E}_{r}^{\nu }{T}_{\nu }}}\exp \left(-\tfrac{{(\omega -{E}_{r}^{\nu })}^{2}}{4{E}_{r}^{\nu }{T}_{\nu }}\right)$, which describes bath-induced rates in the spin-boson model at high temperature. Specifically, performing a short time expansion of equation (48), we get ${Q}_{\nu }^{{\prime\prime} }(t)=\tfrac{2}{\pi }t\int {\rm{d}}\omega {g}_{\nu }(\omega )\omega$ and ${Q}_{\nu }^{{\prime} }(t)=\tfrac{2{T}_{\nu }}{\pi }{t}^{2}\int {\rm{d}}\omega {g}_{\nu }(\omega )\omega$. Defining ${E}_{r}^{\nu }\equiv \tfrac{2}{\pi }\int {\rm{d}}\omega {g}_{\nu }(\omega )\omega$, we arrive at ${Q}_{\nu }^{{\prime\prime} }(t)={E}_{r}^{\nu }t$ and ${Q}_{\nu }^{{\prime} }(t)={E}_{r}^{\nu }{T}_{\nu }{t}^{2}$. The energy scale ${E}_{r}^{\nu }$ is referred to as the 'reorganization energy', and it represents the strength of the system-bath interaction. The unimodal Gaussian function satisfies the detailed balance relation. It is important to note that the width of this Gaussian function is determined by the temperature. We will return to this point in equation (B4).

The current from the cold bath is given by equation (11), which solves to

$$\begin{eqnarray}{J}_{c} & = & -\displaystyle \frac{1}{2}{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{\omega }_{1}{M}_{c}({\omega }_{1}){M}_{h}({\omega }_{2}){M}_{w}(-{\omega }_{1}-{\omega }_{2})\\ & & +\displaystyle \frac{1}{2}{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\rm{d}}{\omega }_{2}{\omega }_{1}{M}_{c}(-{\omega }_{1}){M}_{h}(-{\omega }_{2}){M}_{w}({\omega }_{1}+{\omega }_{2})\\ & = & -\displaystyle \frac{1}{2}{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{1}{\omega }_{1}{M}_{c}({\omega }_{1}){\displaystyle \int }_{-\infty }^{\infty }{\rm{d}}{\omega }_{2}{M}_{h}({\omega }_{2}){M}_{w}({\omega }_{1}+{\omega }_{2})[{{\rm{e}}}^{-{\beta }_{w}({\omega }_{1}+{\omega }_{2})}-{{\rm{e}}}^{-{\beta }_{c}{\omega }_{1}}{{\rm{e}}}^{-{\beta }_{h}{\omega }_{2}}]\\ & = & 4{E}_{r}^{c}{\pi }^{\tfrac{5}{2}}\displaystyle \frac{{E}_{r}^{h}({T}_{c}-{T}_{h})+{E}_{r}^{w}({T}_{c}-{T}_{w})}{{({E}_{r}^{c}{T}_{c}+{E}_{r}^{h}{T}_{h}+{E}_{r}^{w}{T}_{w})}^{\tfrac{3}{2}}}\exp \left(-\displaystyle \frac{{({E}_{r}^{c}+{E}_{r}^{h}+{E}_{r}^{w})}^{2}}{4({E}_{r}^{c}{T}_{c}+{E}_{r}^{h}{T}_{h}+{E}_{r}^{w}{T}_{w})}\right).\end{eqnarray} \tag{ B1 }$$

Since this expression is always negative, the system cannot act as a QAR. We can rationalize this result as follows. Let us define an effective temperature

$$\begin{eqnarray}&&{T}^{* }\equiv \displaystyle \frac{{E}_{r}^{c}{T}_{c}+{E}_{r}^{h}{T}_{h}+{E}_{r}^{w}{T}_{w}}{{E}_{r}^{c}+{E}_{r}^{h}+{E}_{r}^{w}},\end{eqnarray} \tag{ B2 }$$

and the total interaction-reorganization energy, ${E}_{r}\equiv {\sum }_{\nu }{E}_{r}^{\nu }$. We can then recast J_c of equation (B1) as follows,

$$\begin{eqnarray}&&{J}_{c}\propto \displaystyle \frac{1}{\sqrt{{T}^{* }{E}_{r}}}\displaystyle \frac{({T}_{c}-{T}^{* })}{{T}^{* }}{{\rm{e}}}^{-{E}_{r}/4{T}^{* }}.\end{eqnarray} \tag{ B3 }$$

Since ${T}^{* }\gt {T}_{c}$, cooling is prohibited. We now understand that, since the reservoirs can absorb and emit all frequency components, the system (qubit) can in fact be characterized by an effective interaction energy and an effective temperature, the latter is always greater than T_c due to the cumulative effect of the other baths. As a result, it is impossible to extract energy from the cold bath. This argument suggests that finite-range hard-cutoff functions are beneficial for the design of a QAR, as we explain in the main body of the paper.

One should note that we cannot derive the Dirac delta function model (13) from the unimodal Gaussian functions. Recall that ${M}_{\nu }(\omega )=\sqrt{\tfrac{\pi }{{E}_{r}^{\nu }{T}_{\nu }}}\exp \left(-\tfrac{{(\omega -{E}_{r}^{\nu })}^{2}}{4{E}_{r}^{\nu }{T}_{\nu }}\right)$. Let us first consider the low temperature limit,

$$\begin{eqnarray}&&\sqrt{\displaystyle \frac{\pi }{{E}_{r}^{\nu }{T}_{\nu }}}\exp \left(-\displaystyle \frac{{(\omega -{E}_{r}^{\nu })}^{2}}{4{E}_{r}^{\nu }{T}_{\nu }}\right)\mathop{\longrightarrow }\limits^{{T}_{\nu }\to 0}2\pi \delta (\omega -{E}_{r}^{\nu }).\end{eqnarray} \tag{ B4 }$$

This form cannot realize a QAR: as the temperature is approaching zero, the negative frequency (detailed balance) component of the spectral function is shrinking. Nevertheless, this component is essential to realize a QAR. Recall that we must absorb energy from the cold and work baths, and release it into the hot bath. Another limit approaching a Dirac delta function is weak interaction,

$$\begin{eqnarray}&&\sqrt{\displaystyle \frac{\pi }{{E}_{r}^{\nu }{T}_{\nu }}}\exp \left(-\displaystyle \frac{{(\omega -{E}_{r}^{\nu })}^{2}}{4{E}_{r}^{\nu }{T}_{\nu }}\right)\mathop{\longrightarrow }\limits^{{E}_{r}^{\nu }\to 0}2\pi \delta (\omega ).\end{eqnarray} \tag{ B5 }$$

In this case, the three delta functions of the hot, cold and work baths collapse at $\omega =0$. Thus, they fully overlap in frequency, and a QAR operation is lost. We conclude that a unimodal Gaussian-function model as considered in this appendix does not reduce to the delta function model as used in the main part of the paper.

We now suggest a more complex, bimodal Gaussian function,

$$\begin{eqnarray}&&{M}_{\nu }(\omega )=\displaystyle \frac{1}{{\sigma }_{\nu }}[{{\rm{e}}}^{-{(\omega -{\theta }_{\nu })}^{2}/2{\sigma }_{\nu }^{2}}+{{\rm{e}}}^{-{(\omega +{\theta }_{\nu })}^{2}/2{\sigma }_{\nu }^{2}}{{\rm{e}}}^{{\beta }_{\nu }\omega }].\end{eqnarray} \tag{ B6 }$$

Unlike the unimodal function, the width of this Gaussian is introduced here as an additional parameter, independent of temperature. By construction, this FBCF satisfies the detailed balance relation, ${M}_{\nu }(\omega )/{M}_{\nu }(-\omega )={{\rm{e}}}^{{\beta }_{\nu }\omega }$. In figure 6 we exemplify the behavior of the bimodal Gaussian model for the parameters of figure 3. We find that for small width, ${\sigma }_{\nu }/{\theta }_{\nu }\ll 1$, we recover the cooling condition of the ideal design. The QAR performance deteriorates and eventually dissolves as we increase ${\sigma }_{w}$. This is because ${M}_{w}({\theta }_{c}+{\theta }_{h})$ should be sufficiently small here, so as to suppress cooperative energy emission processes, from the work reservoir to the cold and hot thermal baths.

Zoom In Zoom Out Reset image size