Optimal cycles for low-dissipation heat engines

TL;DR

This paper analyzes how to optimize finite-time Carnot engines with small dissipation, showing that small cycles around an optimal point maximize power, which scales with heat capacity, enabling efficient many-body engines.

Contribution

It introduces a method to identify optimal cycle strategies for low-dissipation engines and demonstrates the potential for scalable many-body Carnot engines with finite power.

Findings

Optimal strategy involves small cycles around a specific working point.

Power output is proportional to the heat capacity of the working substance.

Many-body engines can achieve maximum efficiency at finite power per constituent.

Abstract

We consider the optimization of a finite-time Carnot engine characterized by small dissipations. We bound the power with a simple inequality and show that the optimal strategy is to perform small cycles around a given working point, which can be thus chosen optimally. Remarkably, this optimal point is independent of the figure of merit combining power and efficiency that is being maximized. Furthermore, for a general class of dynamics the power output becomes proportional to the heat capacity of the working substance. Since the heat capacity can scale supra-extensively with the number of constituents of the engine, this enables us to design optimal many-body Carnot engines reaching maximum efficiency at finite power per constituent in the thermodynamic limit.

Tables

Table 1. Table 1: Guide table for the Supplemental Material.

Section I	Description of a quasistatic, reversible Carnot cycle for a general quantum system.
Section II	Introduction of finite-time Carnot cycle and thermodynamic length.
IIA	Thermodynamic metric for exponential relaxations of the observables.
IIB	Thermodynamic metric from microscopical models.
IIB1	Example: Qubit with bosonic baths.
IIB2	Example: Harmonic oscillator with bosonic baths.
IIB3	Example: 3-level system with detailed balance evolution.
Section III	Generalizations of cycle optimization to general metrics and protocols.
IIIA	Generalization for asymmetric dissipation and symmetric protocols.
IIIA1	Optimization for systems with point-dependent thermalization timescale.
IIIA2	Optimization for general metrics.
IIIA3	Solution of the optimization for the models presented in Section IIB1-2-3.
IIIB	Bounds on completely asymmetric protocols.
Section IV	Maximization of the heat capacity for systems with different degree of control.
Section V	Time-tuning optimization of a low-dissipation Carnot engine.
VA	Critical scalings of power, efficiency and fluctuations of optimal protocols.
VB	Comparison with Otto cycle of Ref.Holubec and Ryabov (2017)
Section VI	Explicit analytical control protocol for a $D-1$ degenerate model.

Equations

$$G(s{\tau }_c)=\sum _j{\lambda }_j(s)X_j,s\in [0,1]$$

$$Q_x=T_x\Delta S_x-W_x^{\mathrm{diss}},x=(h,c)$$

$$W_x^{\mathrm{diss}}=\frac{T_x}{{\tau }_x}{\int }_0^1\sum _{ij}{\dot{\lambda }}_i{m}_{ij}{\dot{\lambda }}_{j}ds\equiv T_x\frac{{\Sigma }_x}{{\tau }_x},$$

$$m_{ij}={\tau }_{\mathrm{eq}}\frac{{\partial }^2\ln \mathcal{Z}}{\partial {\lambda }_i\partial {\lambda }_j}.$$

$$P^{(\max )}=\frac{(\Delta S{)}^2}{4\Sigma }(\sqrt{T_h}-\sqrt{T_c}{)}^2$$

$$\eta =\gamma {\eta }_C,\gamma \le 1,$$

$$P_{\gamma }^{(\max )}=\frac{(\Delta S{)}^2}{4\Sigma }\frac{(T_h-T_c{)}^2\gamma (1-\gamma )}{\gamma T_c+(1-\gamma )T_h},$$

$$\Sigma ={\tau }_{\mathrm{eq}}{\int }_0^{1}ds{\mathrm{cov}}_{\omega }(\dot{G},\dot{G})$$

$$\Delta S=-{\int }_0^{1}ds{\mathrm{cov}}_{\omega }(G,\dot{G}).$$

$$\frac{(\Delta S{)}^2}{\Sigma }\le \frac{1}{{\tau }_{\mathrm{eq}}}{\int }_0^{1}ds{\mathrm{cov}}_{\omega }(G,G)\equiv \frac{1}{{\tau }_{\mathrm{eq}}}{\int }_0^{1}ds\mathcal{C}$$

$$\mathcal{C}=-{\beta }^2\frac{{\partial }^2\ln \mathcal{Z}}{\partial {\beta }^2}=\mathrm{Tr}(\omega G^2)-\mathrm{Tr}(\omega G{)}^2$$

$$G({\tau }_c)=G(0)(1+ϵ)ϵ\ll 1,$$

$$G({\tau }_c+{\tau }_h)=G(0),$$

$$P_{\gamma }^{(\max )}=\frac{\mathcal{C}}{4{\tau }_{\mathrm{eq}}}\frac{(T_h-T_c{)}^2\gamma (1-\gamma )}{\gamma T_c+(1-\gamma )T_h}.$$

$${\tau }_{\mathrm{eq}}/\tau \sim {ϵ}^{-1}(T_h/T_c-1)(1-\gamma ).$$

$$P_{\gamma }^{(\max )}=\frac{c_0(T_c-T_h{)}^2}{4{\tau }_0{T}_c}{N}^{1+a-b-\xi }$$

$${\tau }_c=\frac{2ϵT_c}{T_h-T_c}{\tau }_0{N}^{b+\xi },{\tau }_h={\tau }_c,$$

$$W=(T_h-T_c)c_0ϵN^{1+a},{\sigma }_w^2=2(T_h-T_c)\frac{W}{ϵ}.$$

$$\alpha -z\nu -1>0.$$

$$U=\text{Tr}[\rho H].$$

$$dW=\text{Tr}[\rho dH],$$

$$dQ=\text{Tr}[d\rho H].$$

$${\omega }_{\beta }(\vec{l})=e^{-\beta H}/{\mathcal{Z}}_{\beta }(H)$$

$$G_{t^{\prime }}:=\beta H_{t^{\prime }\tau }$$

$${\eta }_{Carnot}=\frac{-\Delta W^{(0)}}{Q_{abs}^{(0)}}=1+\frac{Q_c^{(0)}}{Q_h^{(0)}}=1-\frac{T_c}{T_h},$$

$$\frac{Q_c^{(0)}}{T_c}={\int }_X^Y\text{Tr}[d\omega G]=-{\int }_Y^X\text{Tr}[d\omega G]=-\frac{Q_h^{(0)}}{T_h}.$$

$$\frac{d}{dt}\rho =\mathcal{L}[\rho ]\text{s.t.}\mathcal{L}[{\omega }_{\beta }(H(t))]=0.$$

$${\rho }_s={\omega }_s+\frac{1}{\tau }{\mathcal{L}}_s^{-1}[{\dot{\omega }}_s]+...$$

$$Q=T\Delta S-\frac{T}{\tau }\Sigma$$

$$\Sigma =\beta {\int }_0^{1}ds\text{Tr}({\mathcal{L}}_s^{-1}[{\dot{\omega }}_s]{\dot{H}}_s)$$

$$\dot{\omega }=-{\int }_0^{1}dx{\omega }^{1-x}(\dot{G}-\mathrm{Tr}(\omega \dot{G})){\omega }^x\equiv {\mathfrak{J}}_{\omega }[\dot{G}],$$

$$\Delta S=-{\int }_0^{1}ds\text{Tr}[{\mathfrak{J}}_{{\omega }_s}[{\dot{G}}_s]\ln {\omega }_s]={\int }_0^{1}ds\text{Tr}[G_s{\mathfrak{J}}_{{\omega }_s}[{\dot{G}}_s]],$$

Full text

The Carnot engine has a pivotal role in thermodynamics, both from a fundamental and applied perspective, being the reference point for other engines in terms of efficiency Carnot (1824); Callen (1998). It is thus of paramount importance to understand its limits and strategies for its best utilization. In this article, we consider the optimization of a finite-time Carnot cycle within the so called low-dissipation (LD) regime Esposito et al. (2010a); Guo et al. (2013); den Broeck (2013); Hernández et al. (2015); Holubec and Ryabov (2015, 2016); Ma et al. (2018a); de Tomás et al. (2012); de Tomas et al. (2013); Wang et al. (2012); Hu et al. (2013); Cavina et al. (2017), where the driving of the control parameters is slow but finite. Previous studies of Carnot engines in the LD regime have considered bounds on the reachable efficiencies Esposito et al. (2010a), tradeoffs between efficiency and power Holubec and Ryabov (2015, 2016); Shiraishi et al. (2016); Ma et al. (2018a), the coefficient of performance of refrigerators Wang et al. (2012); Hu et al. (2013), the impact of the spectral density of the thermal baths Cavina et al. (2017), and other thermodynamic figures of merit de Tomás et al. (2012); de Tomas et al. (2013). Despite this remarkable progress, the following crucial question has remained unaddressed: given a certain level of control on the working substance (e.g. some parameters of the Hamiltonian, or some macroscopic variables such as volume or pressure), what is the optimal cyclic modulation of the control parameters to maximise the power output (or, more generally, any figure of merit involving power and efficiency Holubec and Ryabov (2015, 2016); Shiraishi et al. (2016); Ma et al. (2018a)) of a finite-time Carnot engine? Such an optimal cycle has been designed for a single-qubit engine in Esposito et al. (2010b); Abiuso and Giovannetti (2019), but a general understanding is lacking.

Using recent insights on a geometrical approach to thermodynamics Weinhold (1975); Ruppeiner (1979); Salamon and Berry (1983); Schlögl (1985); Nulton et al. (1985); Crooks (2007); Sivak and Crooks (2012); Bonança and Deffner (2014); Zulkowski and DeWeese (2015); Mandal and Jarzynski (2016); Scandi and Perarnau-Llobet (2019); Miller et al. (2019), we show that, given any reasonable figure of merit involving power and efficiency, the optimal control strategy is to perform infinitesimal cycles around a fixed point. Furthermore, when the thermalization of the relevant quantities can be described by a single time-scale $\tau_{\rm eq}$ (see details below), the optimal power output becomes proportional to $\mathcal{C}/\tau_{\rm eq}$, where $\mathcal{C}$ is the heat capacity of the working substance (WS). Hence, the optimization of the heat engine cycle becomes intimately related to the maximization of $\mathcal{C}$ of the WS (interestingly, maximizing $\mathcal{C}$ is also crucial in thermometry Correa et al. (2015); Potts et al. (2019); Hovhannisyan and Correa (2018); Mehboudi et al. (2019)).

We then use these insights to design many-body heat engines that can operate at Carnot efficiency with finite power per constituent of the WS through a supraextensive scaling of $\mathcal{C}/\tau_{\rm eq}$ (e.g. in a phase transition), in the spirit of Allahverdyan et al. (2013); Campisi and Fazio (2016) (see also Campisi and Fazio (2016); Holubec and Ryabov (2017); Ma et al. (2017); Chand and Biswas (2018)). Despite differences w.r.t. previous proposals, which were based on Otto engines Allahverdyan et al. (2013); Campisi and Fazio (2016), we find the same asymptotic scalings for performance, while proving by construction their optimality in the slow-driving regime. Other recent proposals towards the possibility of reaching Carnot efficiency at finite power have been developed in  Mintchev et al. (2013); Campisi and Fazio (2016); Allahverdyan et al. (2013); Holubec and Ryabov (2017); Polettini et al. (2015); Koning and Indekeu (2016); Polettini and Esposito (2017); Johnson (2018); Holubec and Ryabov (2018); Erdman et al. (2018) (see also Shiraishi (2015); Shiraishi et al. (2016); Shiraishi and Tajima (2017) for no-go results 111Note that these results are not in contradiction with the present analysis, as we are concerned with reaching Carnot efficiency asymptotically.).

Finite-time Carnot cycle. We consider a quantum working substance (WS) with a driven Hamiltonian $H(t)$, which interacts alternatively with a cold ($B_{c}$) and a hot ($B_{h}$) heat bath at temperature $T_{c}$ and $T_{h}$, respectively (the results presented in this article can be naturally extended to classical systems). The Carnot cycle consists of four steps: (i) While being coupled to $B_{c}$, $H(t)$ is modified continuously from $H(0)=H^{(X)}$ to $H(\tau_{c}^{-})=H^{(Y)}$ for a time $\tau_{c}$. (ii) With the system isolated from the reservoirs, an adiabatic process is performed taking $H^{(Y)}\rightarrow H^{(Y)}T_{h}/T_{c}$. During this process $H(t)$ satisfies $[H(t),H(t^{\prime})]=0$ $\forall t,t^{\prime}$ and commutes with the initial state (a Gibbs state w.r.t $H(\tau_{c}^{-})$), so it is possible to perform it arbitrarily quickly without affecting the state (hence corresponding to Hamiltonian quenches). (iii) While being coupled to $B_{h}$, $H(t)$ is modified back from $H(\tau_{c}^{+})=H^{(Y)}T_{h}/T_{c}$ to $H(\tau_{c}+\tau_{h})=H^{(X)}T_{h}/T_{c}$ in a time $\tau_{h}$. (iv) Finally a second quench is performed to restore $H^{(X)}T_{h}/T_{c}\rightarrow H^{(X)}$, closing the cycle.

It is convenient to introduce the adimensional Hamiltonian $G(t):=\beta H(t)$, where $\beta=1/k_{B}T$ is the inverse temperature of the bath that the WS is coupled to (we shall set the Boltzmann constant $k_{B}$ equal to $k_{B}=1$). Note that the Carnot cycle becomes smooth with respect to $G(t)$, and in what follows, we take the driving protocol to be time-reversal symmetric, more precisely that $G(s\tau_{c})=G(\tau_{c}+\tau_{h}(1-s))$ with $s\in[0,1]$. This property is always satisfied by optimal heat engines in the LD regime as long as the two baths have the same spectral density Cavina et al. (2017), which is the subject of interest of this work. By expressing $G(t)$ as

[TABLE]

where $\lambda_{j}(s)$ are the control parameters and $X_{j}$ the conjugated forces, the cycle control can be characterised by $\tau_{c}$, $\tau_{h}$ and its shape $\vec{\lambda}(s)$ (notice that $\tau_{c}\neq\tau_{h}$ in general; the time-reversal symmetry is intended in the adimensional unit $s$). We will not write explicitly the dependence on $s$, which will be clear from the context, but will indicate with a dot the time derivative w.r.t. $s$, i.e. $\dot{G}\equiv\frac{\partial}{\partial s}G\equiv\tau_{x}\frac{d}{dt}G$, $x=(h,c)$.

We now assume the slow driving or low dissipation approximation, which is crucial in this work. That is, $\frac{d}{dt}{G}$ is finite but small, so that we can expand the relevant quantities and keep only leading terms (formally, the small parameter of the expansion is $\tau_{\rm eq}/\tau_{x}$, where $\tau_{\rm eq}$ is the relaxation time of the dynamics). In this regime, the state of the WS is always close to thermal equilibrium, and the heat exchanged during the isotherms (steps (i) and (iii)) can be divided as

[TABLE]

where $\Delta S_{x}$ is the reversible contribution obtained in the quasistatic limit $\tau_{x}\rightarrow\infty$, which is given by $\Delta S\equiv\Delta S_{h}=-\Delta S_{c}=S(\omega(\tau_{c}))-S(\omega(0))$ with $\omega(t)=e^{-G(t)}/\mathrm{Tr}(e^{-G(t)})$ and where $S(\rho)$ is the Von Neumann entropy. The irreversible term $W^{\rm diss}_{x}$ can be described in this regime by the so called thermodynamic length Schlögl (1985); Salamon and Berry (1983); Nulton et al. (1985); Crooks (2007); Sivak and Crooks (2012); Bonança and Deffner (2014); Scandi and Perarnau-Llobet (2019)

[TABLE]

where $m_{ij}$ is given, when the driven observables $\langle X_{j}\rangle$ relax to equilibrium with the same time-scale $\tau_{\rm eq}$, by Schlögl (1985); Salamon and Berry (1983); Nulton et al. (1985); Crooks (2007); Sivak and Crooks (2012); Bonança and Deffner (2014); Scandi and Perarnau-Llobet (2019):

[TABLE]

where $\mathcal{Z}=\mathrm{Tr}{(e^{-G})}$ is the partition function. Given the time-reversal symmetry of the driving protocol, it follows that $\Sigma_{h}=\Sigma_{c}$ (symmetric low-dissipation regime SLD). Importantly, while the results presented in the main text are based upon the standard thermodynamic metric (4), generalizations (including $\tau_{\rm eq}$ depending on the trajectory Sivak and Crooks (2012); Bonança and Deffner (2014), the possibility of having several relaxation time-scales, general Lindbladian dynamics Sivak and Crooks (2012); Scandi and Perarnau-Llobet (2019), and protocols in which $\Sigma_{h}\neq\Sigma_{c}$) are developed in the Supp. Mat. Abiuso and Perarnau-Llobet .

Optimisation of the cycle. We now optimise the power (and efficiency) of the Carnot engine over $\tau_{c}$, $\tau_{h}$ and its shape $\lambda_{j}(s)$, which are all the possible degrees of freedom. This enables us to obtain a fundamental upper bound on the power in the slow-driving regime, as well as the corresponding optimal control.

The work extracted during a cycle is given by $W=Q_{h}+Q_{c}$, and the total time is $\tau=\tau_{c}+\tau_{h}$. The power hence reads $P=(Q_{h}+Q_{c})/\tau$, and the efficiency $\eta=(Q_{h}+Q_{c})/Q_{h}$. By substituting the expressions (2)-(3) and appropriately setting $\tau_{c}$ and $\tau_{h}$, one can maximize the power of the engine ($\partial P/\partial\tau_{j}=0$) obtaining Schmiedl and Seifert (2007a); Esposito et al. (2010a)

[TABLE]

and the corresponding efficiency at maximum power (EMP) is given by the Curzon-Ahlborn EMP, $\eta_{CA}=1-\sqrt{T_{c}/T_{h}}$ Curzon and Ahlborn (1975). In the most general case one might seek, in order to not sacrifice completely the efficiency optimization over the power, to maximize a hybrid figure of merit Holubec and Ryabov (2015, 2016); Shiraishi et al. (2016); Ma et al. (2018a). The maximum efficiency for any given power output of the engine has been derived in Ma et al. (2018a) (see also Holubec and Ryabov (2015)). Analogously, we can express the best power for a given efficiency, fixed to be a fraction of the maximum one:

[TABLE]

where $\eta_{C}=1-T_{c}/T_{h}$ is the Carnot efficiency. In the SLD regime this leads to a maximum power (cf. Supp.Mat. Abiuso and Perarnau-Llobet ),

[TABLE]

obtained by setting $\tau_{c}=2\Sigma T_{c}/(\Delta S\left(T_{h}-T_{c}\right)(1-\gamma))$ and $\tau_{h}=\tau_{c}(T_{h}T_{c}^{-1}(1-\gamma)+\gamma)$. Essentially, by tuning $\gamma$ in $\tau_{c}$ and $\tau_{h}$, one can interpolate between a maximally powerful engine with power (5) at efficiency $\eta_{CA}$, and a Carnot engine with maximal efficiency and null power.

At this point, we note a crucial observation: after the optimisation of $P$ over $\tau_{c}$ and $\tau_{h}$, the remaining figure of merit is always ${(\Delta S)^{2}}/{\Sigma}$, independently of the value of $\gamma$. In fact, this is a property that can be argued to be general given any figure of merit combining power and efficiency Abiuso and Perarnau-Llobet . We now show how to maximize it by an opportune use of a Cauchy-Schwarz inequality. First, using the formula for the derivative of an exponential Hiai and Petz (2014): $\partial e^{-G}/\partial\lambda_{j}=-\int_{0}^{1}e^{-G}X_{j}e^{-(1-s)G}$, we can express again $\Sigma$ in (3) in the more compact form

[TABLE]

where ${\rm cov}_{\omega}(A,B)$ is the Kubo-Mori-Boguliobov inner product: ${\rm cov}_{\omega}(A,B)=\int_{0}^{1}ds\hskip 2.84526pt\mathrm{Tr}(A\omega^{1-s}(B-\mathrm{Tr}(\omega B))\omega^{s})$. Next, we note that the Von Neumann entropy $S=-\text{Tr}[\omega\ln\omega]$ satisfies: $\dot{S}=-\mathrm{Tr}(\dot{\omega}\ln\omega)=-{\rm cov}_{\omega}(\ln\omega,\dot{G})={\rm cov}_{\omega}(G,\dot{G})$ where we used again the formula $\dot{\omega}=-\int_{0}^{1}dx\hskip 2.84526pt\omega^{1-x}(\dot{G}-\mathrm{Tr}(\omega\dot{G}))\omega^{x}$, as well as $\mathrm{Tr}(\dot{\omega})=0$ and ${\rm cov}_{\omega}(A,\mathds{1})=0$. Hence, we can write the total change in entropy $\Delta S$ as:

[TABLE]

Crucially, both $\Delta S$ and $\Sigma$ can be expressed through an infinite-dimensional scalar product given by: $\langle A,B\rangle_{\omega}=\int_{0}^{1}ds\hskip 2.84526pt{\rm cov}_{\omega}(A,B)$, that depends on the path $\{\lambda_{j}(s)\}$. Using the Cauchy-Schwarz (C-S) inequality $|\langle A,B\rangle|^{2}\leq\langle A,A\rangle\langle B,B\rangle$, the ratio ${(\Delta S)^{2}}/{\Sigma}$ can then be bounded as:

[TABLE]

where $\mathcal{C}$ is the heat capacity of the WS,

[TABLE]

i.e. the variance of the adimensional Hamiltonian $G$ for its thermal state. The C-S inequality is saturated for $\dot{G}\propto G$, which means that creating quantum coherence cannot favour the power output in the slow driving regime, in agreement with Refs. Brandner et al. (2017); Ma et al. (2018b). More importantly, it shows that optimal thermodynamic protocols take the simple form $G(t)=\lambda(t)G(0)$.

We can further maximize (63) by noting that $\int_{0}^{1}ds\hskip 2.84526pt\mathcal{C}\leq\max_{s}\mathcal{C}$. To saturate this inequality in practice, one needs to consider cycles where the modulation of $G$ is small

[TABLE]

around an optimal point $G(0)$ where $\mathcal{C}/\tau_{\rm eq}$ is maximised (as long as $\epsilon$ is small enough so that $G(s\tau_{c})$ with $s\in[0,1]$ does not change substantially along the cycle, the precise form of $G(s\tau_{c})$ is not important; see the SM for examples of explicit cycles). In this case, in the limit $\epsilon\rightarrow 0$ the maximal power (7) for a given efficiency $\gamma\eta_{C}$ is given by

[TABLE]

where $\mathcal{C}$ is the heat capacity at $G(0)$. We stress that (13) has been obtained after maximising $P$ (for a fixed $\eta$) over all degrees of freedom: $\tau_{c}$, $\tau_{h}$ and the protocol $\{\lambda_{j}(s)\}$. This result shows a fundamental link between maximal power of a finite-time Carnot cycle and the heat capacity of the WS.

The simplicity of (13) can be contrasted to exact optimisations of heat engines Schmiedl and Seifert (2007b); Allahverdyan et al. (2010); Cavina et al. (2018); Erdman et al. (2018); Menczel et al. (2019), where the full solution easily becomes too complex or not even computable with the size of the WS; and with other geometric optimisations which require solving geodesic equations to design optimal paths in the parameter space Sivak and Crooks (2012); Zulkowski et al. (2013); Zulkowski and DeWeese (2015); Rotskoff and Crooks (2015); Rotskoff et al. (2017); Scandi and Perarnau-Llobet (2019); Miller et al. (2019). In our approach, from the point of view of optimization all that is left to do in (13) is to maximise $\mathcal{C}$ over the control parameters $\{\lambda_{j}\}$ to identify the optimal working point $G(0)$ in (12). In Fig. 6 explicit results are reported for the value of maximum $\mathcal{C}$, for different paradigmatic levels of control on the same system of $N$ qubits.

Crucially, this approach can be generalized to any metric $m_{ij}$ in (3) describing dissipation: in the SM Abiuso and Perarnau-Llobet we show that the optimal control problem is always reduced to infinitesimal cycles, and the optimal working point can be found by a scalar maximisation problem. We work out as well examples of standard microscopical dynamics in open quantum systems: a qubit, a qutrit, or an harmonic oscillator as a WS in contact with bosonic thermal baths with different spectral densities.

It is important to point out that (13) should be understood as a theoretical ultimate upper bound on power, obtained by taking $\epsilon\rightarrow 0$ in (12). In practice, any experimental or realistic protocol will have finite $\epsilon$. In this case, as long as $\mathcal{C}$ is sufficiently smooth along the thermodynamic cycle (12), the power output $P$ (given by the integrated $\mathcal{C}$ in (63)) will not change considerably, so that realistic cycles will provide a similar $P$ than the theoretically maximal one. In practice, keeping $\epsilon$ finite is also important to ensure the consistency of the slow-driving approximation $\tau_{\rm eq}/\tau\ll 1$, given that for the optimized protocols that lead to the power (7), one has $\tau\propto\epsilon$, more precisely

[TABLE]

Note that this can always be guaranteed for high efficiencies $\gamma\rightarrow 1$. From the same equation we note incidentally that engines whose maximum efficiency is constrained to be low ($T_{c}/T_{h}\lesssim 1$), i.e. arguably those engines that mostly need optimization in the high efficiency regime, show better convergence to the absolute bound.

Reaching Carnot efficiency at finite power. As an example of application of the previous results, we now use the designed optimal finite-time Carnot cycles (Eqs. (12) where $G(0)$ will depend on each model of interest) to explore the possibility of reaching Carnot efficiency $\eta_{C}$ at finite power in the macroscopic limit. We follow the approach put forward in Refs. Allahverdyan et al. (2013); Campisi and Fazio (2016): considering a $N$-particle WS, we aim at approaching Carnot efficiency in the macroscopic limit $N\rightarrow\infty$ without giving up power per constituent.

To reach Carnot efficiency, we need $\gamma=1$ in (6), and hence we take $1-\gamma=N^{-\xi},$ where $\xi>0$ can be chosen at will. On the other hand, the maximal power $P_{\gamma}^{({\rm max})}$ in (13) depends only on $\mathcal{C}$ and $\tau_{\rm eq}$; we then assume $\mathcal{C}=c_{0}N^{1+a}$ and $\tau_{\rm eq}=\tau_{0}N^{b}$, where the meaning of the different constants will later be described for each model of interest. Expanding the relevant quantities for $N\gg 1$, we obtain at leading order in $N$:

[TABLE]

where $\sigma_{w}^{2}=\langle w^{2}\rangle-(\langle w\rangle)^{2}$ is the variance of the work distribution, which measures the work fluctuations per cycle of the engine (see Supp.Mat. Abiuso and Perarnau-Llobet for details on the calculation). Let us now discuss two separate cases, inspired by Allahverdyan et al. (2013) and Campisi and Fazio (2016), respectively.

(a) Control on the engine and the engine-bath interaction. We first assume full control over the engine Hamiltonian with $2^{N}$ levels: that is, all levels can be modified at will by the experimentalist. While this is extremelly challenging in practice, it is useful to obtain fundamental upper bounds on the maximal power. The optimal Hamiltonian maximising $\mathcal{C}$ then consists in a ground level and a $2^{N}-1$ degenerate level (see Reeb and Wolf (2014); Correa et al. (2015)) which, as shown in Fig.6, leads to $\mathcal{C}\propto N^{2}$, i.e., $a=1$. This supralinear scaling is obtained in an increasingly smaller region of the parameter space, which requires taking $\epsilon\propto 1/N$ in (12) (see details in Supp. Mat. Abiuso and Perarnau-Llobet ), and this constraints from Eq. (14) also $1-\gamma$ to scale accordingly, i.e. $\xi=1$. Furthermore, it is possible to reach in realistic collisional scenarios $\tau_{\rm eq}\propto\sqrt{N}$ (i.e. $b=1/2$), or constant $\tau$ ($b=0$) if one is allowed to fine-tune the interaction between the WS and the baths Allahverdyan et al. (2013)222We note that because of the degeneracies in the Hamiltonian of the WS the steady state might not be unique (as the population in each degenerate energy state might depend on the initial state). However, this is is irrelevant as we are considering a driving $\dot{G}\propto G$, where the power output depends only on the total population in each degenerate energy level. , in agreement with Ref. Allahverdyan et al. (2013). In the Supp.Mat. Abiuso and Perarnau-Llobet , we solve exactly this proposal for a feasible driving protocol close to the optimal one.

(b) Engine working on a phase transition point. A promising platform to obtain supralinear scaling of power with realistic control is by choosing the engine to work in a phase transition point of the many-body WS. For a finite WS operating close to the critical point, finite size scaling theory tells us that $\mathcal{C}$ develops a peak of height $C\propto N^{1+\alpha/(\nu d)}$ and width $\delta\propto N^{-1/(d\nu)}$, while $\tau_{\rm eq}\propto N^{z/d}$ (here $\alpha$, $\nu$ and $z$ correspond to the specific heat, correlation length and dynamical critical exponents, while $d$ is the spatial dimension of the engine Fisher and Barber (1972); Suzuki (1977)). In order to exploit the critical scaling of the WS, we need to perform the cycle (12) where $\mathcal{C}$ becomes peaked, and hence $\epsilon\propto\delta\propto N^{-1/(d\nu)}$, which implies $\xi=1/(d\nu)$ from Eq. (14). Then, from (Optimal cycles for low-dissipation heat engines) with $a=\alpha/(\nu d)$ and $b=z/d$, supralinear scaling of $P_{\gamma}^{(\rm max)}$ is possible if

[TABLE]

This condition is the same found for the Otto cycle proposed in Campisi and Fazio (2016). Examples of physical systems where (16) is satisfied are also provided in Campisi and Fazio (2016), particularly in the presence of critical speed-ups of thermalisation ($z<0$ Boukari et al. (1990); Grams et al. (2014); Tavora et al. (2014)).

Besides efficiency and power, another crucial aspect of a heat engine is its reliability, i.e. the fluctuations in the output power. In fact, it has been recently pointed out in Holubec and Ryabov (2017) that the Otto-cycle of Campisi and Fazio (2016) suffers from macroscopic fluctuations in the thermodynamic limit. For the Carnot-cycle considered here, from (Optimal cycles for low-dissipation heat engines) the relative work fluctuations read $f_{w}=\sigma_{w}/W=\sqrt{2(T_{h}-T_{c})/\epsilon W}.$ First of all, in the case (a) where $a=1$ and $\epsilon\propto N^{-1}$ in (Optimal cycles for low-dissipation heat engines), one has $f_{w}\sim\mathcal{O}(1)$ in the macroscopic limit. A similar situation takes place for the critical heat engine (b) as one simultaneously has $\epsilon\propto N^{-1/(d\nu)}$ and $a=\alpha/(d\nu)$, and hence $f_{w}\propto N^{(-2+d\nu+\alpha)/(d\nu)}$. Using the relation $d\nu=2-\alpha$ Huang (2009), we hence obtain $f_{w}\sim\mathcal{O}(1)$. Therefore, for both proposals $f_{w}\sim\mathcal{O}(1)$ in the thermodynamic limit, hence hindering their reliability, which is the same result found in Holubec and Ryabov (2017) for the Otto cycle. Despite $f_{w}\sim\mathcal{O}(1)$, these fluctuations can be suitably avoided when the number $M$ of cycles is large (the argument below applies to the Otto and Carnot cycle). Given $M$ cycles, we have that $f_{w}\propto 1/\sqrt{M}$ as the average work $W\propto M$ whereas the work fluctuations $\sigma_{w}\propto\sqrt{M}$ (think of a biased random walk). Therefore, the ratio between the fluctuations per unit time and the power goes to zero as $M$ grows even when the fluctuations per single cycle are large. Since we have that $M\propto\tau_{\rm tot}/\tau_{\rm eq}$, for a total time $\tau_{\rm tot}$ of observation, fluctuations can be suitably avoided e.g. for critical speed-ups where $\tau_{\rm eq}\propto N^{b}$ with $b<0$.

In actual implementations, the technical requirements to realise such optimised Carnot cycles are: global control of the WS, $H(t)=\lambda(t)H(0)$, and enough precision to engineer small cycles in the region where $\mathcal{C}/\tau_{\rm eq}$ has supralinear scaling with $N$. Since the width $\epsilon$ of this region shrinks with $N$, in a realistic implementation the supralinear scaling will be eventually lost as the control precision is limited 333Other possible implementation-dependent limitations may be the cost of turning on/off the interaction with the baths Perarnau-Llobet et al. (2016); Newman et al. (2017), and a non-zero time for the quenches.. We remark that even when the experimental control may be limited, our considerations provide upper bounds on the maximal power of finite-time Carnot engines.

Conclusions. We have characterized the optimal cycle of a finite-time Carnot engine in the low-dissipation regime. The dissipation has been characterized by the thermodynamic metric (4), which is justified when the thermalization of the working substance (WS) is well described by a single time-scale $\tau_{\rm eq}$. In this case, the optimal cycle turns out to be remarkably simple: it consists of modulations in the form $\lambda(t)H(0)$, where $H(0)$ is the Hamiltonian of the WS. The power output is then proportional to the heat capacity $\mathcal{C}$ of the WS, linking the optimal performance to the nature of the WS: as an application we showed how the critical scaling of $\mathcal{C}$ can enable the design of optimal engines with extensive power reaching Carnot efficiency. These results have been generalised to general metrics in the Supp. Mat. Abiuso and Perarnau-Llobet , which we have used to derive the optimal cycle and corresponding power output of different WS (qubit, 3-level system, or harmonic oscillator) interacting with a bosonic thermal bath. Putting everything together, our reults provide a general framework to efficiently optimise the control of slowly driven Carnot engines.

We hope this work stimulates further investigations in the interplay between many-body physics and heat engines Mascarenhas et al. (2014); Campisi and Fazio (2016); Holubec and Ryabov (2017); Ma et al. (2017); Chand and Biswas (2018); Perarnau-Llobet et al. (2016); Lekscha et al. (2018); Yunger Halpern et al. (2019); Kloc et al. (2019); Jaramillo et al. (2016); Niedenzu and Kurizki (2018), as well as connections between performance, fluctuations, and degree of control Bauer et al. (2016). In particular our results hint at the answer for two open problems: 1) small cycles are optimal for engine performance in all regimes Cavina et al. (2018); Abiuso and Giovannetti (2019); Erdman et al. (2018), and 2) the performance of the proposals Allahverdyan et al. (2013); Campisi and Fazio (2016) cannot be improved.

Acknowledgments. We thank Harry J. D. Miller and Karen Hovhannisyan for insightful discussions. We also thank the anonymous referees for constructive criticisms and useful feedback. P.A. is supported by “la Caixa” Foundation (ID 100010434, fellowship code LCF/BQ/DI19/11730023), the Spanish MINECO (QIBEQI FIS2016-80773-P, and Severo Ochoa SEV-2015-0522), Generalitat de Catalunya (SGR1381 and CERCA Programme), Fundacio Privada Cellex.

Supplemental Material

This Supplemental Material (SM) contains a generalisation of the results presented in the main text to open quantum systems, as well as technical details. We start in Sec. I with a review the quantum Carnot cycle in the quasistatic limit. In Sec. II we introduce finite-time corrections for open quantum systems dynamics and show to characterise them geometrically. In Sec. III we generalise the main result (the optimal cycle and corresponding maximal power) for general open quantum system dynamics by giving the solution as the maximisation of a scalar function. In Sec. IV, we present the maximisation of the heat capacity discussed in the main text. In Sec. IV we discuss the technical details of the asymptotic expansions of power and efficiency presented in the main text. Finally, in Sec. VI we solve analytically an illustrative solvable case. To help the reader find information in this Supplemental Material we introduced Table 1.

I Quasistatic quantum Carnot cycle

For completeness, in this first section we review the quantum Carnot cycle in the quasistatic limit (i.e. without including finite-time corrections).

The internal energy of a system with Hamiltonian $H$ in the state $\rho$ is defined as

[TABLE]

Considering the variation $dU$, it is possible to identify Alicki (1979); Anders and Giovannetti (2013); Kieu (2004); Vinjanampathy and Anders (2016) the work and heat contributions

[TABLE]

To simulate sensible restraints on the system, external control is assumed on dynamical parameters of the local Hamiltonian of the system $H_{t}=H(\vec{l}(t))$. When in contact with a reservoir at temperature $T=1/\beta$ (we use units in which the Boltzmann’s constant is $k_{B}=1$), such a system relaxes to the the Gibbs state

[TABLE]

(here $\mathcal{Z}_{\beta}(H)=\text{Tr}[e^{-\beta H}]$ is the partition function of the system in the canonical ensemble).

A Carnot Cycle Vinjanampathy and Anders (2016); Abiuso and Giovannetti (2019); Cavina et al. (2018, 2017) is identified with a 4 steps process, that is two isothermal strokes alternated with two isoentropic (adiabatic) strokes (cf. Fig. 2). Consider a system with a controlled Hamiltonian $H_{t}$ which can be coupled independently to two reservoirs with temperature $T_{h}>T_{c}$. In the ideal quasistatic limit the operations are performed slowly enough to allow the system to be in thermal equilibrium $\rho(t)\equiv\omega_{\beta}(\vec{l}(t))$ at every instant. The 4 steps are:

while being coupled to the cold reservoir, the Hamiltonian is modified continuously from $H^{(X)}$ to $H^{(Y)}$ such that $\text{Tr}[\dot{\omega}_{\beta}H_{t}]$ is negative, in order for heat to be released to the cold source.

• 2)

with the system isolated from the reservoirs, a quench is performed taking $H^{(Y)}\rightarrow H^{(Y)}\frac{T_{h}}{T_{c}}$.

• 3)

while being coupled to the hot reservoir, the Hamiltonian is modified continuously from $H^{(Y)}\frac{T_{h}}{T_{c}}$ to $H^{(X)}\frac{T_{h}}{T_{c}}$.

• 4)

again isolating the system a quench is performed to restore $H^{(X)}\frac{T_{h}}{T_{c}}\rightarrow H^{(X)}\frac{T_{h}}{T_{c}}\frac{T_{c}}{T_{h}}=H^{(X)}$.

Note the factors $\frac{T_{h}}{T_{c}}$ and $\frac{T_{c}}{T_{h}}$ are chosen in order for the state to be continuous during the quenches. In fact the thermal state uniquely depends on $\beta H_{T}$ (cf. (20)); i.e. for example during the quench 2) the relation $\beta_{c}H^{(Y)}=\beta_{h}H^{(Y)}\frac{T_{h}}{T_{c}}$ guarantees $\omega_{\beta_{c}}(H^{(Y)})=\omega_{\beta_{h}}(H^{(Y)}\frac{T_{h}}{T_{c}})$. We shall thus define then the adimensional Hamiltonian at temperature $1/\beta$ (note that the temperature takes only two values, that depend on the respective isotherms)

[TABLE]

so that the thermal state is $\omega_{t}^{\prime}=e^{-G_{t^{\prime}}}/\text{Tr}[e^{-G_{t^{\prime}}}]$ on both the cold and hot isotherm. The time reparametrization is conceived in order to isolate the shape of the control $G_{t^{\prime}}$, with $0\leq t^{\prime}\leq 1$. Note that $G$ (and hence $\omega$) is continuous also on the quenches; that is, we can consider the cold isotherm consisting in a transformation $(\omega^{(X)},G^{(X)})\rightarrow(\omega^{(Y)},G^{(Y)})$ with duration $\tau_{c}$ and the hot isotherm the opposite $(\omega^{(Y)},G^{(Y)})\rightarrow(\omega^{(X)},G^{(X)})$, with duration $\tau_{h}$. In a Carnot cycle heat is exchanged only during the 1), 3) steps (absorbed from the hot source, released to the cold one), hence we can compute the efficiency using the observation that over a cycle $\Delta Q+\Delta W=0$ and the heat is absorbed from the hot bath $Q^{(0)}_{abs}=Q_{h}^{(0)}$ (we use the superscript*(0)* to indicate quantities in the quasistatic regime)

[TABLE]

where we have used in the last step that in the quasistatic limit the above mentioned considerations together with Eq. (19) give

[TABLE]

II Finite-time Carnot cycle and thermodynamic length

We now consider finite-time corrections on the above quasistatic Carnot cycle through a quantum open system approach, where the dissipation in linear response can be described in a geometric form by following Scandi and Perarnau-Llobet (2019). We consider an isothermal process where the Hamiltonian $H(t)$ of the WS is driven for $t\in[0,\tau]$ in contact with a thermal bath at temperature $T=1/\beta$. In order to characterise the process beyond the quasistatic limit, we need to assume some structure on the thermalization processes of the working substance (WS) induced by the reservoirs. In a rather generic scenario we consider that the relaxation of the WS can be described by a master equation with the thermal state as a unique fixed point:

[TABLE]

where we note that $\mathcal{L}\equiv\mathcal{L}(t)$ is also time-dependent. Following Cavina et al. (2017), the solution of this master equation can be found perturbatively around $1/\tau$. Using the renormalised time $s\in[0,1]$ with the convection $H_{s}\equiv H(s\tau)$ and $\rho_{s}\equiv\rho(s\tau)$, the solution is given by (indicating with a dot $\dot{}$ the time derivative w.r.t. $s$):

[TABLE]

with $\omega_{s}=e^{-\beta H_{s}}/\mathrm{Tr}(e^{-\beta H_{s}})$, and where $\mathcal{L}_{s}^{-1}$ is the inverse of $\mathcal{L}_{s}$ within the traceless subspace of density matrices, the so called Drazin inverse (see e.g. Crooks (2018); Scandi and Perarnau-Llobet (2019); Abiuso and Giovannetti (2019); Cavina et al. (2017)). Plugging this expression into $Q=\int_{0}^{1}ds\hskip 2.84526pt\text{Tr}(\dot{\rho}H)$ and using integration by parts we obtain:

[TABLE]

with the first order correction to the quasistatic limit given by:

[TABLE]

In terms of the adimensional Hamiltonian $G_{s}\equiv\beta H_{s}$, and using the formula for the derivative of an exponential,

[TABLE]

we can write the variation of entropy as

[TABLE]

while $\Sigma$ can be reexpressed in the convenient quadratic form:

[TABLE]

Now, expanding $G_{s}$ as

[TABLE]

we can conveniently write $\Delta S$ and $\Sigma$ as

[TABLE]

with

[TABLE]

The matrix $m_{ij}$ is symmetric, positive-definite due to the second law $d\Sigma\geq 0$, and it depends smoothly on the base point $\omega$; hence it defines a metric.

Finally, we note that the linear expansion (32) can also be obtained in exact treatments of the system-bath Hamiltonian dynamics through linear-response theory Campisi et al. (2012); Bonança and Deffner (2014); Acconcia et al. (2015); Ludovico et al. (2016). As a rule, the expansion (32) can be argued to be general for any system with dissipations that are linear (at the lowest order) in the speed of the driving. Suppose indeed that the dissipation along an infinitesimal segment of the trajectory depends only on the local point and the local driving. As a consequence it must be in the form $d\Sigma=d\lambda_{i}f_{i}(\lambda,\dot{\lambda},\ddot{\lambda},...)$, but the $1/\tau$ scaling implies that the first derivative terms enter linearly in the product while higher orders are suppressed, which implies $d\Sigma=d\lambda_{i}m_{ij}(\lambda)\dot{\lambda}_{j}$, which is equivalent to  (32).

II.1 Thermodynamic metric for single or multiple time scales

We first show how to obtain the Kubo-Mori-Boguliobov metric used in the main text (Eqs. (4) and (8)). This can be easily done by taking an exponential relaxation of $\rho$ to equilibrium with a single time-scale, as described by the Lindbladian:

[TABLE]

which has the Drazin inverse $\mathcal{L}_{s}^{-1}[.]=\tau_{\rm eq}(\omega_{s}\text{Tr}(.)-\mathbb{I})$. In this case, using that $\text{Tr}(\mathfrak{J}_{\omega}[\dot{G}])=0$ one finds that $\Sigma$ in (30) is given by:

[TABLE]

that is Eq. (8) from the main text, as the generalised covariance is given by ${\rm cov}(A,B)=\text{Tr}[A[\mathfrak{J}_{\omega_{s}}[B]]]$. We hence see that the standard thermodynamic length (Eq. (4) in the main text) can be obtained by an heuristic model of thermalisation with a single time-scale.

It is important to keep in mind that these considerations are only relevant close to equilibrium, where the metric (4) in the main text becomes a good thermodynamic description. Furthermore, although we have assumed that the whole state $\rho$ converges to equilibrium as in (34), strictly speaking it is only necessary that the driven observables (the $X_{j}$ in (31)) converge to equilibrium with a single time-scale, i.e., $\langle\dot{X}_{i}\rangle_{\rho}=\tau_{\rm eq}^{-1}(\langle X_{i}\rangle_{\omega}-\langle X_{i}\rangle_{\rho})$, with $\langle X\rangle_{\rho}=\mathrm{Tr}(X\rho)$. This is especially relevant in complex systems (e.g. many-body systems), where the full dissipative dynamics can be extremely complex but the equilibration of some macroscopic observables can be well described by an exponential relaxation with a suitable time-scale. In this sense, it is also worth pointing out that if each generalised observable $X_{i}$ decays with a different time-scale $\tau_{i}$

[TABLE]

then the metric Eq. (4) of the main text can be easily generalised as Scandi and Perarnau-Llobet (2019):

[TABLE]

where we have absorbed the dependence on $\tau_{i}$ in the metric and where $\tau_{i}$ can in principle depend on the point of the trajectory.

As a final remark, we note that while we have derived the metrics (35) and (37) with an heuristic approach based on exponential relaxation near equilibrium, one can derive them from a microscopic derivation based upon linear-response Sivak and Crooks (2012); Bonança and Deffner (2014), in which case $\tau_{i}$ will in general depend on the point of the trajectory, i.e., $G$. This case will be treated in Sec. III.1.2.

II.2 Thermodynamic metric on standard microscopical models

In this section, we apply our general considerations to derive thermodynamic metrics for systems described by quantum master equations. Before, let us discuss here some generic properties. Following Cavina et al. (2017); Breuer et al. (2002), we have in standard scenarios where a quantum system is coupled to a bosonic bath:

[TABLE]

where $\gamma_{0}\nu^{\alpha}$ is the spectral density of the bath ($\alpha$ defines the Ohmicity),

[TABLE]

and

[TABLE]

is the Bose-Einstein distribution (one can consider ferimonic baths by replacing $N(\beta\nu)$ by the Fermi distribution). The general form (38) of $\mathcal{L}$ depends on $H_{s}$, $\beta$ and the spectral density $J(\nu)=\gamma_{0}\nu^{\alpha}$.

Consider now a Carnot cycle with two baths at temperature $\beta_{c}$ and $\beta_{h}$, and the action of each bath being described by a $\mathcal{L}(H_{S},\beta,J(\nu))$. Assuming that the spectral density of both baths is the same, then both Lindbladians are related by the transformation $\beta\rightarrow\lambda^{-1}\beta$ and $H_{S}\rightarrow\lambda H_{S}$, with $\beta\equiv\beta_{c}$ and $\lambda\equiv\beta_{c}/\beta_{h}$. In terms of the Lindbladian, note that Cavina et al. (2017):

[TABLE]

which shows how the generators of the dynamics are related between the cold and hot isotherm, i.e.

[TABLE]

After noticing that the dissipation is related to $\mathcal{L}$ through Eq. (30), we can write a simple proportionality relation for the metric describing the dissipation on the hot and cold isotherm,

[TABLE]

which implies, for time-reversal symmetric protocols described in the main text,

[TABLE]

While from this relation we see that the considerations of the main text (for symmetric protocols $\Sigma_{c}=\Sigma_{h}$) in principle only apply for flat spectral densities $\alpha=0$, we will show in Sec. III.1 that the same figure of merit $\Delta S^{2}/\Sigma_{c}$ (or equivalently $\Delta S^{2}/\Sigma_{h}$) needs to be maximised to obtain maximal power.

An important comment is now in order. Whereas the dissipator (38) is usually derived for time-independent Hamiltonians Breuer et al. (2002), here we are interested in slowly driven Hamiltonians. Nevertheless, the same form for the dissipator is justified as long as the bath dynamics are fast compared to the driving rate of the system Hamiltonian, which leads to the well-known adiabatic master equation Albash et al. (2012); Dann et al. (2018); Yamaguchi et al. (2017).

II.2.1 Two-level system with a bosonic bath

The well-known optical master equation Breuer et al. (2002) can describe a two-level system with Hamiltonian $G\equiv\beta H=\frac{w}{2}\sigma^{z}$ relaxing in a bosonic thermal bath at temperature $1/\beta$, and from (38) and working in the interaction picture it takes the form

[TABLE]

where the rate $\Gamma\propto w^{\alpha}$ depends on the ohmicity of the bath (for a standard atomic-optical field interaction $\alpha=3$ Breuer et al. (2002)) and $\sigma^{\pm}=(\sigma^{x}\pm i\sigma^{y})/2$. It is easy to translate this equation on the single Bloch-vector components

[TABLE]

This equations being in the form (36), imply that the thermodynamic metric is indeed in the form (37). Moreover the covariances of different cartesian components decouple as ${\rm cov}(\sigma^{i},\sigma^{j})\propto\delta_{ij}$ (with $i,j,=x,y,z$) hence implying

[TABLE]

Consider now the external control on the qubit, i.e. ($\hat{n}(t)$ is a unit vector)

[TABLE]

Assuming that the driving is much slower than the internal dynamics of the bath (i.e. the adiabatic master equation Albash et al. (2012); Dann et al. (2018); Yamaguchi et al. (2017)), the same master equation (45) and metric (47) instantaneously holds in the rotated basis where $\sigma^{z}(t)=\hat{n}(t)\cdot\vec{\sigma}$ and with $\omega\rightarrow\omega(t)$ becoming time dependent. The same metric is expressed in polar coordinates in Scandi and Perarnau-Llobet (2019).

II.2.2 Harmonic oscillator with a bosonic bath

As in the previous example, the optical master equation can be written similarly for an harmonic oscillator with gap control $\beta H(t)=w(t)a^{{\dagger}}a$ (again working on the interaction picture and assuming an adiabatic master equation)

[TABLE]

Considering a modulation of the level splitting $w$, we have for the average occupation number

[TABLE]

Thus the metric for the single control parameter $w(t)$ is in this case simply

[TABLE]

Note that a frequency-controlled harmonic oscillator can be realized in a different set up, i.e.

[TABLE]

(to see how this control is not equivalent to (49) it is sufficient to notice $[\dot{H},H]\neq 0$ in general). The metric arising with this control is being worked out in Miller et al. .

II.2.3 Three-level system with a bosonic bath

For a final example we consider here a three-level system satisfying a detailed balance master equation. We assume control on each of the energy levels, while keeping the basis fixed (this is motivated by the observation that creating coherence does not increase power in the linear-response regime Bauer et al. (2016)). Without loss of generality we can assume the ground state energy to be zero ($\beta E_{0}\equiv 0$) and the two excited states with energies $\beta^{-1}{E_{1}(t)}\leq\beta^{-1}E_{2}(t)$, hence the thermal state being characterized by a probability vector

[TABLE]

The evolution of the probability vector is described by the Markovian master equation Breuer et al. (2002)

[TABLE]

The rate matrix $\Gamma_{ij}$ can be obtained from (38), yielding

[TABLE]

where $\Gamma$ is a rate, $N(w)=1/(e^{w}-1)$ for bosonic baths, and $\alpha$ defines the spectral density of the bath. Note that these dynamics cannot be described by the simple exponential decay (36), so that we need to use the more general approach (33) to find the metric. Indeed, we use Eq.(33), with a simplified version of $\mathfrak{J}_{\omega}[\dot{G}]$ due to commuting operators inside it (i.e. in this case Eq.(28) corresponds to $-\int_{0}^{1}dx\hskip 2.84526pt\omega^{1-x}(\dot{G}-\mathrm{Tr}(\omega\dot{G}))\omega^{x}=-\dot{G}\omega+\omega\text{Tr}[\dot{G}\omega]$), obtaining

[TABLE]

where $\tilde{\Gamma}^{-1}_{ij}$ is the Drazin inverse of ${\Gamma}_{ij}$ (see e.g. Crooks (2018); Scandi and Perarnau-Llobet (2019); Abiuso and Giovannetti (2019); Cavina et al. (2017) for details on the Drazin inverse) while $\partial_{j}\omega_{k}$ is the variation of $\omega_{k}$ due to $\dot{E_{j}}$, i.e. $\partial_{j}\omega_{k}=-\omega_{k}\delta_{jk}+\omega_{j}\omega_{k}$. The analytic form of $m_{ij}$ can be computed from this expression even for larger systems with more than 3 levels using symbolic computation software.

III Generalisation of the main result: Optimal cycles for general metrics and asymmetric dissipations

We now have the necessary tools to generalise the optimisation to more generic metrics and thermodynamic processes. In the most general case, we can consider the Carnot cycle where the cold isotherm is characterised by the control parameters $\{\lambda_{j}^{(c)}(s)\}$ and the hot one by $\{\lambda_{j}^{(h)}(s)\}$ in adimensional time units (recall $s\in[0,1]$), as well as the durations $\tau_{c},\tau_{h}$. For time-reversal cycles considered in the main text, $\{\lambda_{j}^{(c)}(s)\}=\{\lambda_{j}^{(h)}(1-s)\}$ (but remind that in general $\tau_{c}\neq\tau_{h}$).

We first note that under optimization of a very general class of processes where $m_{ij}^{(h)}$ and $m_{ij}^{(c)}$ are simply proportional to each other (this class includes e.g. all the standard microscopic scenarios with same spectral density of the baths, as showed in II.2), the time-reversal property is not an assumption but simply follows from the optimization process Cavina et al. (2017). This implies, as we show in the following III.1, that the relevant figure of merit to be maximized is still $\Delta S^{2}/\Sigma$. Therefore, we show how to accomplish this task for different structures of the metric $m_{ij}^{h,c}$, up to the most general case where the dependence on the temperature, spectral density, etc. is encoded in $\mathcal{L}^{-1}$ in the expression (33).

Finally in Section III.2 we will consider generalizations to cases where $m_{ij}^{(h)}$ and $m_{ij}^{(c)}$ differ significantly, so that $\{\lambda_{j}^{(c)}(s)\}\neq\{\lambda_{j}^{(h)}(1-s)\}$.

III.1 Simple asymmetric dissipations

We consider here cases in which the dissipations along the hot isotherm are proportional to the dissipations along the cold isotherm via a constant factor that does not depend on the specific control protocol:

[TABLE]

where $\sigma$ is a number. Importantly this class of cases includes the scenario where the two baths have the same (non-flat) spectral density, as explained in Sec. II.2 and in Cavina et al. (2017), where we obtained that $m^{(c)}_{ij}=({T_{h}}/{T_{c}})^{\alpha}m^{(h)}_{ij}$ (here $m^{(c)}_{ij}$ ($m^{(h)}_{ij}$) is the metric associated to the dissipation of the cold (hot) bath), which implies $\Sigma_{c}=({T_{h}}/{T_{c}})^{\alpha}\Sigma_{h}$.

Starting from (57), the power to maximize is given by

[TABLE]

under the efficiency constraint

[TABLE]

As the only time unit is $\Sigma/\Delta S\equiv\tau^{*}$ and the temperature ratio $r=\frac{T_{c}}{T_{h}}$ the problem can be rephrased as

[TABLE]

where $x_{j}=\tau_{j}/\tau^{*}$. Here the maximization with its constraint is expressed in full adimensional terms, meaning that after solving it the resulting power will be in the form $\frac{T_{h}(\Delta S)^{2}}{\Sigma}f(\gamma,r,\sigma)$ for some scalar function $f$. The general form of the function $f$ can be obtained analytically but it is in general non-trivial and very lengthy (so we do not show it here), but in the limit of high efficiency it simplifies to

[TABLE]

At this point the optimisation of the power boils down to the maximisation of $(\Delta S)^{2}/\Sigma$, as in the main text,. In the following we show how to maximise $\Delta S^{2}/\Sigma$ in general scenarios, and prove that asymptotically infinitesimal cycles are optimal.

Before, let us note that given any reasonable figure of merit between efficiency and power, the relevant figure of merit will always be $(\Delta S)^{2}/\Sigma$. Given as objective any figure of merit $f(\eta,P)$ expressed in terms of the efficiency and the power, after optimization on $\tau_{h},\tau_{c}$ the resulting value $\mathaccentV{bar}016{f}:=f(\mathaccentV{bar}016{\eta},\mathaccentV{bar}016{P})$ (we indicate with $\mathaccentV{bar}016{}$ quantities after time optimization) will be given, by $\mathaccentV{bar}016{\eta}=\mathaccentV{bar}016{\eta}(T_{h},T_{C},\sigma)$ and $\mathaccentV{bar}016{P}=\mathaccentV{hat}05E{P}(T_{h},T_{c},\sigma)(\Delta S)^{2}/\Sigma$ for some adimensional function $\mathaccentV{bar}016{\eta}$ and homogeneous function $\mathaccentV{hat}05E{P}$ of the temperatures, by dimensional analysis. Moreover, any reasonable $f$ will be monotonously increasing in both its parameters singularly (i.e. for a given efficiency we wish to enhance the power and vice-versa), therefore after speed optimization it will still be possible to improve the engine performance by increasing the ratio $(\Delta S)^{2}/\Sigma$. The role of the above mentioned ratio as a characteristic scale defining the performance of both engines and refrigerators was already noticed in Hernández et al. (2015) without further analysis.

III.1.1 Dependence of the time scale in $G$

As a simple first extension of our results, we consider the case of the standard thermodynamic metric with a time-scale that can depend on the point of the parameter space $G=\beta H$. In other words, we consider the heuristic model of exponential thermalisation (34) with a $G$-dependent $\tau_{\rm eq}$ (see also Bonança and Deffner (2014) for a microscopic derivation). Following the reasoning of the main text, we can define an infinite-dimensional scalar product given by

[TABLE]

with $\omega_{s}=e^{-G_{s}}/\mathrm{Tr}(e^{-G_{s}})$, and where we note that the scalar product is defined for a fixed trajectory $G_{s}$ with $s\in[0,1]$. We then have,

[TABLE]

so we conclude that the results of the main text can be naturally extended to time-scales that depend on $\beta H$, meaning that the performance is upper-bounded by small cycles with proportional modulation of the Hamiltonian $\dot{H}\propto H$ performed on the point where the ratio between heat capacity and relaxation time of the system is maximum.

III.1.2 General metrics

Let us consider the maximisation of $\Delta S^{2}/\Sigma$ for more general protocols where the thermodynamic metric is given by the general expression (33) which depends on the particular Lindbladian describing the thermalisation dynamics. Following the considerations of the main text, to maximise $(\Delta S)^{2}/\Sigma$, we can expand them as

[TABLE]

with $s_{i}=\text{Tr}[(\text{Tr}[G\omega]\omega-G\omega)X_{i}]$, and

[TABLE]

In contrast to the main text, where $m_{ij}$ is given by Eq. (4) in the main text, here we a general metric of the form (33). The Cauchy-Schwarz inequality can be applied by considering the following:

Inequality.

Consider two vectors $\vec{a},\vec{b}$ and a quadratic invertible form $g>0$ defined on their vector space. Then the standard C-S inequality applied to $g^{1/2}\vec{a}$ and $g^{-1/2}\vec{b}$ tells

[TABLE]

If we now consider this inequality applied to two vectors

[TABLE]

and the metric

[TABLE]

where $m(t)$ is a positive time-dependent quadratic form in $\mathbb{R}^{k}$ , we have

[TABLE]

and thus it is possible to write the C-S inequality as

[TABLE]

This inequality allows us to bound the figure of merit for the maximum power as

[TABLE]

where we recall that $s_{i}=\text{Tr}[(\text{Tr}[G\omega]\omega-G\omega)X_{i}]$, and the $X_{i}$ are given by the expansion (64). Hence, the initial optimisation over all possible trajectories $\{\lambda_{j}(s)\}$ boils down to maximising a single real function given the control parameters. In other words, our result shows that in order to design optimal finite-time Carnot cycles one does not need to optimise over all trajectories $\{\lambda_{j}(s)\}$ with $s\in[0,1]$, but it is enough to search a suitable point in the thermodynamic space $\{\lambda_{j}\}$. To saturate it in practice, one needs to maximise $\vec{s}m^{-1}\vec{s}$ over the control parameters $\vec{\lambda}$,

[TABLE]

and consider infinitesimal variations around this optimal point:

[TABLE]

where we have used the more compact vector notation: $\vec{s}=\{s_{1},s_{2},...\}$, etc. The direction of the modulations is defined, from the Cauchy-Schwarz saturation condition, by the vector $m^{-1}\vec{s}$, while the normalization of the modulation is in principle to be taken infinitesimal. The abstract expression for $\vec{\mu}$ corresponds, in the case of the analysed in the main text, to $\vec{\lambda}^{*}$ itself, i.e. the modulations of the control parameters are proportional to the control themselves in that case.

III.1.3 Examples

Qubit with optical master equation.

We first consider here a slowly driven qubit with full hamiltonian control (here $G$ is the adimensional Hamiltonian in temperature energy units)

[TABLE]

in contact with a bosonic bath, so that the evolution is described by the adiabatic master equation (in the interaction picture)

[TABLE]

Here, $\sigma^{\pm}\equiv\sigma^{\pm}(t)=\sigma^{x}(t)\pm\sigma^{y}(t)$ in the basis where $\sigma^{z}(t)=\hat{n}(t)\cdot\vec{\sigma}$, and $\omega\equiv\omega(t)$, as in Section II.2.1. Given the metric (47) and considering that the variation of entropy is non-zero only along the instantaneous $z$ direction $\hat{n}(t)$ it is evident how, to maximize $\Delta S/\Sigma^{2}$, changing the direction of $\hat{n}$ is useless (as noted in Scandi and Perarnau-Llobet (2019), generation of coherences is detrimental, and the eigenvectors of $m_{ij}$ (47) are larger along the tangential directions $x$ and $y$), hence the optimal trajectories can be recognized to be in the form of simple gap modulation $w(t)$, reducing the metric to

[TABLE]

Therefore we fall into the category of systems analysed in Section III.1.1, for which the optimal control consists in modulations on the point where the ratio between the heat capacity ${\rm Var}_{\omega}(G)$ and the relaxation timescale is maximum. In this case $\tau_{\rm eq}\equiv\Gamma^{-1}(2N(w)+1)^{-1}$ and $\mathcal{C}_{Q}=w^{2}e^{-w}/(1+e^{-w})-(we^{-w}/(1+e^{-w}))^{2}$. In Fig. 3 we show the maximization of $\mathcal{C}_{Q}\Gamma(2N(w)+1)$, corresponding to the maximum value of $(\Delta S)^{2}/\Sigma$ in this case, choosing $\Gamma$ constant (i.e. flat spectral density, $\alpha=0$). Allowing $\alpha\neq 0$ would not change significantly the difficulty of the maximization problem, here $\alpha=0$ has been chosen to allow a comparison with the case studied in the main text, where the thermalization timescale does not depend on the value of the Hamiltonian.

To sum up, by making use of equation (61), the full maximization of power for a qubit in contact with bosonic thermal baths, leads to

[TABLE]

where the maximization result depends on $\alpha$.

Harmonic oscillator

The above analysis for the qubit can be replicated thoroughly for the controlled harmonic oscillator described in II.2.2, and the final result

[TABLE]

where heat capacity of the harmonic oscillator is equal to $\mathcal{C}_{HO}=(2\sinh(w/2))^{-2}$. In Fig. 4 an example is shown for $\alpha=3$.

3-level system with master equation

Finally we consider a 3-level system as described in II.2.3, controlled by the energy levels of the two excited states $E_{1}$ and $E_{2}$, and whose evolution is given by a standard detailed balanced master equation (54)-(II.2.3). By following the general approach described in III.1.2, we bound $\Delta S^{2}/\Sigma$ by $\sum_{ij}s_{i}m_{ij}^{-1}s_{j}$, see (71). We symbolically compute the metric (56) and the vector $\{s_{i}\}$ (see details in Sec. II.2.3) in order to express $\sum_{ij}s_{i}m_{ij}^{-1}s_{j}$ as a function of $\{E_{1},E_{2}\}$, and then maximise it over $\{E_{1},E_{2}\}$ numerically. An example is shown in Fig. 5, where we plot $\sum_{ij}s_{i}m_{ij}^{-1}s_{j}$ for a linear spectral density $\alpha=1$. Note that the maximum is obtained for $E_{1}=E_{2}$. We also find that in this case the vector $\vec{\mu}$ (73) is proportional to $\begin{pmatrix}1\\ 1\end{pmatrix}$ on the bisector, meaning the modulation on the optimal working point are symmetric on $E_{1}$ and $E_{2}$.

The maximum power results in such a case

[TABLE]

Furthermore, in Fig. 4 (right), we also compare $\sum_{ij}s_{i}m_{ij}^{-1}s_{j}$ with the heat capacity $\mathcal{C}$ divided by the characteristic timescale $\tau_{\rm eq}=\Gamma^{-1}w^{-1}$ (but remember the thermalization process is not simply exponential in this case), and it can be seen how this ratio provides a rather good approximation of the maximal power. That is, while $\mathcal{C}/\tau_{\rm eq}$ is not exactly the correct figure of merit in this case, it can be used as a relevant figure of merit to determine where the optimal power is. This example also illustrates that our approach can be used to deal with more complex dissipative dynamics than the standard exponential decays discussed in the main text.

III.2 General asymmetric protocols

Finally, we consider the most general setting (time-reversal asymmetric protocol, different baths, generic dependence on the temperature of $\mathcal{L}$): $\Sigma_{c}\neq\Sigma_{h}$. Then, the maximum power expression after isotherms-duration tuning (with no efficiency constraints) is given by Ma et al. (2018a)

[TABLE]

which can be thus written as from Eq.s (30)-(33)

[TABLE]

where now two different metrics $m^{(j)}$ appear in the denominator. In this case we lose the simple structure (scalar product$)^{2}/($quadratic form), hence it is not possible to apply directly the C-S inequality; nevertheless it is possible to upper-bound Eq. (81) using simple inequalities e.g.

[TABLE]

These can be now optimized using C-S as we did in the previous section, and depending on the relative size (in the interval $\frac{1}{3}\leq\frac{a}{b}\leq 3$ the former is tighter, otherwise the latter) they will give a bound (not tight in general). Note that while inequalities (82-83) are useful to give upperbounds to the maximum power theoretically obtainable, in practical terms it is also possible to give lower bounds, which are useful to certify that it is possible to reach at least a given value of the power, and maximizing it by the same methods we used in the main work. Namely the bounds can be written

[TABLE]

As a final observation, we note how all of the above can be applied also to the power in the high efficiency regime, which is (cf. Eq. (61))

[TABLE]

IV Maximizing the heat capacity vs. degree of control

The optimal cycles described in the main text (induced by the standard thermodynamic metric (4)) make the maximal power proportional to the heat capacity (13) of the chosen working point of the infinitesimal Carnot cycle. The task of maximizing the power is then translated in finding the point in the control parameter space that maximizes the heat capacity, i.e. the variance of the adimensional Hamiltonian $G$ for a thermal state $\omega=e^{-G}/\text{Tr}[e^{-}G]$. We show in the following 3 paradigmatic cases that are summed up in Fig.6 showing how different degrees of control offer different possible performance on the same system, which is made of $N$ qubits.

(a) Full control over the spectrum. We first assume total control on the Hamiltonian of the $N$ qubits, which is made up of $D=2^{N}$ energy levels. That is, any desired (possibly long-range) interaction can be engineered, so that the $D$-level spectrum can be controlled at will. While this might be extremely challenging to realise in practice, it is useful to consider this situation to obtain a fundamental upper bound on the maximal power. Indeed, the maximization of heat capacity $\mathcal{C}$ of a general $D$-dimensional system at thermal equilibrium has been carried out in (Reeb and Wolf, 2014; Correa et al., 2015). The optimal Hamiltonian consists of a ground level and a $D-1$ degenerate level, with an optimal gap $x$ in adimensional units (i.e. rescaled by the temperature) defined by $e^{x}=(D-1)(x+2)/(x-2)$ and the corresponding $\mathcal{C}$ is $\mathcal{C}_{\rm max}=x^{2}e^{x}(D-1)/(D-1+e^{x})^{2}$. This expression gives in the asymptotic regime ($D\rightarrow\infty$) $x\simeq\ln D$, hence $\mathcal{C}_{\rm max}\simeq(\ln D)^{2}/4$; which in terms of the particle number $N\propto\ln D$ means $\mathcal{C}_{\rm max}\simeq N^{2}/4$ for $N\gg 1$, i.e., a quadratic scaling.

On the robustness of result (a). Suppose the control of the optimal energy gap $x$ has some imprecision (or simply needs to be modulated to perform the cycle). We show here that as far as the modulation/imprecision doesn’t scale with the dimension, the heat capacity behaves smoothly. Suppose indeed that $x=(1+\epsilon)\ln(D-1)$. Then the adimensional variance of the flat Hamiltonian with $d\equiv D-1$ degenerate excited states and gap $x$ is

[TABLE]

which, to keep the scaling as $(\ln d)^{2}\sim N^{2}$ needs the rest to stay finite, i.e. $\epsilon$ to scale as $1/\ln d=1/N$.

(b) N independent qubits. As another extreme case, corresponding to almost no control on the spectrum, we can consider $N$ independent qubits, that is an Hamiltonian $H^{(N)}=\sum_{i=1}^{N}\lambda_{i}(t)\sigma^{z}_{i}$. Note that given the qubits do not interact the thermal heat capacity will result additive, and the optimal gap will be the same for all the qubits $\lambda_{i}=\lambda_{j}$. This means it is enough to solve the previous case for $D=2$ (i.e. $N=1$), to find $C_{\rm max}\approx 0.439N$ (with an optimal gap $\lambda^{*}_{i}\approx 2.40$), in agreement with Abiuso and Giovannetti (2019).

(c) Ising chain. Finally, we consider an Ising chain ($H^{(N)}=-\lambda_{1}(t)\sum_{i=1}^{N}\sigma^{z}_{i}\sigma^{z}_{i+1}-\lambda_{2}(t)\sum_{i=1}^{N}\sigma^{z}_{i}$), and assume control over $\lambda_{1,2}(t)$. Numerical results in Fig. 6 show how the interactions allow for substantially improving on (b), although asymptotically we obtain a linear scaling in $N$, $C_{\rm max}\simeq 0.59N$ for $N\gg 1$. This linear scaling is in fact expected for systems away of a phase transition point (e.g., a rigorous proof that $\mathcal{C}$ is extensive for translationally invariant gapped systems in the low-temperature regime is provided in Appendix A of Hovhannisyan and Correa (2018)). The optimal controls are found to be $\lambda^{*}_{2}\approx 5,\lambda^{*}_{1}\approx 8$.

V Asymptotic expansions and critical scaling

Inspired by Ref. Ma et al. (2018a) we can find the best power for a given fixed efficiency, i.e. \eta=\gamma\eta_{C}=\gamma\big{(}1-\frac{T_{c}}{T_{h}}\big{)}\, which in the symmetric low-dissipation regime means to fix (we identify $\Delta S\equiv\Delta_{h}S$)

[TABLE]

which relates $\tau_{c}$ and $\tau_{h}$, after which it is possible to maximize the power

[TABLE]

by enforcing $\partial_{\tau_{x}}P=0$ (here $\tau_{x}$ can be either $\tau_{c}(\tau_{h})$ or $\tau_{h}(\tau_{c})$ via Eq. (88)), to obtain

[TABLE]

by choosing

[TABLE]

The correspondent work is

[TABLE]

V.1 Critical scaling

Now we consider the optimal engine described in the main text whose cycle consists of $G(s\tau_{c})=(1+\epsilon g(s))G(0)$ with $s\in(0,1)$, $\epsilon\ll 1$ but finite, and $G(\tau_{c}+s\tau_{h})=G((1-s)\tau_{c})$. $G$ is chosen appropriately to maximize the heat capacity given the allowed control. In the limit $\epsilon\ll 1$ the shape $g(s)$ becomes irrelevant as far as it is a smooth function with $g(0)=0$ and $g(1)=1$, which implies $\langle\dot{g}\rangle=1$. This leads to

[TABLE]

We also consider that we scale up the engine with $N$ while approaching $\gamma\rightarrow 1$, by setting (using the notation for critical exponents of a second order phase transition as in Campisi and Fazio (2016); Holubec and Ryabov (2017))

[TABLE]

Then, expanding the relevant quantities for $N\gg 1$, we obtain at leading order in $N$:

[TABLE]

Now we look at the fluctuations of work. The cycle consists of four processes: Two (quasi)-isothermal processes and two quenches. We note that as $N\rightarrow\infty$, we have that $W\rightarrow(T_{h}-T_{c})\Delta S$, so that the isothermal processes become exact at leading order in $N$ which implies that they become fluctuationless Miller et al. (2019); note that this is expected as $\tau_{c}\rightarrow\infty$ with $N\rightarrow\infty$, which makes the low-dissipation assumption also more and more exact as we increase $N$. Secondly, regarding the two quenches: $H\rightarrow HT_{h}/T_{c}$ and the inverse one $H\leftarrow HT_{h}/T_{c}$; it is easy to see that the work fluctuations are given by

[TABLE]

and thus

[TABLE]

Hence we have that

[TABLE]

As explained in the main text, if we want to exploit the critical scaling in the phase transition of the system $\epsilon$ must satisfy $\epsilon\leq N^{-1/(2-\alpha)}$. Then, using $d\nu=2-\alpha$ Huang (2009), we obtain

[TABLE]

which is the same result found in Holubec and Ryabov (2017) for the Otto engine proposal. More details on the necessary critical exponents, and on the fluctuations of the engine, are provided in the main text.

V.2 Comparison with Otto cycle

Let us check in more detail how our results compare to Ref. Campisi and Fazio (2016), where an Otto cycle is considered. Following Campisi and Fazio (2016), let us define the internal energy

[TABLE]

and the heat capacity

[TABLE]

The work output of a Otto cycle working between Hamiltonians $\lambda_{h}H\leftarrow\lambda_{c}H$ for a fixed efficiency $\eta=\gamma\eta_{C}$, $\Delta\eta=\eta_{C}(1-\gamma)$ is given by

[TABLE]

Expanding for low $\Delta\eta$ which corresponds to $\lambda_{c}H_{c}\approx\lambda_{h}H_{h}$, using (102) and keeping only leading terms in $\Delta\eta$ we obtain

[TABLE]

where we defined

[TABLE]

It is important to note that the linear expansion (104) of (103) is only justified close to the phase transition point (105), whose width scales as $\delta\propto N^{-1/(d\nu)}$. This sets an extra requirement on $\Delta\eta$:

[TABLE]

The power is simply $W/\tau$ where $\tau$ is the time of the cycle, which involves two thermalization processes Campisi and Fazio (2016). Let us take $\tau=2\kappa\tau_{\rm eq}$, where $\kappa$ measures how exact the thermalisation process is (the error being exponentially small with $\kappa$ if one assumes a standard exponential relaxation). Then the power reads

[TABLE]

Let us now consider the Carnot cycle. After maximisation over $(\Delta S)^{2}/\Sigma$ we have shown in the main text that

[TABLE]

Expanding around $\gamma\approx 1$ , we have

[TABLE]

which is correct for a fixed $\gamma$ close to 1.

VI Explicit protocol for a $2^{N}-1$ degenerate system

For the sake of completeness, we show here how the results of the paper apply explicitly for a feasible driving protocol close to the optimal one, in the case of a $D$-level system with a ground state and $D-1$ degenerate excited states. As we explained in Sec. IV case (a), this Hamiltonian is motivated as it maximises the heat capacity given a $D$-level system. For comparison purposes, we will sometimes use $N$ satisfying $D=2^{N}$, so that the Hamiltonian can be thought of $N$ suitably interacting particles. Note that the case $N=1$ corresponds to the driving of a qubit Cavina et al. (2018); Abiuso and Giovannetti (2019).

Consider a $D$-level system in a thermal state at temperature $1/\beta$ with an engineered Hamiltonian with a ground state and a $d\equiv D-1$ degenerate excited states with gap $E/\beta$, i.e.

[TABLE]

such that the thermal state is

[TABLE]

where we defined the driving ground state population

[TABLE]

Consider the dynamics given by the equation

[TABLE]

(which induces the metric (4) of the main text in the context of continuous time dynamics, as explained in Section II.1 or Ref. Scandi and Perarnau-Llobet (2019)). We can then consider a solution in the form

[TABLE]

characterized in terms of the ground state probability, which by the slow-driving approximation Cavina et al. (2017) can be solved to be

[TABLE]

that is, expressing quantities in terms of the adimensional time unit $\rho(t)=\rho(s\tau)$ (here $\tau$ is the duration of the driving), and using the notation $\dot{A}=\partial_{s}A=\tau\frac{d}{dt}A$,

[TABLE]

The heat exchange can be computed at all the orders

[TABLE]

The point of maximum heat capacity that optimizes the power output of this system, as shown in Sec. IV case (a), is asymptotically $E\sim\ln d$ which implies $q\sim\frac{1}{2}$. We consider thus a modulation near this point, i.e. a driving protocol in the form q(t)=\frac{1}{2}\big{(}1+\varepsilon\cos(\pi t/\tau)\big{)}, with a small abuse of notation

[TABLE]

or equivalently $E(s)=\ln(\frac{dq(s)}{1-q(s)})$. Notice here that we chose to use the letter $\varepsilon$ instead of $\epsilon$ used in the main text and in Eq. (87), because the latter multiplies the Hamiltonian, which scales linearly in $N$ on the optimal point, while here the modulation of $E(s)$ is of order $\varepsilon$ with no pre-factor. Hence for scaling comparisons it is important to remember $\epsilon\sim\varepsilon/N$.

The total heat can be explicitly computed from (117)

[TABLE]

The low-dissipation we considered in the text consists in neglecting all the terms of order $O(\frac{1}{\tau^{2}\Gamma^{2}})$ on (which in this specific case would consist only in a renormalization of the total result), hence identifying $\beta Q=\Delta S-\Sigma/\tau$:

[TABLE]

In the high efficiency regime we can then write from the asymptotic expansions of Eq. (90) and (91)

[TABLE]

In Figure 7 we show the evolution of the ground state probability $p(s)$ (which characterizes the whole state solution) according to different approximations, as well as the actual performance of the engine compared to the maximal one (which is obtained in the limit $\varepsilon\rightarrow 0$), for feasible choices of the parameters. The slow-driving approximation gives results that almost coincide with the full analytical series for this particular protocol, and the performance is close to the optimal one.