Maximum power and corresponding efficiency for two-level heat engines and refrigerators: optimality of fast cycles

TL;DR

This paper demonstrates that the maximum power in a two-level quantum heat engine or refrigerator is universally achieved by a fast cycle with rapid control switching, and the efficiency at maximum power can approach Carnot's limit.

Contribution

It establishes a universal optimal control protocol for maximum power in two-level thermal machines, regardless of microscopic details or operation mode.

Findings

Maximum power is achieved by a fast Otto-cycle like protocol.

Derived a closed-form formula for maximum power.

Efficiency at maximum power can approach Carnot's efficiency.

Abstract

We study how to achieve the ultimate power in the simplest, yet non trivial, model of a thermal machine, namely a two-level quantum system coupled to two thermal baths. Without making any prior assumption on the protocol, via optimal control we show that, regardless of the microscopic details and of the operating mode of the thermal machine, the maximum power is universally achieved by a fast Otto-cycle like structure in which the controls are rapidly switched between two extremal values. A closed formula for the maximum power is derived, and finite-speed effects are discussed. We also analyse the associated efficiency at maximum power (EMP) showing that, contrary to universal results derived in the slow-driving regime, it can approach Carnot's efficiency, no other universal bounds being allowed.

Equations

$$P_{[\mathrm{E}]}=⟨J_{\mathrm{H}}⟩+⟨J_{\mathrm{C}}⟩,$$

$$P_{[\mathrm{A}]}=-⟨J_{\mathrm{C}}⟩,$$

$$\frac{d}{dt}\hat{\rho }=-\frac{i}{\hslash }[\hat{\mathcal{H}},\hat{\rho }{]}_{-}+\sum _{\alpha =\mathrm{H},\mathrm{C}}{\mathcal{D}}_{\alpha }\left[\hat{\rho }\right],$$

$${\mathcal{D}}_{\alpha }\left[\cdots \right]:=\sum _{i=\pm }{\lambda }_{\alpha }(t){\Gamma }_{\alpha }^{(i)}(ϵ(t))({\hat{\sigma }}_i\cdots {\hat{\sigma }}_i^{†}-\frac{1}{2}[{\hat{\sigma }}_i^{†}{\hat{\sigma }}_i,\cdots {]}_{+})$$

$$p_{\mathrm{eq}}^{(\alpha )}(ϵ):=\frac{{\Gamma }_{\alpha }^{(+)}(ϵ)}{{\Gamma }_{\alpha }^{(+)}(ϵ)+{\Gamma }_{\alpha }^{(-)}(ϵ)}=\frac{1}{1+e^{{\beta }_{\alpha }ϵ}}$$

$${\Gamma }_{\alpha }(ϵ):={\Gamma }_{\alpha }^{(+)}(ϵ)+{\Gamma }_{\alpha }^{(-)}(ϵ)$$

$$J_{\alpha }=\mathrm{tr}[\hat{\mathcal{H}}{\mathcal{D}}_{\alpha }\left[\hat{\rho }\right]]=-ϵ(t){\lambda }_{\alpha }(t){\Gamma }_{\alpha }[ϵ(t)](p(t)-p_{\mathrm{eq}}^{(\alpha )}[ϵ(t)]),$$

$$\frac{d}{dt}p(t)=-\sum _{\alpha =\mathrm{H},\mathrm{C}}{\lambda }_{\alpha }(t){\Gamma }_{\alpha }[ϵ(t)](p(t)-p_{\mathrm{eq}}^{(\alpha )}[ϵ(t)]),$$

$$P_{[\nu ]}^{(\max )}=\underset{({ϵ}_{\mathrm{H}},{ϵ}_{\mathrm{C}})\in \mathcal{C}}{\max }\frac{{\Gamma }_{\mathrm{H}}({ϵ}_{\mathrm{H}}){\Gamma }_{\mathrm{C}}({ϵ}_{\mathrm{C}})\left(p_{\mathrm{eq}}^{(\mathrm{H})}({ϵ}_{\mathrm{H}})-p_{\mathrm{eq}}^{(\mathrm{C})}({ϵ}_{\mathrm{C}})\right)}{{\left(\sqrt{{\Gamma }_{\mathrm{H}}({ϵ}_{\mathrm{H}})}+\sqrt{{\Gamma }_{\mathrm{C}}({ϵ}_{\mathrm{C}})}\right)}^2}{\widetilde{ϵ}}_{[\nu ]},$$

$${\Gamma }_{\alpha }^{({\mathrm{F}}_n)}(ϵ)=k_{\alpha }{ϵ}^n,$$

$$\eta (P_{[\mathrm{E}]}^{(\max )})=1-{ϵ}_{\mathrm{C}}^{\ast }/{ϵ}_{\mathrm{H}}^{\ast }=1-\left(1-{\eta }_{\mathrm{c}}\right){ϵ}_{\mathrm{C}}^{\ast }{\beta }_{\mathrm{C}}/{ϵ}_{\mathrm{H}}^{\ast }{\beta }_{\mathrm{H}}.$$

$${\mathrm{C}}_{\mathrm{op}}(P_{[\mathrm{R}]}^{\max })=\frac{{ϵ}_{\mathrm{C}}^{\ast }}{({ϵ}_{\mathrm{H}}^{\ast }-{ϵ}_{\mathrm{C}}^{\ast })}={\left[\frac{{\beta }_{\mathrm{H}}{ϵ}_{\mathrm{H}}^{\ast }}{{\beta }_{\mathrm{C}}{ϵ}_{\mathrm{C}}^{\ast }}(1/{\mathrm{C}}_{\mathrm{op}}^{(\mathrm{c})}+1)-1\right]}^{-1}$$

$${\mathrm{C}}_{\mathrm{op}}(P_{[\mathrm{R}]}^{\max })={\mathrm{C}}_{\mathrm{op}}^{(0)}{\mathrm{C}}_{\mathrm{op}}^{(\mathrm{c})}/(1+{\mathrm{C}}_{\mathrm{op}}^{(0)}+{\mathrm{C}}_{\mathrm{op}}^{(\mathrm{c})}),$$

$$P_{[\mathrm{R}]}^{(\max )}=c_n{k}_{\mathrm{C}}/{\beta }_{\mathrm{C}}^{n+1},{\mathrm{C}}_{\mathrm{op}}(P_{[\mathrm{R}]}^{\max })=0,$$

$$P_{[\mathrm{H}]}^{(\max)}=\frac{k\Delta^{n+1}}{2}\times\left\{\begin{array}[]{lcr}\tanh{\frac{\beta\Delta}{2}},&&\mbox{($\mathrm{F}_{n}$ model)},\\ 1,&&\mbox{($\mathrm{B}_{n}$ model)},\end{array}\right.$$

$$A(t)=-\sum _{\alpha =\mathrm{H},\mathrm{C}}{\lambda }_{\alpha }(t){\Gamma }_{\alpha }[ϵ(t)],B(t)=\sum _{\alpha =\mathrm{H},\mathrm{C}}{\lambda }_{\alpha }(t){\Gamma }_{\alpha }[ϵ(t)]p_{\mathrm{eq}}^{(\alpha )}[ϵ(t)],$$

$$p(t)={\int }_0^t{e}^{{\int }_{t^{\prime }}^{t}A(t^{\prime \prime })d{t}^{\prime \prime }}B(t^{\prime })d{t}^{\prime }+e^{{\int }_0^{t}A(t^{\prime })d{t}^{\prime }}p(0).$$

$$\underset{0}{\overset{t}{\int }}{e}^{{\int }_{t^{\prime }}^{t}A(t^{\prime \prime })d{t}^{\prime \prime }}B(t^{\prime })d{t}^{\prime }=\underset{t-\tau }{\overset{t}{\int }}{e}^{{\int }_{t^{\prime }}^{t}A(t^{\prime \prime })d{t}^{\prime \prime }}B(t^{\prime })d{t}^{\prime }+\underset{0}{\overset{t-\tau }{\int }}{e}^{{\int }_{t^{\prime }}^{t-\tau }A(t^{\prime \prime })d{t}^{\prime \prime }}{e}^{{\int }_{t-\tau }^{t}A(t^{\prime \prime })d{t}^{\prime \prime }}B(t^{\prime })d{t}^{\prime }.$$

$$p(t)=c(t)+d{\int }_0^{t-\tau }{e}^{{\int }_{t^{\prime }}^{t-\tau }A(t^{\prime \prime })d{t}^{\prime \prime }}B(t^{\prime })d{t}^{\prime }+e^{{\int }_0^{t}A(t^{\prime })d{t}^{\prime }}p(0).$$

$$p(t-\tau )\approx {\int }_0^{t-\tau }{e}^{{\int }_{t^{\prime }}^{t-\tau }A(t^{\prime \prime })d{t}^{\prime \prime }}B(t^{\prime })d{t}^{\prime },$$

$$p(t)\approx c(t)+dp(t-\tau ),$$

$${\tau }_{\mathrm{H}}/{\tau }_{\mathrm{C}}=\sqrt{{\Gamma }_{\mathrm{C}}({ϵ}_{\mathrm{C}})/{\Gamma }_{\mathrm{H}}({ϵ}_{\mathrm{H}})}.$$

$$p_{\mathrm{H}}(t)=H{e}^{-{\Gamma }_{\mathrm{H}}t}+p_{\mathrm{eq}}^{(\mathrm{H})},$$

$$p_{\mathrm{H}}({\tau }_{\mathrm{H}})=p_{\mathrm{C}}({\tau }_{\mathrm{H}})$$

$$p_{\mathrm{H}}(0)=p_{\mathrm{C}}({\tau }_{\mathrm{H}}+{\tau }_{\mathrm{C}}).$$

$$\frac{1}{p_{\mathrm{eq}}^{(\mathrm{C})}-p_{\mathrm{eq}}^{(\mathrm{H})}}\left(\matrix{e^{-\Gamma_{\mathrm{H}}\tau_{\mathrm{H}}}&-e^{-\Gamma_{\mathrm{C}}\tau_{\mathrm{H}}}\cr 1&-e^{-\Gamma_{\mathrm{C}}(\tau_{\mathrm{H}}+\tau_{\mathrm{C}})}}\right)\left(\matrix{H\cr C}\right)=\left(\matrix{1\cr 1}\right),$$

$$\left(\text{matrix}H\text{cr}C\right)=\frac{\left(p_{\mathrm{eq}}^{(\mathrm{C})}-p_{\mathrm{eq}}^{(\mathrm{H})}\right)}{\sinh \left[({\Gamma }_{\mathrm{H}}{\tau }_{\mathrm{H}}+{\Gamma }_{\mathrm{C}}{\tau }_{\mathrm{C}})/2\right]}\left(\text{matrix}{e}^{{\Gamma }_{\mathrm{H}}{\tau }_{\mathrm{H}}/2}\sinh \left({\Gamma }_{\mathrm{C}}{\tau }_{\mathrm{C}}/2\right)\text{cr}-e^{{\Gamma }_{\mathrm{C}}{\tau }_{\mathrm{H}}}{e}^{{\Gamma }_{\mathrm{C}}{\tau }_{\mathrm{C}}/2}\sinh \left({\Gamma }_{\mathrm{H}}{\tau }_{\mathrm{H}}/2\right)\right),$$

$$⟨J_{\mathrm{H}}⟩:=\frac{1}{{\tau }_{\mathrm{H}}+{\tau }_{\mathrm{C}}}\underset{{\mathcal{I}}_{\mathrm{H}}}{\int }{J}_{\mathrm{H}}dt=\frac{{ϵ}_{\mathrm{H}}}{{\tau }_{\mathrm{H}}+{\tau }_{\mathrm{C}}}\underset{{\mathcal{I}}_{\mathrm{H}}}{\int }{\dot{p}}_{\mathrm{H}}dt=\frac{{ϵ}_{\mathrm{H}}H}{{\tau }_{\mathrm{H}}+{\tau }_{\mathrm{C}}}\left[e^{-{\Gamma }_{\mathrm{H}}{\tau }_{\mathrm{H}}}-1\right],$$

$$⟨J_{\mathrm{C}}⟩:=\frac{1}{{\tau }_{\mathrm{H}}+{\tau }_{\mathrm{C}}}\underset{{\mathcal{I}}_{\mathrm{C}}}{\int }{J}_{\mathrm{C}}dt=\frac{{ϵ}_{\mathrm{C}}}{{\tau }_{\mathrm{H}}+{\tau }_{\mathrm{C}}}\underset{{\mathcal{I}}_{\mathrm{C}}}{\int }{\dot{p}}_{\mathrm{C}}dt=\frac{{ϵ}_{\mathrm{C}}C}{{\tau }_{\mathrm{H}}+{\tau }_{\mathrm{C}}}{e}^{-{\Gamma }_{\mathrm{C}}{\tau }_{\mathrm{H}}}\left[e^{-{\Gamma }_{\mathrm{C}}{\tau }_{\mathrm{C}}}-1\right],$$

$$⟨J_{H/C}⟩=\pm \frac{{ϵ}_{H/C}}{{\tau }_{\mathrm{H}}+{\tau }_{\mathrm{C}}}\frac{{\Gamma }_{\mathrm{H}}{\tau }_{\mathrm{H}}{\Gamma }_{\mathrm{C}}{\tau }_{\mathrm{C}}}{{\Gamma }_{\mathrm{H}}{\tau }_{\mathrm{H}}+{\Gamma }_{\mathrm{C}}{\tau }_{\mathrm{C}}}\left(p_{\mathrm{eq}}^{(\mathrm{H})}-p_{\mathrm{eq}}^{(\mathrm{C})}\right)\frac{{\left({\Gamma }_{\mathrm{H}}{\tau }_{\mathrm{H}}/2\right)}^{-1}+{\left({\Gamma }_{\mathrm{C}}{\tau }_{\mathrm{C}}/2\right)}^{-1}}{\coth \left({\Gamma }_{\mathrm{H}}{\tau }_{\mathrm{H}}/2\right)+\coth \left({\Gamma }_{\mathrm{C}}{\tau }_{\mathrm{C}}/2\right)}.$$

$$⟨J_{H/C}⟩=\pm {ϵ}_{H/C}\frac{{\Gamma }_{\mathrm{H}}\theta {\Gamma }_{\mathrm{C}}(1-\theta )}{{\Gamma }_{\mathrm{H}}\theta +{\Gamma }_{\mathrm{C}}(1-\theta )}\left(p_{\mathrm{eq}}^{(\mathrm{H})}-p_{\mathrm{eq}}^{(\mathrm{C})}\right).$$

Full text

• June 2019

1 Introduction

Two thermal baths in contact through a working fluid that can be externally driven represent the prototypical setup that has been studied from the origin of thermodynamics up to our days. The energy balance can be described in terms of three quantities: the work extracted from the fluid and the heat exchanged with the hot/cold baths. The fundamental limitations to the inter-conversion of heat into work stem from the concept of irreversibility and are at the core of the second law of thermodynamics. A working medium in contact with two baths at different temperatures is also significant from a practical point of view, since it is the paradigm behind the following specific machines: the heat engine, the refrigerator [1, 2, 3, 4], the thermal accelerator [5], and the heater [5].

Quantum thermodynamics [6, 7, 8] has emerged both as a field of fundamental interest, and as a potential candidate to improve the performance of thermal machines [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. The optimal performance of these systems has been discussed within several frameworks and operational assumptions, ranging from low-dissipation and slow driving regimes[24, 25, 26, 27, 28], to shortcuts to adiabaticity approaches [29, 30, 31, 32], to endoreversible engines [33, 34]. Several techniques have been developed for the optimal control of two-level systems for achieving a variety of goals: from optimizing the speed [35, 36, 37], to generating efficient quantum gates [38, 39], to controlling dissipation [41, 40], and to optimizing thermodynamic performances [42, 43, 44, 45, 46, 47].

The aim of the present paper is to find the optimal strategy to deliver maximum power in all four previously mentioned machines. We perform this optimization in the simplest, yet non trivial, model of a machine which, in the spirit of quantum thermodynamics, is based on a two-level quantum system as working fluid. As opposed to current literature, we explicitly carry out the power maximization without making any assumptions on the operational regime, nor on the speed of the control parameters, nor on the specific coupling between the working fluid and the bath. We find that, if the evolution of the working medium is governed by a Markovian master equation [49, 48], the optimal driving takes a universal form: an infinitesimal Otto-cycle-like structure in which the control parameters must be varied between two extremal values as fast as possible. This is our first main results, described in Eq. (8). Surprisingly, the optimal solution is achieved in the “fast-driving” regime, i.e. when the driving frequency is faster than the typical dissipation rate induced by the baths, which has received little attention in literature [50, 51, 52].

By applying our optimal protocol to heat engines and refrigerators, we find new theoretical bounds on the efficiency at maximum power (EMP). Many upper limits to the EMP, strictly smaller than Carnot’s efficiency, have been derived in literature, such as the Curzon-Ahlborn and Schmiedl-Seifert efficiencies. The Curzon-Ahlborn efficiency emerges in various specific models [53, 55, 54], and it has been derived by general arguments from linear irreversible thermodynamics [56]. The Schmiedl-Seifert efficiency has been proven to be universal in cyclic Brownian heat engines [57] and for any driven system operating in the slow-driving regime [24]. By studying the efficiency of our system at the ultimate power, i.e. in the fast-driving regime, we prove that there is no fundamental upper bound to the EMP. Indeed, we show that the Carnot efficiency is reachable at maximum power through a suitable engineering of the bath couplings. This is our second main results, illustrated in Figs. 2b, 2c and 3. In view of experimental implementations, we assess the impact of finite-time effects on our optimal protocol, finding that the maximum power does not decrease much if the external driving is not much slower than the typical dissipation rate induced by the baths [58, 59]. Furthermore, we apply our optimal protocol to two experimentally accessible models, namely photonic baths coupled to a qubit [60, 61, 62, 22, 63] and electronic leads coupled to a quantum dot [65, 64, 58, 59, 21, 23].

2 Maximum Power.

The setup we consider consists of a two-level quantum system S with energy gap $\epsilon(t)$ that can be externally modulated [66]. As schematically shown in Fig. 1a, the system is placed in thermal contact with two reservoirs, the hot bath H at inverse temperature $\beta_{\rm H}$ and the cold bath C at inverse temperature $\beta_{\rm C}$, respectively characterized by coupling constants $\lambda_{\mathrm{H}}(t)$ and $\lambda_{\mathrm{C}}(t)$ that can be modulated in time. The system can operate in four different modes: i) the heat engine mode [E], where S is used to produce work by extracting heat from H while donating it to C; ii) the refrigerator mode [R], where S is used to extract heat from C; iii) the thermal accelerator mode [A], where S operates to move as much heat as possible to C; iv) the heater mode [H], where we simply use S to deliver as much heat as possible to both H and C. Assuming cyclic modulation of the controls (i.e. of $\epsilon(t)$, $\lambda_{\mathrm{H}}(t)$ and $\lambda_{\mathrm{C}}(t)$) we are interested in maximizing the corresponding averaged output powers of each operating mode, i.e. the quantities

[TABLE]

where $J_{\mathrm{H}}$ and $J_{\mathrm{C}}$ are the instantaneous heat fluxes entering the hot and cold reservoirs respectively, and where the symbol $\braket{\cdots}$ stands for temporal average over a modulation cycle of the controls. To tackle the problem we adopt a Markovian Master Equation (MME) approach [60], namely we write

[TABLE]

where $\hat{\rho}$ is the density matrix of the two-level system at time $t$, $\hat{\mathcal{H}}:=\epsilon(t)\hat{\sigma}_{+}\hat{\sigma}_{-}$ its local Hamiltonian, and

[TABLE]

is the Gorini-Kossakowski-Sudarshan-Lindblad dissipator [48, 49] associated with the bath $\alpha={\rm H},{\rm C}$. We have denoted with $\hat{\sigma}_{+}$ and $\hat{\sigma}_{-}$ the raising and lowering operators of S and with the symbol $[\cdots,\cdots]_{\mp}$ the commutator ($-$) and anti-commutator ($+$) operations. $\cal{D}_{\alpha}$ is characterized by dissipation rates $\Gamma_{\alpha}^{(i=\pm)}(\epsilon)$ and by the dimensionless coupling constant $\lambda_{\alpha}(t)\in[0,1]$ that plays the role of a “switch” control parameter. It is worth noticing that, since $[\hat{\mathcal{H}}(t),\hat{\mathcal{H}}(t^{\prime})]=0$, the MME we employ is valid also in the fast-driving regime, provided that the correlation time of the bath is the smallest timescale in our problem [67]. Therefore, the fast-driving regime is characterized by a control frequency which is faster than the typical dissipation rate, but slower than the inverse correlation time of the bath. Furthermore, we assume that the Hamiltonain $\hat{\mathcal{H}}_{\mathrm{int}}$, describing the system-bath interaction, is such that its expectation value on the Gibbs state of the baths is zero (this is true, for example, for tunnel-like Hamiltonians, where the number of creation/annihilation operators of the bath entering $\hat{\mathcal{H}}_{\mathrm{int}}$ is odd). Such assumption guarantees that no work is necessary to switch on and off the coupling between the system and the baths.

Without assigning any specific value to the dissipation rates, we only require them to obey the detailed balance equation ${\Gamma^{(+)}_{\alpha}(\epsilon)}/{\Gamma^{(-)}_{\alpha}(\epsilon)}=e^{-\beta_{\alpha}\epsilon}$. This ensures that, at constant level spacing $\epsilon$, the system S will evolve into a thermal Gibbs state characterized by an excitation probability

[TABLE]

when in contact only with heat bath $\alpha$. For simplicity, we consider the system to be coupled to one heat bath at the time, i.e. we assume that $\lambda_{\mathrm{H}}(t)+\lambda_{\mathrm{C}}(t)=1$, and that $\lambda_{\alpha}(t)$ can take the values [math] or $1$. As we shall see in the following, this constraint, as well as the possibility of controlling the coupling constants $\lambda_{\alpha}(t)$, is not fundamental to derive our results, at least for those cases where the effective dissipation rate

[TABLE]

of each bath is sufficiently peaked around distinct values. The instantaneous heat flux leaving the thermal bath $\alpha$ can now be expressed as [62]

[TABLE]

where $p(t):=\mathrm{tr}[\hat{\sigma}_{+}\hat{\sigma}_{-}\hat{\rho}(t)]$ is the probability of finding S in the excited state of ${\cal H}$ at time $t$ which obeys the following differential equation

[TABLE]

according to the MME specified above. By explicit integration of (7) we can hence transform all the terms in Eq. (2) into functionals of the controls which can then be optimized with respect to all possible choices of the latter.

As shown in App. A, we find that the protocols which maximize the average power of a fixed physical setup, i.e. at fixed dissipation rates, are cycles performed in the fast-driving regime, i.e. when the driving frequency is faster than the typical dissipation rate. More precisely, the optimal cycle is such that $\epsilon(t)$ instantaneously jumps between two values $\epsilon_{\mathrm{H}}$ and $\epsilon_{\mathrm{C}}$, see Fig. 1b, while being in contact, respectively, only with the hot and cold bath for infinitesimal times $\tau_{\mathrm{H}}$ and $\tau_{\mathrm{C}}$ fulfilling the condition ${\tau_{\mathrm{H}}}/{\tau_{\mathrm{C}}}=\sqrt{{\Gamma_{\mathrm{C}}(\epsilon_{\rm C})}/{\Gamma_{\mathrm{H}}(\epsilon_{\rm H})}}$ [68]. As in Otto cycles considered in literature (see the extensive literature on this topic, e.g. [10, 69, 12, 70, 71, 20]), no heat is transferred during the jumps and no work is done while the system is in contact with the baths. The resulting maximum power averaged over one period can then be cast into the following compact expression (see App. B for details)

[TABLE]

where $\nu=$ E,R,A,H and the quantity $\widetilde{\epsilon}_{[\nu]}$ is given by $\widetilde{\epsilon}_{[E]}=\epsilon_{\mathrm{H}}-\epsilon_{\mathrm{C}}$, $\widetilde{\epsilon}_{\mathrm{[R]}}=-\epsilon_{\mathrm{C}}$, $\widetilde{\epsilon}_{\mathrm{[A]}}=\epsilon_{\mathrm{C}}$, and $\widetilde{\epsilon}_{[H]}=\epsilon_{\mathrm{C}}-\epsilon_{\mathrm{H}}$. In Eq. (8) $\mathcal{C}$ is the range over which the energy gap $\epsilon(t)$ of S is allowed to be varied according to the possible technical limitations associated with the specific implementation of the setup.

Equation (8), which stems from the optimality of the fast-driving regime, is the first main result of the present work. We emphasize that, as opposed to current literature, our closed expression for the maximum power holds for any dissipation rate function $\Gamma_{\mathrm{H/C}}(\epsilon)$. In the following we will apply our result to specific implementations which are relevant experimentally and compute their associated efficiencies at maximum power. In particular we shall consider the case of fermionic ($\mathrm{F}_{n}$) and bosonic ($\mathrm{B}_{n}$) baths with associated effective rates of the form

[TABLE]

with $n\geq 0$ integer and with $k_{\alpha}$ being a coupling strength constant. The fermionic rate (the first of Eq. (9)) for instance can describe two electronic leads, with density of states depending on $n$, tunnel coupled to a single-level quantum-dot [72, 64, 73]; the bosonic one instead is applied in the study of two-level atoms in a dispersive quantum electromagnetic cavity [74].

3 Heat engine mode $[E]$.

It is common belief that the efficiency of a heat engine (work extracted over heat absorbed from the hot bath H), driven at maximum power (EMP), should exhibit a finite gap with respect to the Carnot efficiency $\eta_{\mathrm{c}}:=1-\beta_{\mathrm{H}}/\beta_{\mathrm{C}}$. Indeed, this is corroborated by various results on EMP bounds: the Curzon-Ahlborn EMP $\eta_{\mathrm{CA}}:=1-\sqrt{1-\eta_{\mathrm{c}}}$   emerges in various specific models[53, 57, 27], and it has been derived by general arguments from linear irreversible thermodynamics [56], while the Schmiedl-Seifert EMP $\eta_{\mathrm{SS}}:=\eta_{\mathrm{c}}/(2-\eta_{\mathrm{c}})$ has been proven to be universal for any driven system operating in the slow-driving regime [24]. However, the completely out-of-equilibrium and optimal cycles associated with the values of $P^{(\max)}_{[E]}$ reported in Eq. (8), do not fulfill such assumptions. As a matter of fact, by choosing particular “energy filtering” dissipation rates $\Gamma_{\alpha}(\epsilon)$ (instead of the regular ones given e.g. in Eq. (9)), we can produce configurations which approach Carnot’s efficiency with arbitrary precision while delivering maximum power, proving the lack of any fundamental bound to the EMP. Before discussing this highly not trivial effect, it is worth analyzing the performances associated with the baths models of Eq. (9).

We remind that the efficiency of an Otto cycle heat engine working between the internal energies $\epsilon_{\mathrm{C}}$ and $\epsilon_{\mathrm{H}}$ is given by $\eta=1-{\epsilon_{\mathrm{C}}}/{\epsilon_{\mathrm{H}}}$. Accordingly, indicating with $\epsilon_{\mathrm{H}}^{*}$ and $\epsilon_{\mathrm{C}}^{*}$ the values of the gaps that lead to the maximum of the r.h.s term of Eq. (8), we write the EMP of our scheme as

[TABLE]

In Fig. 2a we report the value of $\eta(P^{(\max)}_{[\mathrm{E}]})$ obtained from (10) for the rates of Eq. (9) for $n=0,1$. By a direct comparison with $\eta_{\mathrm{CA}}$ and $\eta_{\mathrm{SS}}$, one notices that while the second is always respected by our optimal protocol, the first is outperformed at least for the baths $\mathrm{F}_{0}$ and $\mathrm{B}_{0}$, confirming the findings of Refs. [65, 73, 45]. For small temperature differences between the baths, the EMP can be expanded as a power series in $\eta_{\mathrm{c}}$ of the form $a_{1}\eta_{\mathrm{c}}+a_{2}\eta_{\mathrm{c}}^{2}+\cdots.$ It has been shown that $a_{1}=1/2$ is a universal property of low dissipation heat engines [24] and, in this context, $a_{2}=1/8$ is associated with symmetric dissipation coefficients. As explicitly discussed in App. C, we find that also our protocol delivers an efficiency at maximum power with a first order expansion term $a_{1}=1/2$ and with a second order correction $a_{2}=1/8$ achieved if we assume that the two leads are symmetric, i.e. $\Gamma_{\mathrm{H}}(\epsilon,\beta)=\Gamma_{\mathrm{C}}(\epsilon,\beta)$, or if the rates are constants.

We now turn to the possibility of having $\eta(P^{(\max)}_{[\mathrm{E}]})$ arbitrarily close to $\eta_{\mathrm{c}}$. By a close inspection of the second identity of Eq. (10) we notice that one can have $\eta(P^{(\max)}_{[\mathrm{E}]})\simeq\eta_{\mathrm{c}}$ for all those models where the maximum power (see Eq. (8)) is obtained for values of the gaps fulfilling the condition $\epsilon^{*}_{\mathrm{C}}\beta_{\mathrm{C}}\approx\epsilon^{*}_{\mathrm{H}}\beta_{\mathrm{H}}$. Consider hence a scenario where the rates $\Gamma_{\alpha}(\epsilon_{\alpha})$ are such that the power is vanishingly small for all values of $\epsilon_{\alpha}$ except for a windows of width $\sigma$ around a given value $\bar{\epsilon}_{\alpha}$, a configuration that can be used to eliminate the presence of the activation controls $\lambda_{\alpha}(t)$ from the problem. Under the assumption of small enough $\sigma$, we expect the maximization in Eq. (8) to yield $\beta_{\mathrm{C}}\epsilon^{*}_{\mathrm{C}}\approx\beta_{\mathrm{H}}\epsilon^{*}_{\mathrm{H}}$ when the inverse temperature ratio is $\beta_{\mathrm{C}}/\beta_{\mathrm{H}}\approx\bar{\epsilon}_{\mathrm{H}}/\bar{\epsilon}_{\mathrm{C}}$, so that $\eta(P^{(\max)}_{[\mathrm{E}]})\approx\eta_{\mathrm{c}}$. This is indeed evident from Fig. 2b and 2c, where we report the value $\eta(P^{(\max)}_{[\mathrm{E}]})$ as a function of $\eta_{\mathrm{c}}$ (which represents the temperature of the baths) for rates having a Lorentzian shape dependence: by decreasing $\sigma$, the EMP approaches Carnot’s efficiency at $\bar{\eta}_{\mathrm{c}}:=1-\bar{\epsilon}_{\mathrm{C}}/\bar{\epsilon}_{\mathrm{H}}=1/2$ (Fig. 2b), while by tuning the position of the peak of the Lorentzian rates, the EMP can approach Carnot’s efficiency at any given bath temperature configuration $\bar{\eta}_{\mathrm{c}}$ (Fig. 2c). We emphasize that even our system with Lorentzian shaped rates would exhibit an EMP bounded by $\eta_{\mathrm{SS}}$ if operated in the slow-driving regime. The possibility of reaching Carnot’s efficiency at maximum power is thus a characteristic which emerges thanks to the fast-driving regime. Conversely, as $\sigma$ decreases and $\eta(P^{(\max)}_{[\mathrm{E}]})\to\eta_{\mathrm{c}}$, the corresponding maximum power tends to zero (see the inset of Fig. 2b where the maximum power, at $\bar{\eta}_{\mathrm{c}}=1/2$, is plotted as a function of $\sigma$).

4 Refrigerator mode $[{\rm R}]$.

The efficiency of a refrigerator is quantified by the coefficient of performance (COP), i.e. the ratio between the heat extracted from the cold bath and the work done on the system. For an Otto-cylce the COP is given by $\mathrm{C}_{\mathrm{op}}={\epsilon_{\mathrm{C}}}/({\epsilon_{\mathrm{H}}-\epsilon_{\mathrm{C}}})$ which, by replacing the values $\epsilon^{*}_{\mathrm{C}}$, $\epsilon^{*}_{\mathrm{H}}$ that lead to the maximum $P^{(\max)}_{[R]}$ of Eq. (8), yields an associated COP at maximum power (CMP) equal to

[TABLE]

where $\mathrm{C}_{\mathrm{op}}^{(\mathrm{c})}:=\beta_{\mathrm{C}}^{-1}/(\beta_{\mathrm{H}}^{-1}-\beta_{\mathrm{C}}^{-1})$ is the maximum COP dictated by the second law. Remarkably, as in the heat engine case, we can produce configurations which approach $\mathrm{C}_{\mathrm{op}}^{(\mathrm{c})}$ with arbitrary precision while delivering maximum power exploiting the same “energy filtering” dissipation rates. Before discussing this effect we present some universal properties of the CMP and we analyze the performance of the baths models of Eq. (9).

Assuming that the rates depend on the energy and on the temperature through the product $\beta\epsilon$, i.e. $\Gamma_{\alpha}(\epsilon_{\alpha})=\Gamma_{\alpha}(\beta_{\alpha}\epsilon_{\alpha})$ (e.g. the models (9) satisfy this hypothesis for $n=0$, while they do not for $n>0$), we find that the COP at maximum power reduces to the universal family of curves

[TABLE]

where $\mathrm{C}_{\mathrm{op}}^{(0)}$ represents the COP when $\beta_{\mathrm{H}}=\beta_{\mathrm{C}}$. It thus follows that for these models the knowledge of $\mathrm{C}_{\mathrm{op}}(P^{\max}_{[\mathrm{R}]})$ at a single bath temperature configuration identifies unambiguously the COP for all other temperature differences. This feature is in contrast with the heat engine mode since, under the same hypothesis, the EMP at arbitrary temperatures depends on the details of the system.

Consider next the maximum power for the models described Eq. (9). We find that the maximization in Eq. (8) yields $\epsilon^{*}_{\mathrm{H}}\to+\infty$ (and a finite value of $\epsilon^{*}_{\mathrm{C}}$), which implies

[TABLE]

where $c_{n}$ is a dimensionless number which only depends on $n$ for $n>0$, while it is a function of $k_{\mathrm{H}}/k_{\mathrm{C}}$ if $n=0$ (see App. D for details). The fact that the corresponding COP is equal to zero is a direct consequence of the divergent value of $\epsilon^{*}_{\mathrm{H}}$: physically it means that the maximum cooling power [which is finite, see Eq. (13)] is obtained by performing an infinite work, thus by releasing an infinite amount of heat into the hot bath. In the more realistic scenario where there are limitations on our control of the gaps, say $|\epsilon_{\alpha}|\leq\Delta$, the resulting value of $P^{(\max)}_{[\mathrm{R}]}$ will be smaller than in Eq. (13) but the associated COP will be non-zero with a scaling that for large enough $\Delta$ goes as $\mathrm{C}_{\mathrm{op}}(P^{\max}_{[\mathrm{R}]})\propto 1/(\beta_{\mathrm{C}}\Delta)$ (see App. D for details). Equation (13) shows that in all models the maximum cooling power only depends on the temperature $1/\beta_{\mathrm{C}}$ of the cold lead as a simple power law, and it vanishes as $1/\beta_{\mathrm{C}}\to 0$. Intuitively this makes sense since it is harder to refrigerate a colder bath and at $1/\beta_{\mathrm{C}}=0$ there is no energy to extract from the bath. Furthermore, for $n>0$ the properties of the hot bath (i.e. temperature and coupling constant) do not enter the $P^{(\max)}_{[\mathrm{R}]}$ formula.

We now return to the possibility of having the CMP arbitrarily close to $\mathrm{C}_{\mathrm{op}}^{(c)}$. As in the heat engine case, from the second equality of Eq. (11) we see that, if the maximization in Eq. (8) yields values of $\epsilon^{*}_{\mathrm{H}}$ and $\epsilon^{*}_{\mathrm{C}}$ such that $\epsilon_{\mathrm{H}}^{*}\beta_{\mathrm{H}}\approx\epsilon_{\mathrm{C}}^{*}\beta_{\mathrm{C}}$, then $\mathrm{C}_{\mathrm{op}}(P^{\max}_{[\mathrm{R}]})\approx\mathrm{C}_{\mathrm{op}}^{(c)}$. Indeed, as we can see in Fig. 3, we are able to have a CMP close to $\mathrm{C}_{\mathrm{op}}^{(\mathrm{c})}$ at any desired temperature configuration $\bar{C}_{\mathrm{op}}^{(\mathrm{c})}$ by considering appropriately tuned Lorentzian rates (described in Fig. 2).

5 Thermal accelerator $[\rm A]$ and heater $[\rm H]$ modes.

For the physical models described in Eq. (9) it turns out that in order to maximize the heat entering the cold bath, it is more convenient to release heat into both baths ($J_{\mathrm{H}},J_{\mathrm{C}}<0$), rather than extracting heat from the hot bath H and releasing it into the cold bath ($J_{\mathrm{H}}>0,J_{\mathrm{C}}<0$). The thermal accelerator mode $[\rm A]$ thus appears to be useless if we are just interested in maximizing the heat delivered to the cold bath. Accordingly, in the following we shall focus on the heater $[\rm H]$ mode only with a single bath (or equivalently with two baths at the same temperature). Assuming to have some physical limit $|\epsilon|\leq\Delta$ on the way we can control the gap, from Eq. (8) we find

[TABLE]

where $k$ is the coupling constant appearing in Eq. (9). Equation (14) shows that the maximum power diverges as $\Delta\to+\infty$, the exponent of $\Delta$ depending on the density of states associated with the rates. Interestingly, the maximum power that can be delivered to the bath vanishes for high temperatures ($\beta\Delta\ll 1$) in the fermionic models, while it is finite and insensitive to temperature in the bosonic models. This is due to the peculiar rates of the bosonic models which diverge for $\beta\epsilon\ll 1$. On the contrary, for low temperatures ($\beta\Delta\gg 1$) both models yield the same value of $P_{[\mathrm{H}]}^{(\max)}$.

6 Finite-Time Corrections.

The derivation of our main equation  (8) was obtained under the implicit assumption that one could implement infinitesimal control cycles. Yet this hypothesis is not as crucial as it may appear. Indeed the feasibility of an infinitesimal Otto cycle relies on the ability of performing a very fast driving with respect to the typical time scales of the dynamics, a regime that can be achieved in several experimental setups [58, 59]. Furthermore by taking the square-wave protocol shown in Fig. 1b characterized by finite time intervals $\tau_{\mathrm{H}}$ and $\tau_{\mathrm{C}}$ still fulfilling the ratio $\tau_{\mathrm{H}}/\tau_{\mathrm{C}}=\sqrt{{\Gamma_{\mathrm{C}}(\epsilon_{\rm C})}/{\Gamma_{\mathrm{H}}(\epsilon_{\rm H})}}$, we find that, at leading order in $dt$, the maximum power $P^{(\max)}_{[\nu]}(dt)$ only different from the ideal value $P^{(\max)}_{[\nu]}$ of Eq. (8) by a quadratic correction, i.e. $P^{(\max)}_{[\nu]}(dt)\approx(1-\widetilde{\Gamma}_{\mathrm{H}}\widetilde{\Gamma}_{\mathrm{C}}dt^{2}/12)P^{(\max)}_{[\nu]}$, where $\widetilde{\Gamma}_{\alpha}=(\widetilde{\Gamma}\Gamma_{\alpha})^{1/2}$, $\widetilde{\Gamma}=\Gamma_{\mathrm{H}}\Gamma_{\mathrm{C}}/(\sqrt{\Gamma_{\mathrm{H}}}+\sqrt{\Gamma_{\mathrm{C}}})^{2}$, and all rates are computed for $\epsilon_{\mathrm{H}}$ and $\epsilon_{\mathrm{C}}$ that maximize Eq. (8). Besides, even in the regime where $\widetilde{\Gamma}_{\mathrm{H}}dt,\widetilde{\Gamma}_{\mathrm{C}}dt\gg 1$, $P^{(\max)}_{[\nu]}(dt)$ can be shown (see App. B for details) to only decrease as $(\widetilde{\Gamma}_{\mathrm{H}}dt/2)^{-1}+(\widetilde{\Gamma}_{\mathrm{C}}dt/2)^{-1}$, implying that a considerable fraction of $P^{(\max)}_{[\nu]}$ can still be achieved also in this case (e.g. see Fig. 1c where we report the $dt$ dependence of $P^{(\max)}_{[\mathrm{H}]}(dt)$ in the heater mode). On the contrary deviations from Eq. (8) due to finite time corrections in the quenches turns out to be more relevant. These last are first order in the ration between the duration of the quench (now different from [math]) and the period of the protocol $dt\,$ (see App. B for details).

7 Conclusions.

We proved that a cycle switching between two extremal values in the fast-driving regime achieves universally the maximum power and the maximum cooling rate (respectively for a working medium operating as a heat engine or as a refrigerator), regardless of the specific dissipation rates, and we found a general expression for the external control during the cycle. The power advantage of modulating the control fields with rapid adiabatic transformations has been observed in the literature [51, 70, 76] for some specific model and this intuition is in agreement with our general results. We also found that the first coefficient of the expansion in power of $\eta_{C}$ of the EMP is universal while the second one is linked to the symmetry of the dissipation coefficients. This paper enlights that the features mentioned above are valid also strongly out of equilibrium, while already proven in low dissipation [64] and steady state [3] heat engines. If the bath spectral densities can be suitably tailored through energy filters (as for instance in [75]) our protocol allows to reach the Carnot bound at maximum power both operating as a heat engine or refrigerator, although at the cost of a vanishing power. This observation proves the lack of universal upper bounds to the efficiency at maximum power. Finally, a new scaling for the COP of a bath with flat spectral density is shown and a clear dependence of the EMP and the COP at maximum power on the spectral densities of the two thermal baths is established. The results are discussed in detail for some specific models, from flat bosonic and fermionic baths to environments with more complicated spectral densities, and finite driving speed effects are analyzed.

8 Acknowledgments.

We thank G. M. Andolina for useful discussions. This work has been supported by SNS-WIS joint lab “QUANTRA”, by the SNS internal projects “Thermoelectricity in nano-devices”, and by the CNR-CONICET cooperation programme “Energy conversion in quantum, nanoscale, hybrid devices”.

Appendix A. Optimality of infinitesimal Otto cycles

In this appendix we present explicit proof that infinitesimal Otto cycles are optimal for reaching maximum power performances for our two-level setting.

As a preliminary result we clarify that under periodic modulations of the control parameters, the master equation (Eq. (4) of the main text) produces solutions which asymptotically are also periodic. For this purpose let us write Eq. (4) of the main text as $\dot{p}(t)=A(t)p(t)+B(t)$ where, for ease of notation, we introduced the functions

[TABLE]

and consider periodical driving forces such that $A(t+\tau)=A(t)$, $B(t+\tau)=B(t)$ for all $t$. By explicitly integration we get

[TABLE]

Decompose then the integral on the right hand side as

[TABLE]

Notice now that, since $A(t)$ and $B(t)$ are periodic, the quantity $c(t)=\int_{t-\tau}^{t}e^{\int_{t^{\prime}}^{t}A(t^{\prime\prime})dt^{\prime\prime}}B(t^{\prime})dt^{\prime}$ is also periodic with period $\tau$, while $d=e^{\int_{t-\tau}^{t}A(t^{\prime\prime})dt^{\prime\prime}}$ is constant in time. Substituting the previous definitions in Eq. (3) and then in Eq. (2) we find

[TABLE]

In the asymptotic limit where the initial condition $p(0)$ has been completely forgotten (since $A(t)\leq 0$ at all times, the contribution of the initial condition decays exponentially), Eq. (2) gives

[TABLE]

which substituted in Eq. (4) allows us to write

[TABLE]

where we neglected again the contribution coming from the initial condition. Equation (6) defines a recursive succession, with limit point equal to $c(t)/(1-d)$, the periodicity of $c(t)$ concludes the proof. This result can also be framed in the general context of Floquet theory [77]. The Floquet theorem states that a fundamental matrix solution of a first order differential equation with periodically driven coefficients is quasi-periodical, namely can be written as $y(t)=P(t)e^{Mt}$ where $P(t)$ is a periodic matrix function (with the same period of the coefficients) and $e^{Mt}$ is the so called monodromy matrix. The real parts of the eigenvalues of $M$ are responsible of the asymptotic behavior of the solutions and are known as Lyapunov exponents, a stable cyclic solution is characterized by their negativity. In the case of Eq. (4) of the main text, our calculations reveal that the monodromy matrix is given by the constant $d$, the sign of the Lyapunov exponent is given by $\log d<0$, confirming our predictions about the stability.

In the above paragraph we showed that the asymptotic solution of Eq. (4) of the main text is periodic with the same period of the external driving $\epsilon(t)$. Notice that in the equilibrium scenario the previous result is trivial, since the population instantly relaxes to the Gibbs state that is a monotonic function of the control parameter $\epsilon$. In our case we can establish only that $p(t)$ and $\epsilon(t)$ share the same period, although finding the proper functional relation between the two is absolutely non trivial (cfr. for example [43]). However we don’t need any additional information to prove that any cycle that is not infinitesimal, namely a square wave protocol in which the controls jump at a time much faster that the typical dynamical scale $\Gamma$, cannot achieve the maximum power. The proof is outlined in the following: since $\epsilon(t)$ and $p(t)$ share the same periodicity, a cycle can be represented in the $(p,\epsilon)$ plane as a closed curve. Let us suppose that the optimal cycle $\mathcal{T}$ is not infinitesimal, for example as in Fig. 4. Thus, it is possible to perform an instantaneous quench, for example, in the middle (where the probability is halfway between the minimum and maximum value), and divide the transformation in two smaller sub-cycles $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$ (cfr. Fig 4). Since the quenches are instantaneous, they don’t contribute to the heat exchanged and to the time duration of the process. Furthermore, performing the two sub-cycles in series effectively builds a transformation with the same average power of the original cycle, a property that in symbols we can exemplify as $P(\mathcal{T})=P(\mathcal{T}_{1}\circ\mathcal{T}_{2})$. Simple calculations reveal that the power of the single sub-cycles cannot be both greater or smaller than the power of the original one, thus we are left with two possibilities, $P(\mathcal{T}_{1})\leq P(\mathcal{T}_{1}\circ\mathcal{T}_{2})\leq P(\mathcal{T}_{2})$ or $P(\mathcal{T}_{2})\leq P(\mathcal{T}_{1}\circ\mathcal{T}_{2})\leq P(\mathcal{T}_{1})$. In both cases the original cycle is sub-optimal, that is absurd, unless $P(\mathcal{T}_{1})=P(\mathcal{T}_{1}\circ\mathcal{T}_{2})=P(\mathcal{T}_{2})$ but even in this case we can choose one of the two sub-cycles still preserving optimality.

The previous argument shows that the only candidates for power maximization are those cycles that cannot be divided with a quench as done in the above proof, thus being infinitesimal. Notice that the previous proof strongly relies on the possibility of performing effectively instantaneous quenches, a characteristic that is better analyzed in the next appendix. At last, by using Pontryagin’s minimum principle, it can be shown that if coupling constants $\lambda_{\mathrm{H}}(t)$ and $\lambda_{\mathrm{C}}(t)$ fulfill a “trade-off relation” (i.e. if one increases, the other one decreases), then the optimal cycle will have $\lambda_{\mathrm{C}}(t)=0$ and $\lambda_{\mathrm{H}}(t)=1$, or $\lambda_{\mathrm{C}}(t)=0$ and $\lambda_{\mathrm{H}}(t)=1$ at all times [43]. This implies that the coupling to the baths must be switched during the quenches of the infinitesimal Otto-cycle.

Appendix B. Maximum Power Formula and Finite-Time Corrections

In this appendix we prove Eq. (5) of the main text and discuss the finite-time corrections.

As as shown in the previous appendix, the optimal cycle must be an infinitesimal Otto cycle, so we consider a protocol (depicted in Fig. 1b of the main text) where $\epsilon(t)=\epsilon_{\mathrm{H}}$, $\lambda_{\mathrm{H}}=1$ and $\lambda_{\mathrm{C}}=0$ for $t\in[0,\tau_{\mathrm{H}}]$, and $\epsilon(t)=\epsilon_{\mathrm{C}}$, $\lambda_{\mathrm{H}}=0$ and $\lambda_{\mathrm{C}}=1$ for $t\in[\tau_{\mathrm{H}},\tau_{\mathrm{H}}+\tau_{\mathrm{C}}]$. The optimal cycle and corresponding power will then be found by taking the limit $dt=\tau_{\mathrm{H}}+\tau_{\mathrm{C}}\to 0$ and maximizing over the free parameters $\epsilon_{\mathrm{H}}$, $\epsilon_{\mathrm{C}}$ and $\tau_{\mathrm{H}}/\tau_{\mathrm{C}}$.

We proceed the following way: first we perform an exact calculation, for arbitrary $\tau_{\mathrm{H}}$ and $\tau_{\mathrm{C}}$, of the heat rates $\braket{J_{\mathrm{H}}}$, $\braket{J_{\mathrm{C}}}$, averaged over one period, flowing out of the hot and cold bath respectively. Then, in the limit $dt\to 0$, we find the ratio $\tau_{\mathrm{H}}/\tau_{\mathrm{C}}$ that maximizes the power and we find the corresponding expression of the maximum power, proving Eq. (5) of the main text and the optimal ratio condition

[TABLE]

The instantaneous heat currents can be written in terms of the probability $p(t)$ by plugging the solution of Eq. (4) of the main text into Eq. (3) of the main text. We denoted with $p_{\mathrm{H}}(t)$ and $p_{\mathrm{C}}(t)$ the solution of Eq. (4) of the main text respectively in the time intervals $\mathcal{I}_{\mathrm{H}}=[0,\tau_{\mathrm{H}}]$ and $\mathcal{I}_{\mathrm{C}}=[\tau_{\mathrm{H}},\tau_{\mathrm{H}}+\tau_{\mathrm{C}}]$. Since the control parameters (i.e. $\epsilon(t)$, $\lambda_{\mathrm{H}}(t)$ and $\lambda_{\mathrm{C}}(t)$) are constant in each interval, we have that

[TABLE]

where $H$ and $C$ are two constants and where, for ease of notation, we introduced the symbols $\Gamma_{\alpha}:=\Gamma_{\alpha}(\epsilon_{\alpha})$ and $p_{\mathrm{eq}}^{(\alpha)}:=p_{\mathrm{eq}}^{(\alpha)}(\epsilon_{\alpha})$ (for $\alpha=\mathrm{H},\mathrm{C}$). We determine the two constants $H$ and $C$ by imposing that the probability $p(t)$ is continuous in $t=\tau_{\mathrm{H}}$, i.e.

[TABLE]

and that $p(t)$ is periodic with period $\tau_{\mathrm{H}}+\tau_{\mathrm{C}}$, i.e.

[TABLE]

We impose periodic boundary conditions because, as discussed in the previous appendix, a periodic protocol produces a periodic $p(t)$ after an initial transient time, and we are indeed interested in the “asymptotic” regime. Equations (3) and (4) reduce to the following linear-algebra problem for the constants $H$ and $C$:

[TABLE]

with solution

[TABLE]

which, via Eq. (2), completely determine $p(t)$. By substituting Eq. (2) into Eq. (3) of the main text, we can write the averaged heat rates $\braket{J_{\mathrm{H}}}$ and $\braket{J_{\mathrm{C}}}$ as

[TABLE]

where we use the fact that $\epsilon(t)$ is constant in each $\mathcal{I}_{\alpha}$ and the fact that, since the two-level system is coupled to one bath at a time, $\dot{p}_{\alpha}=\dot{p}$ with $p(t)=p_{\mathrm{H}}(t)$ during $\mathcal{I}_{\mathrm{H}}$ and $p(t)=p_{\mathrm{C}}(t)$ during $\mathcal{I}_{\mathrm{C}}$. Using the expressions for $H$ and $C$ given in Eq. (6), we can rewrite Eq. (7) as

[TABLE]

We now impose that $dt=\tau_{\mathrm{H}}+\tau_{\mathrm{C}}$ by setting $\tau_{\mathrm{H}}=\theta dt$ and $\tau_{\mathrm{C}}=(1-\theta)dt$, for $\theta\in[0,1]$, in Eq. (8). Taking hence the infinitesimal cycle limit $dt\to 0$ we get

[TABLE]

The maximization over $\theta$ of the above expression yields the condition

[TABLE]

which, multiplying by $dt$ the numerator and the denominator of the left hand side of Eq. (10), proves Eq. (1). Solving hence Eq. (10) for $\theta$ and plugging the result into Eq. (9) yields

[TABLE]

which replaced into Eq. (1) of the main text, and maximizing with respect to the only two free parameters left, i.e. $\epsilon_{\mathrm{H}}$ and $\epsilon_{\mathrm{C}}$, allows us to derive Eq. (5) of the main text for all four thermal machine modes. An additional comment has to be made for the accelerator mode [A], that aims at maximizing the heat released into the cold bath while extracting heat from the hot bath. By definition, we must restrict the maximization in Eq. (5) of the main text to guarantee $\langle J_{\mathrm{H}}\rangle\geq 0$, e.g. by forcing $\mathcal{C}$ to be $(\epsilon_{\mathrm{H}}\geq 0\cap\beta_{\mathrm{C}}\epsilon_{\mathrm{C}}\geq\beta_{\mathrm{H}}\epsilon_{\mathrm{H}})\cup(\epsilon_{\mathrm{H}}\leq 0\cap\beta_{\mathrm{C}}\epsilon_{\mathrm{C}}\leq\beta_{\mathrm{H}}\epsilon_{\mathrm{H}})$. On the other hand, the heater mode consists of heating a single reservoir whose interaction with the two-level system is described by a rate $\Gamma(\epsilon)$ and equilibrium probability $p_{\mathrm{eq}}(\epsilon)$. So in this case the maximization must be performed taking $\Gamma_{\alpha}(\epsilon)=\Gamma(\epsilon)$ and $p_{\mathrm{eq}}^{(\alpha)}(\epsilon)=p_{\mathrm{eq}}(\epsilon)$ (for $\alpha=$ H, C). If we also require that $\Gamma(\epsilon)=\Gamma(-\epsilon)$, which physically means that the rates do not distinguish which one of the two energy levels is the ground and excited state, we find that Eq. (5) can be simplified to

[TABLE]

and the corresponding optimal cycle is given by an Otto cycle where $\tau_{\mathrm{H}}=\tau_{\mathrm{C}}$ and the value $\epsilon$ that maximizes Eq. (12) determines $\epsilon_{\mathrm{H}}=-\epsilon_{\mathrm{C}}=\epsilon$. Thus the optimal cycle in the heater case corresponds to attempting continuous population inversions.

B.1 Finite-Time Corrections part one

Setting $\tau_{\mathrm{H}}=\theta dt$ and $\tau_{\mathrm{C}}=(1-\theta)dt$ in Eq. (8), and plugging in the expression of $\theta$ that satisfies Eq. (10), we find that the average heat rate for an arbitrary period $dt$ is given by

[TABLE]

where $\widetilde{\Gamma}_{\alpha}=(\widetilde{\Gamma}\Gamma_{\alpha})^{1/2}$ and $\widetilde{\Gamma}=\Gamma_{\mathrm{H}}\Gamma_{\mathrm{C}}/(\sqrt{\Gamma_{\mathrm{H}}}+\sqrt{\Gamma_{\mathrm{C}}})^{2}$. Plugging this results into Eq. (1) of the main text and maximizing over $\epsilon_{\mathrm{H}}$ and $\epsilon_{\mathrm{C}}$ yields the expression

[TABLE]

which provides the finite time version of Eq. (5) of the main text. On one hand, as anticipated in the main text, by expanding Eq. (14) for small $dt$, we find the following quadratic correction

[TABLE]

On the other hand for $\widetilde{\Gamma}_{\mathrm{H}}dt,\widetilde{\Gamma}_{\mathrm{C}}dt\gg 1$, we get

[TABLE]

implying that a considerable fraction of $P^{(\max)}_{[\nu]}$ can be achieved even if the driving frequency is slower than the typical rate. Notice that Eq. (14) is a strictly decreasing function of $dt$; this is consistent with the fact that an infinitesimal cycle is indeed the optimal solution.

We conclude by observing that we can simplify Eq. (14) for the heater mode where a single reservoir is coupled to the two-level system. Under the hypothesis leading to Eq. (12), we find that

[TABLE]

where $\Gamma$ is computed in the value of $\epsilon$ that maximizes Eq. (12). Figure 1c of the main text, which is a plot of Eq. (17), shows that $P^{(\max)}_{\mathrm{[H]}}(dt)\approx P^{(\max)}_{\mathrm{[H]}}$ up to $dt\Gamma\approx 2$, while for $dt\Gamma=10\gg 1$, $P^{(\max)}_{\mathrm{[H]}}(dt)$ is only decreased of a factor two.

B.2 Finite-Time Corrections part two: the Quenches

Finite-time corrections to the power may not only arise from the finite duration of the isothermal transformations (i.e. from a finite value of $\tau_{\mathrm{C}}$ and $\tau_{\mathrm{H}}$), but also from a finite duration $\tau$ of the quenches, i.e. of the transformations during which $\epsilon$ changes between the two extremal values $\epsilon_{\mathrm{C}}$ and $\epsilon_{\mathrm{H}}$. We will thus assume that each quench is carried out in a time $\tau$. The aim of this appendix is to show how these effects could be accounted for, and to estimate the leading order corrections to the maximum power delivered by the heat engine due to this effect; analogous considerations hold also for the other machines. We will thus restrict ourself to the regime $\tau\ll dt\ll\gamma^{-1}$, where $dt=\tau_{\mathrm{C}}+\tau_{\mathrm{H}}$ and $\gamma$ is the characteristic rate of the system during the protocol. The first inequality states that the duration of the quenches is much smaller than the duration of the isothermal transformations. The second inequality implies that the finite-time corrections discussed in the previous subsection are neglected, since they have been previously discussed.

Using the results of App. A, we know that $p(t)$ has a limit cycle with the same period of $\epsilon(t)$. If we further assume that the protocol is much faster than $\gamma$, the probability tends to a fixed value $\bar{p}$ given by

[TABLE]

where $T=dt+2\tau$ is the total duration of the protocol, $\Gamma(s)=\lambda_{\mathrm{H}}(s)\Gamma_{\mathrm{H}}[\epsilon(s)]+\lambda_{\mathrm{C}}(s)\Gamma_{\mathrm{C}}[\epsilon(s)]$ and

[TABLE]

By using again the hypothesis that the protocol is much faster than $\gamma$, we can write the power of the heat engine, averaged over one period, as

[TABLE]

As in the ideal protocol (see Fig 1b of the main text), we will assume that during the two isothermal transformations we respectively have $\epsilon(s)=\epsilon_{\mathrm{H}}$, $\lambda_{\mathrm{H}}(s)=1$, $\lambda_{\mathrm{C}}(s)=0$ and $\epsilon(s)=\epsilon_{\mathrm{C}}$, $\lambda_{\mathrm{H}}(s)=0$, $\lambda_{\mathrm{C}}(s)=1$. This means that we are coupled to one bath at a time. Instead, during the quenches we assume that all three control parameters ($\epsilon(s)$, $\lambda_{\mathrm{H}}(s)$ and $\lambda_{\mathrm{C}}(s)$) vary linearly in time between the corresponding extremal values. We thus divide the integral in Eq. (20) in the four different transformations:

[TABLE]

where $(\dots)$ stands for $ds\,\epsilon(s)\Gamma(s)\left[p_{\mathrm{eq}}(s)-\bar{p}\right]$. In the regime we consider the power, up to leading order corrections in $\tau/dt$, can be written as

[TABLE]

where the first addend is obtained by means of a first order expansion in $\tau/dt$ of the denominator in the r.h.s. of Eq. (21).

We wish to compare Eq. (22) to the power $P_{\mathrm{[E]}}^{(\max)}$ achieved in the ideal protocol, so we will estimate the four terms $\mathcal{W}_{\mathrm{H}}$, $\mathcal{W}_{\mathrm{C}}$, $\mathcal{W}_{\mathrm{H}\to\mathrm{C}}$ and $\mathcal{W}_{\mathrm{C}\to\mathrm{H}}$. First, we notice that $\bar{p}$ depends on the whole protocol, see Eq. (18). We can thus write

[TABLE]

where $\bar{p}^{(0)}$ is the value of $\bar{p}$ in the ideal protocol, and $\delta\bar{p}^{(1)}$ the corrections due to the finite-time quenches. These two terms can be calculated simply by dividing the integrals in the definition of $\bar{p}$ as we did for $P_{\mathrm{[E]}}(\tau)$. It is easy to see that $\delta\bar{p}^{(1)}$ is of the order $\tau/dt$. Since $\mathcal{W}_{\mathrm{H}}$ and $\mathcal{W}_{\mathrm{C}}$ are linear functions of $\bar{p}$, and $\bar{p}=\bar{p}^{(0)}+\delta\bar{p}^{(1)}$, we have that, for $\alpha=\mathrm{H},\mathrm{C}$,

[TABLE]

where $\mathcal{W}_{\alpha}^{(0)}$ is the work extracted in the ideal protocol during the isothermal transformation, while $\mathcal{W}_{\alpha}^{(1)}$ represents the corrections due to the variation in population $\bar{p}$ induced by the finite-time quenches. We have that

[TABLE]

where the last term means that $\mathcal{W}_{\alpha}^{(1)}$ is of the order $\gamma\tau$.

Next, we need to estimate $\mathcal{W}_{\mathrm{H}\to\mathrm{C}}$ and $\mathcal{W}_{\mathrm{C}\to\mathrm{H}}$. By inspecting the definition, we see that

[TABLE]

where $\braket{\epsilon}$ is a characteristic value of the energy gap during the quench.

Now we can return to Eq. (22). Using Eqs. (24), (25) and (26), and noticing that the order of magnitude of $\gamma\braket{\epsilon}$ is the same as $P_{\mathrm{[E]}}^{(\max)}$, we find that all the corrections previously discussed are of the order $\gamma/dt$. We thus conclude that

[TABLE]

where the corrections must be negative by virtue of the theorem proved in App. A. The impact of finite time quenches is thus first order $\tau/dt$.

Appendix C. Efficiency at Maximum Power

For small temperature differences, i.e. for small values of $\eta_{\mathrm{c}}$, we can consider an expansion of the efficiency at maximum power of the kind

[TABLE]

In this appendix we prove that $a_{1}=1/2$, while for symmetric or constant rates we further have $a_{2}=1/8$. The maximum power of a heat engine (without constraints on the control parameters) can be written as [see Eq. (5) of the main text]

[TABLE]

where

[TABLE]

$x_{\alpha}=\epsilon_{\alpha}\beta_{\alpha}$ (for $\alpha$ = H, C), $f(x):=[1+\exp{(x)}]^{-1}$ and, expressing the $\Gamma_{\alpha}$ as a function of the gap $\epsilon$ and of the inverse temperature $\beta_{\alpha}$ of lead $\alpha$,

[TABLE]

In Eq. (4) we decide to express $\beta_{\mathrm{C}}$ as $\beta_{\mathrm{H}}/(1-\eta_{\mathrm{c}})$ because we are interested in performing an expansion in $\eta_{\mathrm{c}}$ around a single inverse temperature $\beta_{\mathrm{H}}$. Let $x^{*}_{\mathrm{H}}$ and $x^{*}_{\mathrm{C}}$ be respectively the values of $x_{\mathrm{H}}$ and $x_{\mathrm{C}}$ that maximize $P_{[\mathrm{E}]}(x_{\mathrm{H}},x_{\mathrm{C}})$. By inspecting Eqs. (3) and (4), we see that $x^{*}_{\alpha}$ is a function of $\eta_{\mathrm{c}}$ (and of $\beta_{\mathrm{H}}$ through $g$), so we can express $x^{*}_{\alpha}$ as a power series in $\eta_{\mathrm{c}}$:

[TABLE]

Both $x^{*}_{\mathrm{H}}$ and $x^{*}_{\mathrm{C}}$ have the same leading order term. This can be seen considering Eq. (3) at $\eta_{\mathrm{c}}=0$: $g(x_{\mathrm{H}},x_{\mathrm{C}};\,0)/\beta_{\mathrm{H}}\geq 0$, while $\left[x_{\mathrm{H}}-x_{\mathrm{C}}\right]\left[f(x_{\mathrm{H}})-f(x_{\mathrm{C}})\right]\leq 0$, so the maximum power is zero (at equal temperatures, the second law forbids the possibility of extracting work). Inspecting Eq. (3), it is easy to see that zero power at $\eta_{\mathrm{c}}=0$ implies $x_{\mathrm{H}}=x_{\mathrm{C}}$. Using Eq. (7) of the main text, we have that

[TABLE]

so plugging Eq. (5) into Eq. (6) and expressing $\eta(P_{[\mathrm{E}]}^{(\max)})$ as a power series in $\eta_{\mathrm{c}}$, we find that

[TABLE]

where

[TABLE]

Thus, the knowledge of $b_{1}$ and $b_{2}$ implies also the knowledge of $a_{1}$ and $a_{2}$. Also the maximum power $P_{[\mathrm{E}]}(x^{*}_{\mathrm{H}},x^{*}_{\mathrm{C}})$ can be written as a power series in $\eta_{\mathrm{c}}$ by plugging the expansion Eq. (5) into Eq. (3). This yields

[TABLE]

where the coefficients $P_{[\mathrm{E}]}^{(n)}$ are functions of $m_{i}$, $n_{i}$ (for $i=0,1,2,\dots$) and of $\beta_{\mathrm{H}}$. We now wish to determine $b_{1}$ and $b_{2}$ by maximizing $P_{[\mathrm{E}]}^{(n)}$, starting from the lowest orders. We find that $P_{[\mathrm{E}]}^{(0)}=P_{[\mathrm{E}]}^{(1)}=0$ and

[TABLE]

where we expressed $n_{1}$ in terms of $b_{1}$. The last fraction in Eq. (9) is positive, so $P_{[\mathrm{E}]}^{(2)}$ is maximized by choosing $b_{1}$ that maximizes the term in square brackets, and $m_{0}$ that maximizes the last fraction. The maximization of the first term yields $b_{1}=-1/2$, which readily implies [see Eq. (7)] $a_{1}=1/2$, as we wanted to prove. The maximization of the second term allows us to find the following implicit expression for $m_{0}$

[TABLE]

where $\partial_{x_{\alpha}}g(m_{0},m_{0};\,0)$ denotes the partial derivative of $g(x_{\mathrm{H}},x_{\mathrm{C}};\,\eta_{\mathrm{c}})$, respect to $x_{\alpha}$, calculated in $x_{\mathrm{H}}=x_{\mathrm{C}}=m_{0}$ and $\eta_{\mathrm{c}}=0$. In order to compute $b_{2}$, we must maximize also higher order terms of the power. It turns out that $P_{[\mathrm{E}]}^{(3)}$ only depends on $m_{0}$ if we impose that $b_{1}=-1/2$ and that $m_{0}$ satisfies Eq. (10). Thus, there is nothing to optimize, so we must analyze the next order. $P_{[\mathrm{E}]}^{(4)}$ is a function of $m_{0}$, $m_{1}$, $n_{1}$, $m_{2}$, $n_{2}$ and $\beta_{\mathrm{C}}$. We write $m_{1}$ in terms of $b_{2}$, which is the only coefficient that determines $a_{2}$. We further express $n_{1}$ in terms of $b_{1}$, and impose $b_{1}=-1/2$. At last, we write $g(m_{0},m_{0};\,0)$ in terms of its partial derivatives using Eq. (10). This leads to an expression of $P_{[\mathrm{E}]}^{(4)}$ as a function of $m_{0}$ (which is implicitly known), $b_{2}$, $m_{2}$, $n_{2}$ and $\beta_{\mathrm{H}}$. We maximize $P_{[\mathrm{E}]}^{(4)}$ by setting to zero both partial derivatives of $P_{[\mathrm{E}]}^{(4)}$ respect to $b_{2}$ and $m_{2}$. We thus find the following expression for $b_{2}$:

[TABLE]

where all partial derivatives of $g$ are computed in $x_{\mathrm{H}}=x_{\mathrm{C}}=m_{0}$ and $\eta_{\mathrm{c}}=0$. This is, in principle, a closed expression for $b_{2}$, thus for $a_{2}$, since $m_{0}$ is defined in Eq. (10), and Eq. (11) only depends on $m_{0}$. Eq. (11) shows that in general $b_{2}$, thus $a_{2}$, will depend on the specific rates. However, if $\partial_{x_{\mathrm{H}}}g=\partial_{x_{\mathrm{C}}}g$, the first term in Eq. (11) vanishes, while the second one reduces to a number, yielding $b_{2}=-3/4$. Indeed, plugging this value of $b_{2}$ into Eq. (7) yields precisely $a_{2}=1/8$. We conclude the proof by noticing that if the rates are symmetric, i.e. $\Gamma_{\mathrm{H}}(\epsilon,\beta)=\Gamma_{\mathrm{C}}(\epsilon,\beta)$, $g(x_{\mathrm{H}},x_{\mathrm{C}};\,0)$ is a symmetric function upon exchange of $x_{\mathrm{H}}$ and $x_{\mathrm{C}}$. This implies that $\partial_{x_{\mathrm{H}}}g(m_{0},m_{0};\,0)=\partial_{x_{\mathrm{C}}}g(m_{0},m_{0};\,0)$, so $a_{2}=1/8$. At last, if the rates are constants, also $g(x_{\mathrm{H}},x_{\mathrm{C}};\,\eta_{\mathrm{c}})$ is constant, trivially satisfying $\partial_{x_{\mathrm{H}}}g=\partial_{x_{\mathrm{C}}}g=0$.

Appendix D. COP at Maximum Power

In this appendix we prove Eqs. (9) and (10) of the main text and we derive the scaling of the COP at maximum power for large values of the maximum gap $\Delta$ given by $\mathrm{C}_{\mathrm{op}}(P^{\max}_{[\mathrm{R}]})\propto 1/(\beta_{\mathrm{C}}\Delta)$. The COP at maximum power can be written as [see Eq. (8) of the main text]

[TABLE]

where $\epsilon^{*}_{\mathrm{H}}$ and $\epsilon^{*}_{\mathrm{C}}$ are respectively the values of $\epsilon_{\mathrm{H}}$ and $\epsilon_{\mathrm{C}}$ that maximize [see Eq. (5) of the main text]

[TABLE]

where $f(x):=[1+\exp{(x)}]^{-1}$ and

[TABLE]

We first prove that the COP at maximum power takes the universal form of Eq. (9) of the main text if the rates depend on the energy and on the temperature only through $\beta\epsilon$, i.e. $\Gamma_{\alpha}(\epsilon)=\Gamma_{\alpha}(\beta_{\alpha}\epsilon_{\alpha})$. We rewrite Eq. (1) as a function of $x^{*}_{\alpha}=\beta_{\alpha}\epsilon^{*}_{\alpha}$ (for $\alpha=$ H, C):

[TABLE]

where $\mathrm{C}_{\mathrm{op}}^{(\mathrm{c})}$ is the Carnot COP for a refrigerator (see main text). We can determine $x^{*}_{\alpha}$ by maximizing

[TABLE]

Crucially, given our hypothesis on the rates, there is no explicit dependence on the temperatures in Eq. (5) (except for the prefactor $1/\beta_{\mathrm{C}}$), so the maximization of $P_{[\mathrm{R}]}(x_{\mathrm{H}},x_{\mathrm{C}})$ will simply yield two values of $x^{*}_{\mathrm{H}}$ and $x^{*}_{\mathrm{C}}$ that do not depend on the temperatures. Thus, for all bath temperatures the COP at maximum power will be given by Eq. (4), where $x^{*}_{\mathrm{H}}$ and $x^{*}_{\mathrm{C}}$ are two fixed values. The ratio $x^{*}_{\mathrm{H}}/x^{*}_{\mathrm{C}}$ will depend on the specific rates we consider. By imposing in Eq. (4) that the COP at maximum power of the system for $\beta_{\mathrm{H}}=\beta_{\mathrm{C}}$ (i.e. for $\mathrm{C}_{\mathrm{op}}^{(\mathrm{c})}\to\infty$) is $\mathrm{C}_{\mathrm{op}}^{(0)}$, we can eliminate the ratio $x^{*}_{\mathrm{H}}/x^{*}_{\mathrm{C}}$ in favor of $\mathrm{C}_{\mathrm{op}}^{(0)}$, concluding the proof of Eq. (9) of the main text.

We now prove Eq. (10) of the main text. Since Eq. (2) remains unchanged by sending both $\epsilon_{\mathrm{H}}\to-\epsilon_{\mathrm{H}}$ and $\epsilon_{\mathrm{C}}\to-\epsilon_{\mathrm{C}}$, we can assume without loss of generality that $\epsilon_{\mathrm{C}}\geq 0$ (this is a general property which applies independently of the specific choice of bath models). Furthermore, we must ensure that the system is acting as a refrigerator by imposing $P_{[\mathrm{R}]}(\epsilon_{\mathrm{H}},\epsilon_{\mathrm{C}})\geq 0$. This implies that $f(\beta_{\mathrm{H}}\epsilon_{\mathrm{H}})\leq f(\beta_{\mathrm{C}}\epsilon_{\mathrm{C}})$, thus

[TABLE]

We now show that in the models described by Eq. (6) of the main text, the partial derivative of $P_{[\mathrm{R}]}(\epsilon_{\mathrm{H}},\epsilon_{\mathrm{C}})$ respect to $\epsilon_{\mathrm{H}}$ is non negative for all values of $\epsilon_{\mathrm{H}}$ and $\epsilon_{\mathrm{C}}$ satisfying Eq. (6), which implies that $\epsilon^{*}_{\mathrm{H}}\to+\infty$. Using Eq. (6), the condition $\partial P_{[\mathrm{R}]}(\epsilon_{\mathrm{H}},\epsilon_{\mathrm{C}})/\partial\epsilon_{\mathrm{H}}\geq 0$ can be written as

[TABLE]

Since $\partial\ln{g(\epsilon_{\mathrm{H}},\epsilon_{\mathrm{C}})}/\partial\epsilon_{\mathrm{H}}$ has the same sign as $d\Gamma_{\mathrm{H}}(\epsilon_{\mathrm{H}})/d\epsilon_{\mathrm{H}}$, and since the r.h.s. of Eq. (7) is strictly negative, Eq. (7) is certainly satisfied whenever $\Gamma_{\mathrm{H}}(\epsilon_{\mathrm{H}})$ is a growing function. This proves that $\epsilon^{*}_{\mathrm{H}}\to+\infty$ when the baths are described by the $\mathrm{F}_{n}$ model [see Eq. (6) of the main text] even when the two baths have different powers $n$. The $\mathrm{B}_{n}$ model is more tricky to analyze since the rates are decreasing functions around the origin. Nonetheless, using Eq. (6) it is possible to show that Eq. (7) is satisfied also in the $\mathrm{B}_{n}$ model by plugging $\Gamma^{(\mathrm{B}_{n})}_{\alpha}(\epsilon)$ [see Eq. (6) of the main text] into Eq. (7). This result holds also when the two baths have different powers $n$.

We now know that $\epsilon^{*}_{\mathrm{H}}\to+\infty$ in the $\mathrm{F}_{n}$ and $\mathrm{B}_{n}$ models. Since both $\Gamma^{(\mathrm{F}_{n})}_{\mathrm{H}}(\epsilon)$ and $\Gamma^{(\mathrm{B}_{n})}_{\mathrm{H}}(\epsilon)$ diverge for $n>0$ when $\epsilon^{*}_{\mathrm{H}}\to+\infty$, we have that

[TABLE]

where, as before, $x_{\mathrm{C}}=\beta_{\mathrm{C}}\epsilon_{\mathrm{C}}$ and $h(x):=1$ for the $\mathrm{F}_{n}$ model and $h(x):=\coth{x/2}$ for the $\mathrm{B}_{n}$ model [see Eq. (6) of the main text]. Thus, using $x_{\mathrm{C}}$ instead of $\epsilon_{\mathrm{C}}$, and noting that $f(\epsilon_{\mathrm{H}}\beta_{\mathrm{H}})$ vanishes for $\epsilon_{\mathrm{H}}\to+\infty$, we can write $P_{[\mathrm{R}]}(+\infty,\epsilon_{\mathrm{C}})$ [see Eq. (2)] as

[TABLE]

Equation (9) is non-negarive for all values of $x_{\mathrm{C}}$ and it vanishes in $x_{\mathrm{C}}=0$ and $x_{\mathrm{C}}\to+\infty$ thanks to the exponential decrease of $f(x_{\mathrm{C}})$ for large values of $x_{\mathrm{C}}$. Therefore, Eq. (9) will be maximum for the finite value $x^{*}_{\mathrm{C}}$ that maximizes $x_{\mathrm{C}}^{n+1}h(x_{\mathrm{C}})f(x_{\mathrm{C}})$, and plugging $x^{*}_{\mathrm{C}}$ into Eq. (9) yields the first relation in Eq. (10) of the main text, where $c_{n}=(x^{*}_{\mathrm{C}})^{n+1}h(x_{\mathrm{C}}^{*})f(x^{*}_{\mathrm{C}})$. For $n=0$, we separately analyze the $\mathrm{F}_{0}$ and $\mathrm{B}_{0}$ models. In the $\mathrm{F}_{0}$ model, $g(\epsilon_{\mathrm{H}},\epsilon_{\mathrm{C}})=k_{\mathrm{H}}k_{\mathrm{C}}/(\sqrt{k_{\mathrm{H}}}+\sqrt{k_{\mathrm{C}}})^{2}$, so $P_{[\mathrm{R}]}(+\infty,\epsilon_{\mathrm{C}})$ can be written as

[TABLE]

where $r:=k_{\mathrm{H}}/k_{\mathrm{C}}$. Using the same argument as before, Eq. (10) implies a finite value of $x^{*}_{\mathrm{C}}$ which arises from the maximization of $x_{\mathrm{C}}f(x_{\mathrm{C}})$. We thus proved the first relation in Eq. (10) of the main text for the $\mathrm{F}_{0}$ model, where $c_{0}=r/(\sqrt{r}+1)^{2}\,x^{*}_{\mathrm{C}}f(x^{*}_{\mathrm{C}})$. At last, in the $\mathrm{B}_{0}$ model $g(\epsilon_{\mathrm{H}}\to+\infty,\epsilon_{\mathrm{C}})=k_{\mathrm{H}}k_{\mathrm{C}}\coth{(x_{\mathrm{C}}/2)}/[\sqrt{k_{\mathrm{H}}}+\sqrt{k_{\mathrm{C}}\coth{(x_{\mathrm{C}}/2)}}]^{2}$. Thus, $P_{[\mathrm{R}]}(+\infty,\epsilon_{\mathrm{C}})$ can be written as

[TABLE]

Again, $x^{*}_{\mathrm{C}}$ is a finite value which can be found by maximizing $r\coth{(x_{\mathrm{C}}/2)}/(\sqrt{r}+\sqrt{\coth{(x_{\mathrm{C}}/2)}})^{2}\,x_{\mathrm{C}}f(x_{\mathrm{C}})$. Only in this case, $x^{*}_{\mathrm{C}}$ depends on the ratio $r$. We thus proved the first relation in Eq. (10) of the main text for the $\mathrm{B}_{0}$ model, where $c_{0}=r\coth{(x^{*}_{\mathrm{C}}/2)}/(\sqrt{r}+\sqrt{\coth{(x^{*}_{\mathrm{C}}/2)}})^{2}\,x^{*}_{\mathrm{C}}f(x^{*}_{\mathrm{C}})$.

The second relation in Eq. (10) of the main text stems from the fact that in all models $\epsilon_{\mathrm{H}}^{*}\to+\infty$ while $\epsilon_{\mathrm{C}}^{*}$ is finite. Thus, Eq. (1) implies that the $\mathrm{C}_{\mathrm{op}}(P^{\max}_{[\mathrm{R}]})$ vanishes. At last we want to roughly estimate the behavior of $\mathrm{C}_{\mathrm{op}}(P^{\max}_{[\mathrm{R}]})$ in the presence of a large yet finite constraint on the maximum gap: $|\epsilon(t)|\leq\Delta$. Since $\epsilon_{\mathrm{H}}$ would diverge if there was no constraint, we can assume that, in the presence of $\Delta$, $\epsilon_{\mathrm{H}}^{*}=\Delta$. On the other hand, $\epsilon^{*}_{\mathrm{C}}$ is a finite quantity (which is given by $\epsilon^{*}_{\mathrm{C}}=x^{*}_{\mathrm{C}}/\beta_{\mathrm{C}}$ in the unconstrained case), so if we assume that $\Delta\gg\epsilon^{*}_{\mathrm{C}}$, from Eq. (1) we have that

[TABLE]