*Article* **Geometric Optimisation of Quantum Thermodynamic Processes**

## **Paolo Abiuso 1, Harry J. D. Miller 2, Martí Perarnau-Llobet 3,\* and Matteo Scandi <sup>1</sup>**


Received: 31 August 2020; Accepted:19 September 2020; Published: 24 September 2020

**Abstract:** Differential geometry offers a powerful framework for optimising and characterising finite-time thermodynamic processes, both classical and quantum. Here, we start by a pedagogical introduction to the notion of thermodynamic length. We review and connect different frameworks where it emerges in the quantum regime: adiabatically driven closed systems, time-dependent Lindblad master equations, and discrete processes. A geometric lower bound on entropy production in finite-time is then presented, which represents a quantum generalisation of the original classical bound. Following this, we review and develop some general principles for the optimisation of thermodynamic processes in the linear-response regime. These include constant speed of control variation according to the thermodynamic metric, absence of quantum coherence, and optimality of small cycles around the point of maximal ratio between heat capacity and relaxation time for Carnot engines.

**Keywords:** quantum thermodynamics; finite-time thermodynamics; thermodynamic length; heat engines; cooling

## **1. Introduction**

Quasistatic processes can be successfully characterised by a few simple and universal results: work is given by the equilibrium free energy difference between the endpoints of a transformation, the efficiency of a Carnot engine depends only on the temperatures of the thermal baths, and in general all quantities of interest become state functions [1]. These results are extremely strong, but their applicability to real life situations is hindered by the necessity of performing all protocols in infinite time in order to ensure that the system remains in thermal equilibrium along the process. On the other hand, finite-time thermodynamic processes can become incredibly complex and strongly depend on the particular protocol and system. For this reason, universal results or simple characterisations are rare. A remarkable exception are fluctuation theorems, which are universal results that apply to arbitrary out-of-equilibrium processes under very mild assumptions [2]; however, they provide a few constraints on the statistics, which are far from sufficient for a full characterisation of the out-of-equilibrium process.

Noticeably, the middle ground between the two situations above, i.e., the case in which the protocol is performed in long but finite time, can be characterised by few geometrical quantities. The main ideas were introduced for classical systems in a series of seminal papers in the 80 s by Weinhold and Andresen, Berry and Salamon, among others [3–14]. More recently, the field saw a revival following a series of papers initiated by Crooks in 2007 [15–17], leading to several applications in, e.g., molecular motors [18], small-scale information processing [19], nonequilibrium

steady states [20,21], and many-body systems [22,23]. The same ideas have been generalised to the quantum regime for unitary dynamics using linear response [24–28], and to open system dynamics for Lindbladian systems [29,30]. Recent applications of thermodynamic geometry in quantum systems can be found in quantum heat engines [31–34], equilibration processes [35,36], phase transitions [37], quantum work and heat fluctuations [38–40], thermodynamic uncertainty relations [41,42], and shortcuts to adiabaticity [43]; see also Ref. [44] for a recent perspective on the subject.

The goal of this paper is two-fold: First, we aim to provide a pedagogic introduction to the notion of (quantum) *thermodynamic length*. This is done in Section 2, where we explicitly connect different frameworks where this concept can be derived: adiabatic linear response theory in closed quantum systems [26–28], adiabatic Lindblad master equations [29,30], and discrete processes [7]. Additionally, in Section 3, we use the concept of thermodynamic length to lower bound the dissipation in a finite-time process, generalising to quantum systems the so-called *Horse–Carrot* theorem [6,7]. Notably, the bound is process-independent, being a function of the endpoints and the (smallest) relaxation timescale. Thus, it can be seen as a geometric refinement of the second law of thermodynamics. Second, in Section 4, we apply these ideas to the optimisation of thermodynamic processes, with emphasis on heat engines in the low-dissipation regime [6,45–53]. Building upon previous works, we show how general conclusions can be drawn with analytical tools for a class of thermal machines, and a few principles of common application can be stated for optimal processes, with some examples. Finally, these results are illustrated in detail for the paradigmatic case of a finite-time Carnot engine with a driven two-level system as a working substance in Section 5.

## **2. Overview of Thermodynamic Length in Quantum Systems**

Let us consider a system whose Hamiltonian *Ht* can be externally driven and which is weakly coupled to a thermal bath. Without loss of generality, we will decompose the system Hamiltonian as *Ht* = ∑*<sup>i</sup> λ<sup>i</sup> tXi*, where - *λi t* is a family of time dependent external parameters, and {*Xi*} are the corresponding observables. Moreover, in the following we will assume summation over repeated indexes. In this context the average work performed on the system is given by:

$$w = \int\_{\gamma} dt \,\mathrm{Tr} \left[ \dot{H}\_{\mathrm{I}} \rho\_{\mathrm{I}} \right] = \int\_{\gamma} dt \,\dot{\lambda}\_{\mathrm{I}}^{\mathrm{i}} \,\mathrm{Tr} \left[ X\_{\mathrm{i}} \rho\_{\mathrm{I}} \right],\tag{1}$$

where *γ* is the path in the parameters space, and *ρ<sup>t</sup>* is the evolved system density matrix at time *t* ∈ (0, *τ*). We know from equilibrium thermodynamics that if the process is infinitely slow the system is always at equilibrium. Consequently, the work is given by the difference of free energy at the endpoints of the transformation. Indeed, in this formalism we regain this result:

$$w\_{\rm eq} = \int\_{\gamma} dt \, \text{Tr} \left[ \dot{H}\_t \pi\_t \right] = \int\_{\gamma} dt \, \frac{d}{dt} \left( -\beta^{-1} \log \mathcal{Z}\_t \right) = \Delta \mathcal{F},\tag{2}$$

where we used the notation <sup>Z</sup>*<sup>t</sup>* <sup>=</sup> Tr *e*−*βHt* for the partition function, we denote the thermal state by *<sup>π</sup><sup>t</sup>* :<sup>=</sup> *<sup>e</sup>*−*βHt*/Z*t*, and we used the definition of the free energy *Ft* :<sup>=</sup> <sup>−</sup>*β*−<sup>1</sup> log <sup>Z</sup>*t*, as well as Δ*F* = *F<sup>τ</sup>* − *F*0. Given this result, it is then natural to define the dissipated work as *w*diss := (*w* − *w*eq) = (*w* − Δ*F*), in order to isolate the role of the dissipation arising from finite time effects.

A consequence of the second law is that *w*diss ≥ 0 with equality only in the infinite time limit. Moreover, if the dynamics is divisible (e.g., Markovian) the rate of dissipation is also positive definite, and zero only in the infinite time limit [54]. This suggests that we can expand *<sup>w</sup>*˙ diss in terms of {*λ*˙ *<sup>i</sup> t*} around the quasistatic limit (*λ*˙ *<sup>i</sup> <sup>t</sup>* ≡ 0), and obtain:

$$
\dot{w}\_{\rm diss} = \lambda\_t^i \underline{\partial}\_{\underline{i}} \underline{\psi}\_{\rm diss}|\_{\overline{\lambda}\_t = 0}^{\cdot - \overline{\omega}^-} + \lambda\_t^i \left( \partial\_i \partial\_j \left. \dot{w}\_{\rm diss} \right|\_{\overline{\lambda}\_t = 0} \right) \lambda\_t^j + \mathcal{O}\left( ||\lambda||^3 \right), \tag{3}
$$

where the first derivative cancels since we are expanding around a minimum. For the same reason, we know that the Hessian *gi*,*<sup>j</sup>* = *β∂i∂<sup>j</sup> w*˙ diss ˙ ˘t≡0 is positive definite. From these considerations we see that the dissipated work can be written as:

$$w\_{\rm diss} = \frac{1}{\beta} \int\_{\gamma} dt \,\, \dot{\lambda}\_t^i \, (\mathcal{g}\_{i,j})\_t \dot{\lambda}\_{t'}^j \tag{4}$$

up to higher order corrections. Linear response theory tells us that the matrix *gt* depends smoothly on the thermal state *πt*. Moreover, we can deduce that it is positive definite and symmetric, being the Hessian of a function around its minimum. These are the defining properties of a metric. In fact, we can interpret Equation (4) as the energy functional or the action of the curve *γ* with respect to the metric *g*. This name comes from the formal analogy between Equation (4) and the action of a system of free particles with mass tensor given by *g*.

This interpretation is particularly useful thanks to the following fact. If one defines the length of *γ* as:

$$d\_{\gamma} = \int\_{\gamma} dt \,\sqrt{\lambda\_t^i \left(\mathcal{g}\_{i,j}\right)\_t \lambda\_{t'}^j} \tag{5}$$

we have the Cauchy–Schwarz like expression

$$
\beta w\_{\rm diss} \ge l\_{\gamma}^2 / \pi,\tag{6}
$$

which takes the name of "thermodynamic length inequality" [6]. Among the curves connecting two endpoints, - *λi* 0 and - *λi τ* , we call *γ* geodesic if it minimises the distance between the two points as measured by Equation (5). A geodesic is also characterised by the property that it keeps the product *λ*˙ *i <sup>t</sup>* (*gi*,*j*)*tλ*˙ *<sup>j</sup> <sup>t</sup>* constant along its path, implying that the Cauchy–Schwarz inequality in Equation (6) is saturated if *γ* is a geodesic. Physically, this means that in order to design minimal dissipating protocols in the slow driving regime, it is sufficient to solve a system of differential equations, i.e., the geodesic equations:

$$
\lambda\_t^i + \Gamma\_{j,k}^i|\_{\lambda\_t} \lambda\_t^j \,\dot{\lambda}\_t^k = 0,\tag{7}
$$

where Γ denotes the Christoffel symbols, which are given by:

$$\Gamma^{i}\_{j,k}|\_{\lambda\_l} = \frac{1}{2} \mathcal{g}^{j,l} \left( \partial\_j \mathcal{g}\_{l,k} + \partial\_k \mathcal{g}\_{j,l} - \partial\_l \mathcal{g}\_{j,k} \right)|\_{\lambda\_l}.\tag{8}$$

Here, *<sup>g</sup>i*,*<sup>l</sup>* is the inverse of the metric, and we use the shorthand notation *<sup>∂</sup>igj*,*k*|*λ<sup>t</sup>* <sup>≡</sup> (*∂gj*,*k*/*∂λi*)|*λ*=*λ<sup>t</sup>* . Moreover, the dissipative properties of a driven system can be directly inferred from the spectral properties of *gt* alone. In particular, starting from very general considerations on the nature of the metric tensor, this will allow us to give lower bounds on the rate of dissipation (Section 3) and to conclude that the creation of coherence is always detrimental to the efficiency (Section 3).

Another strength of the formalism presented is that *g* can be explicitly computed in many frameworks. For example, comparing Equations (1) and (2) it can be seen that the metric tensor can be computed from the slow driving approximation of the expectation value of the observables {*Xi*}s. This was explicitly carried out in the context of linear response of an adiabatically driven unitary dynamics in [28] (see also [26,27]), leading to the expansion:

$$\operatorname{Tr}\left[X\_{l}\rho\_{l}\right] = \operatorname{Tr}\left[X\_{l}\pi\_{l}\right] + \chi\_{t}^{\operatorname{ad}}\left[X\_{i},X\_{j}\right]\dot{\lambda}\_{t}^{\dot{j}} + \mathcal{O}\left(||\dot{\lambda}||^{2}\right),\tag{9}$$

where *χ*ad *<sup>t</sup>* is the adiabatic response function given by:

$$\chi\_t^{\text{ad}}[A,B] = -i \int\_0^\infty d\nu \,\left(\nu \, \text{Tr}\left[ [A(\nu), B] \, \pi\_t \right] \right). \tag{10}$$

Here, we set *h*¯ = 1, and the Heisenberg picture *A*(*s*) is defined with respect to the frozen Hamiltonian at time *t*, i.e., *A*(*s*) = *eiHtsAe*−*iHts* . Notice that the upper bound of the integral can be extended to ∞ thanks to the exponential decay of the correlation function Tr[[*A*(*ν*), *B*]*πt*]. Now, if we plug the expansion just obtained in Equation (1) and we recall that the definition of the dissipated work is *w*diss := (*w* − *w*eq), we have the expression:

$$w\_{\rm diss} = \frac{1}{\beta} \int\_{\gamma} dt \,\lambda\_t^i \left(\beta \,\chi\_t^{\rm ad}[X\_i, X\_j]\right) \lambda\_{t'}^j \tag{11}$$

up to higher order in {*λ*˙ *<sup>i</sup>*}. Comparing this equation with Equation (4), we see that in the context of adiabatic linear response the metric tensor is given by *g*<sup>u</sup> *<sup>i</sup>*,*<sup>j</sup>* <sup>=</sup> *<sup>β</sup>* <sup>2</sup> (*χ*ad *<sup>t</sup>* [*Xi*, *Xj*] + *χ*ad *<sup>t</sup>* [*Xj*, *Xi*]) (notice that even if *χ*ad *<sup>t</sup>* is not in general symmetric in its arguments it can always be symmetrised without affecting the result, since the velocities {*λ*˙ *<sup>i</sup> <sup>t</sup>*} enter the integral in a symmetric way). This formalism was recently used to geometrically characterise thermal machines close to Carnot efficiency [33].

Another relevant framework where a thermodynamic length can be derived is open quantum systems [30] (see also [29]). In particular, consider the Lindbladian dynamics:

$$
\phi\_l = \mathcal{X}\_l^\ell[\rho\_l]\_\prime \tag{12}
$$

with the property that each L*<sup>t</sup>* has the real part of all the eigenvalues negative and that there exist a unique instantaneous steady state *πt*. These two conditions ensure that the dynamics asymptotically equilibrates irrespective of the initial conditions:

$$\lim\_{\nu \to \infty} \varepsilon^{\nu, \mathcal{E}\_l} \rho = \pi\_l. \tag{13}$$

In this case, it is possible to expand the state in the slow driving limit as *ρ<sup>t</sup>* ≈ *π<sup>t</sup>* + *δρ<sup>t</sup>* [55], where *δρ<sup>t</sup>* can be expressed up to higher order corrections as [30]:

$$\rho\_t = \pi\_t + \mathcal{E}\_t^{\prime +} [\dot{\pi}\_t] + \mathcal{O}(||\lambda||^2), \tag{14}$$

where L <sup>+</sup> *<sup>t</sup>* is the Drazin inverse of the Lindbladian given by:

$$\left[\mathcal{QC}\_t^+\left[A\right] = \int\_0^\infty d\nu \, e^{\nu \cdot \mathcal{E}\_t^\circ} \left(\pi\_t \text{Tr}\left[A\right] - A\right). \tag{15}$$

As it will be shown explicitly in the following, the eigenvalues of L <sup>+</sup> *<sup>t</sup>* encode the information about the thermalisation timescales. Moreover, we introduce the shorthand notation to indicate the derivative of the state:

$$\dot{\pi}\_{l} = -\beta \,\dot{\lambda}\_{l}^{i} \int\_{0}^{1} dx \,\,\pi\_{l}^{1-x} \,\, \vec{X}\_{i} \pi\_{l}^{x} = -\beta \,\dot{\lambda}\_{l}^{i} \, \mathbb{J}\_{l}[\vec{X}\_{i}],\tag{16}$$

where we denote by *<sup>X</sup>*¯*<sup>i</sup>* :<sup>=</sup> *Xi* <sup>−</sup> Tr[*Xiπt*]. Hence, if we plug in this expansion into the expression of the work, we obtain that the dissipation takes the form:

$$w\_{\rm diss} = -\frac{1}{\beta} \int\_{\gamma} dt \,\,\dot{\lambda}\_t^i \, (\beta^2 \,\text{Tr} \left[ \mathcal{R}\_i \mathcal{L}\_t^{\rho^+} \mathbb{J}\_t [\mathcal{R}\_j] \right]) \dot{\lambda}\_t^j. \tag{17}$$

Again, it should be noticed that the quadratic form *qi*,*<sup>j</sup>* <sup>=</sup> <sup>−</sup>*β*<sup>2</sup> Tr *X*¯*i*L <sup>+</sup> *<sup>t</sup>* J*t*[*X*¯ *<sup>j</sup>*] is in general not symmetric, so that in the definition of the metric we need to explicitly symmetrise the expression:

*gd <sup>i</sup>*,*<sup>j</sup>* :<sup>=</sup> <sup>1</sup> <sup>2</sup> (*qi*,*<sup>j</sup>* + *qj*,*i*). The matrix *<sup>g</sup><sup>d</sup>* so defined can be then interpreted as the metric tensor for open quantum systems [30].

It is interesting to notice that the metric *g<sup>u</sup>* obtained in the unitary setting can be cast in a form resembling the dissipative one *gd*. In fact, explicitly carrying out the integral in the definition of the adiabatic response function *χ*ad *<sup>t</sup>* , we see that the metric can be recast in the form:

$$\chi\_t^{\text{ind}}[X\_i, X\_l] = -i \int\_0^\infty d\nu \,\left(\nu \, \text{Tr} \left[ [X\_i(\nu), X\_l] \pi\_l \right] \right) = -\frac{i}{\mathcal{Z}\_l} \int\_0^\infty d\nu \,\left(\nu \, \epsilon^{j(x\_H - x\_k)\nu} \right) \left(\epsilon^{-\frac{\rho\_{\text{fS}}}{\rho\_{\text{fS}}}} - \epsilon^{-\frac{\rho\_{\text{fS}}}{\rho\_{\text{fS}}}} \right) (X\_i)\_{m,n} (X\_j)\_{n,m} \tag{18}$$

$$\mathcal{I}\_{t} = -\frac{1}{\mathcal{Z}\_{t}} \frac{(\boldsymbol{\varepsilon}^{-\hat{\beta}\varepsilon\_{H}} - \boldsymbol{\varepsilon}^{-\hat{\beta}\varepsilon\_{I}})}{(\boldsymbol{\varepsilon}\_{m} - \boldsymbol{\varepsilon}\_{n})^{2}} (\boldsymbol{X}\_{i})\_{m,n} (\boldsymbol{X}\_{j})\_{n,m} = -i\boldsymbol{\beta} \int\_{0}^{\infty} d\boldsymbol{\nu} \int\_{0}^{1} d\boldsymbol{x} \, \text{Tr} \left[ \pi\_{t}^{1-x} \boldsymbol{\varepsilon}^{i\hat{H}\_{t}\nu} \boldsymbol{X}\_{i} \boldsymbol{e}^{-i\hat{H}\_{t}\nu} \pi\_{t}^{x} \boldsymbol{X}\_{j} \right] \tag{19}$$

$$\dot{\lambda} = -\beta' \text{Tr} \left[ \mathcal{X}\_t \mathcal{U}\_t^+ \left[ \mathbb{I}\_t \left[ \mathcal{X}\_t \right] \right] \right],\tag{20}$$

where we denoted by {*εi*} the eigenvalues of *Ht*, and we defined the operator:

$$\mathcal{M}\_t^+[A] := -i \int\_0^\infty d\nu \, \text{Tr}\_B[e^{-iH\_t\nu} A e^{iH\_t\nu}].\tag{21}$$

We see that the role of L <sup>+</sup> *<sup>t</sup>* is taken in this case by the map <sup>U</sup> <sup>+</sup> *<sup>t</sup>* , so that the dissipation in the unitary case is given in complete analogy to Equation (17).

One last example that one can consider is the case in which the Hamiltonian is changed in a sequence of quenches, followed by a perfect thermalisation of the system [7]. The total duration of the protocol is given by *τ* = *Nτ*eq, where *N* is the number of quenches in which the protocol is realised and *τ*eq is a fixed equilibration time. When the number of steps is large the state at each time *t* = *mτ*eq (*m* = 0, ... , *N* − 1) is approximately given by: *ρ<sup>m</sup> π<sup>m</sup>* − Δ*mπ*, where Δ*mπ* is the difference between the thermal states at two subsequent steps Δ*mπ* := *πm*+<sup>1</sup> − *πm*. This term in the limit *N* 1 is well approximated by *τ*eq*π*˙ *<sup>t</sup>*. We can interpret this contribution as an indication of how much the system lags behind the thermal state. Proceeding as before, the dissipation can be rewritten up to first order in 1/*N* = *τ*eq/*τ* as:

$$w\_{\rm diss} = \frac{1}{2\beta} \int\_{\gamma} dt \,\,\dot{\lambda}\_t^i \, (\tau\_{\rm eq} \beta^2 \,\mathrm{Tr} \left[ \vec{X}\_i \mathbb{J}\_t [\vec{X}\_j] \right]) \dot{\lambda}\_t^j. \tag{22}$$

The metric tensor *g q <sup>i</sup>*,*<sup>j</sup>* can be directly identified with the trace inside the integral, since <sup>J</sup>*<sup>t</sup>* is self-adjoint, making the whole expression symmetric in (*i*, *j*). The metric so obtained can be rewritten as: *g q <sup>i</sup>*,*<sup>j</sup>* = *<sup>τ</sup>*eq **<sup>g</sup>***BKM <sup>i</sup>*,*<sup>j</sup>* , where we implicitly defined **<sup>g</sup>***BKM <sup>i</sup>*,*<sup>j</sup>* <sup>=</sup> *<sup>∂</sup>*<sup>2</sup> ln <sup>Z</sup>/*∂λi∂λj*. This last quantity is known as the Bogoliubov–Kubo–Mori (BKM) statistical distance, which encodes the geometry of the manifold of Gibbs states and has been thoroughly studied in the literature [56–60]. Due to the formal similarity between (22) and (17), it is insightful to study the relation between both metrics. In [30], it was shown that in the particular case in which the observables of interest {*Yα*} are the left eigenoperators of the Lindbladian, meaning that they evolve according to the equation:

$$\frac{d}{dt}\text{Tr}\left[\boldsymbol{\Upsilon}\_{\boldsymbol{a}}\boldsymbol{\rho}\_{l}\right] = \boldsymbol{\pi}\_{\boldsymbol{a}}^{-1}\left(\text{Tr}\left[\boldsymbol{\Upsilon}\_{\boldsymbol{a}}\boldsymbol{\pi}\_{l}\right] - \text{Tr}\left[\boldsymbol{\Upsilon}\_{\boldsymbol{a}}\boldsymbol{\rho}\_{l}\right]\right),\tag{23}$$

where {*τα*} are the different timescales of the system, the expression of the metric for the Lindbladian dynamics takes the simple form:

$$\mathbf{g}\_{a,\emptyset}^d = \frac{\boldsymbol{\pi}\_a + \boldsymbol{\pi}\_{\emptyset}}{2} \mathbf{g}\_{a,\emptyset}^{BKM},\tag{24}$$

in analogy with the classical result [17]. Since, at least for Lindbladians satisfying detailed balance, {*Yα*} is a complete basis of operators, it is possible to rewrite in this case any observable *Xi* as *Xi* = *ui*,*αYα*. That is, the Lindbladian metric for a general family of observables {*Xi*} is given by:

$$\mathbf{g}\_{i,\mathfrak{j}}^d = u\_{i,a} u\_{j,\mathfrak{k}} \frac{\mathfrak{r}\_a + \mathfrak{r}\_{\mathfrak{f}}}{2} \mathbf{g}\_{a,\mathfrak{f}}^{BKM}. \tag{25}$$

This shows that the role of L <sup>+</sup> *<sup>t</sup>* is to encode the thermalisation timescales of the system, while the main geometrical properties are contained in **g***BKM*. Finally, it should be noticed that in the case of a uniformly thermalising dynamics, i.e., *τα* = *τ*eq ∀*α*, the thermodynamic metric is proportional to the BKM one.

## **3. Bounding Dissipation with Thermodynamic Length**

In a wider context, the BKM metric plays a role within quantum information geometry [61], and can be interpreted as a form of quantum Fisher information [62]. Moreover, it belongs to the family of contractive Riemann metrics over the manifold of normalised density operators *<sup>t</sup>* <sup>=</sup> *t*({*λ<sup>i</sup> <sup>t</sup>*}). A theorem by Petz gives a general characterisation of length between neighbouring quantum states [63]:

$$d\ell^2 = \mathbf{g}\_{i\bar{j}}^f d\lambda^i d\lambda^j \implies \mathbf{g}\_{i\bar{j}}^f = \text{Tr}\left[\frac{\partial \varrho\_t}{\partial \lambda^i} c^f(R\_{\varrho\_{l'}}L\_{\varrho\_l}) \frac{\partial \varrho\_l}{\partial \lambda^j}\right],\tag{26}$$

where *c <sup>f</sup>*(*x*, *y*)=(*y f*(*x*/*y*))−<sup>1</sup> and *f*(*t*) is a so-called Morozova–Cencov function which is operator monotone, normalised such that *f*(1) = 1 and fulfils *f*(*t*) = *t f*(1/*t*). Furthermore *L*, *R* represent the left and right multiplication operators defined according to *L*[*A*] = *A* and *R*[*A*] = *A* respectively [63]. For each different metric we have a different notion of distance between density matrices over a path *γ*:

$$\ell^f(\gamma) := \int\_{\gamma} d\ell = \int\_{\gamma} dt \,\sqrt{\mathbf{g}\_{ij}^f \dot{\lambda}^i \dot{\lambda}^j}. \tag{27}$$

For the particular choice *<sup>f</sup>*(*x*)=(*<sup>x</sup>* <sup>−</sup> <sup>1</sup>)/ log *<sup>x</sup>* one obtains the BKM metric **<sup>g</sup>***<sup>f</sup> ij* = **<sup>g</sup>***BKM ij* , namely

$$\mathbf{g}\_{i\rangle}^{BKM} = \int\_0^1 d\mathbf{x} \, \text{Tr} \left[ \left( \frac{\partial \log \varrho\_t}{\partial \lambda^i} \right) \varrho\_t^x \left( \frac{\partial \log \varrho\_t}{\partial \lambda^j} \right) \varrho\_t^{1-x} \right]. \tag{28}$$

Restricting to the manifold of thermal states *<sup>t</sup>* = *π<sup>t</sup>* we indeed recover the thermodynamic metric in (22). In general, any length of the form (27) is lower bounded by a geodesic path. Notably, analytical expressions for the shortest curves on the density operator manifold for each choice of metric are not known, aside from a couple of examples [64,65] excluding the BKM metric. However, for the BKM statistical length a lower bound is known (Corollary 5.1 of [66]) which depends only on the boundary conditions {*λ<sup>i</sup>* 0}→{*λ<sup>i</sup> τ*}:

$$\ell^{BKM}(\gamma) \ge \mathcal{L}(\varrho\_0, \varrho\_\tau),\tag{29}$$

where

$$\mathcal{L}(\rho, \sigma) = 2 \arccos(\text{Tr}\left[\sqrt{\rho}\sqrt{\sigma}\right]),\tag{30}$$

is the quantum Hellinger angle. We stress that while this bound can always be saturated when the initial and final states commute, transitions between non-commuting states cannot typically saturate (29). Note that in the classical commutative regime, all monotone metrics (26) reduce to the classical Fisher–Rao metric, and a unique geodesic length is singled out by the Hellinger angle between

the initial and final probability distribution [65]. For a pair of discrete classical probability distributions *pn* and *qn*, the Hellinger angle is given by

$$\mathcal{L}(p,q) := 2 \arccos \left( \sum\_{n} \sqrt{p\_n \, q\_n} \right). \tag{31}$$

The geodesic bound (29) has an immediate consequence for thermodynamics. For step-equilibration processes, the work dissipation (22) is subsequently lower bounded via the Cauchy–Schwartz inequality (6) combined with (29):

$$w\_{\rm diss} \ge \frac{k\_B T}{2N} \mathcal{L}^2(\pi\_{0\prime}, \pi\_{\pi}). \tag{32}$$

One may interpret this as a geometric refinement to the second law of thermodynamics. Clearly, the bound depends only on the angle between the initial and final equilibrium state rather than the full path *γ*. For open systems undergoing Markovian dynamics, the corresponding dissipation (17) can be bounded in a similar fashion. Consider first the eigendecomposition of the Lindbladian (23) with associated relaxation timescales {*τα*}, which can be achieved for open systems satisfying detailed balance. Denoting *τ*min as the shortest timescale along the curve *γ* and *τ* the total duration, work dissipation is bounded by

$$k\_{\rm diss} \ge k\_B T \left(\frac{\tau\_{\rm min}}{\tau}\right) \mathcal{L}^2(\pi\_{0\prime}, \pi\_{\tau}).\tag{33}$$

Note that, while (32) can always be saturated by following a geodesic, in general (33) is not tight whenever more than one relaxation timescale is present. The bounds (32) and (33) represent quantum generalisations of the so-called *Horse–Carrot* theorem in finite-time thermodynamics [6,7].

## *Considerations on Coherence Creation*

Now we want to investigate the role of coherence in a a thermodynamic transformation whose dissipation can be described by Equation (17), see also Refs. [39,67]. We start by rewriting the expression for the dissipated work assuming full control on the system Hamiltonian

$$\dot{m}\_{\text{diss}} = -\beta \operatorname{Tr} \left[ \dot{H}\_{\text{t}} \mathcal{L}\_{\text{t}}^{\text{+}} \mathbb{J}\_{\text{\pi}\_{\text{t}}} \dot{H}\_{\text{t}} \right] \equiv \langle \dot{H}\_{\text{t}}, \dot{H}\_{\text{t}} \rangle\_{\text{t}} \,. \tag{34}$$

For notation simplicity we omit the explicit time dependence in this section. We split *H*˙ in its diagonal and coherence parts, with respect the Hamiltonian basis of *<sup>π</sup>* <sup>∝</sup> *<sup>e</sup>*−*βH*, <sup>|</sup>*<sup>i</sup>*

$$
\dot{H} = \dot{H}^{(d)} + \dot{H}^{(c)} \qquad \dot{H}^{(d)} = \sum\_{i} |i\rangle\langle i| \left| \dot{H} \left| i\right\rangle\langle i| \,. \tag{35}
$$

Given that for any operator *<sup>A</sup>* we have Tr *A*(*d*)*A*(*c*) = 0, if we are able to prove that J*<sup>π</sup>* and L <sup>+</sup> do not mix the diagonal and coherent subspaces, then we would have

$$
\langle \dot{H}, \dot{H} \rangle = \langle \dot{H}^{(d)}, \dot{H}^{(d)} \rangle + \langle \dot{H}^{(c)}, \dot{H}^{(c)} \rangle \,. \tag{36}
$$

Now, this is always true for J*<sup>π</sup>* as

$$\mathbb{J}\_{\pi}[|i\rangle\langle j|] = \int\_{0}^{1} d\mathbf{x} \pi^{\mathbf{x}} |i\rangle\langle j| \,\pi^{1-\mathbf{x}} \propto |i\rangle\langle j| \,\tag{37}$$

meaning that if |*i j*| is diagonal (i.e., *i* = *j*), it will stay diagonal, and vice versa (i.e., if *i* = *j*).

Is the same true for L <sup>+</sup>? This question can be answered affirmatively, by noting that L <sup>+</sup> can be written as an exponentiation of L (cf. (17)), and that any L satisfying detailed balance does not mix the diagonal and coherent subspaces [68]. More explicitly, standard Markovian thermal Lindbladians (satisfying detailed balance [68,69]) take the form <sup>L</sup> [*ρ*] = <sup>−</sup>*i*[*HLS*, *<sup>ρ</sup>*] + <sup>∑</sup>*<sup>α</sup> γαAαρA*† *<sup>α</sup>* <sup>−</sup> <sup>1</sup> <sup>2</sup> {*A*† *<sup>α</sup>Aα*, *ρ*}, the *A<sup>α</sup>* being jump operators *A<sup>α</sup>* = |*i<sup>α</sup> jα*|, and *HLS* a general Lamb-Shift Hamiltonian [*HLS*, *H*] = 0. This commutation property guarantees that the Hamiltonian term does not mix populations with coherences, while for the dissipative part we note

$$A\_a \left| i \right> \left< j \right| A\_a^\dagger - \frac{1}{2} \{ A\_a^\dagger A\_{a\prime} \left| i \right> \left< j \right| \} = \left| i\_a \right> \left< i\_a \right| \delta\_{j\_a i} \delta\_{j\_a j} - \frac{1}{2} \left| i \right> \left< j \right| \left( \delta\_{j\_a i} + \delta\_{j\_a j} \right) \,. \tag{38}$$

From the expression above, it is easy to see that if *i* = *j* the result will be diagonal as well, while if *i* = *j* the result will be only made of coherences. Equation (36) is thus valid for standard Markovian master equations and

$$w\_{\rm diss} = w\_{\rm diss}^{(d)} + w\_{\rm diss}^{(c)} \tag{39}$$

where *w*(*d*) diss is the term due to the modification of the spectrum of *<sup>H</sup>*, while *<sup>w</sup>*(*c*) diss is due only to the rotation of the basis. Given that both *w*(*d*) diss and *<sup>w</sup>*(*c*) diss are positive, this property immediately implies that *<sup>w</sup>*diss <sup>≥</sup> *<sup>w</sup>*(*d*) diss, and hence we conclude that the creation of coherence is always detrimental when operating a thermal machine in the low-dissipation regime, as we explain more in detail in Section 4.2, and in agreement with recent results [42,67,70]. A similar separation of losses generated by diagonal and coherent parts of the Hamiltonian variation is presented in [32].

## **4. Optimisation of Thermodynamic Processes in the Slow Driving Regime**

In this section, we derive and review generic considerations on the optimisation of finite-time thermal machines in the low-dissipation regime [6,14,31,46]. That is, when the irreversible entropy production is proportional to the inverse time duration. This assumption can be taken as empiric if no information on the system–bath interaction is given, or it can be justified and derived dynamically using the tools examined in Section 2. Part of the results are in agreement with previous literature and we aim here to collect them in a unified exposition that shows the generality and simplicity hidden in earlier works.

More precisely, we consider a thermal machine made up of a working substance (or machine) and several thermal baths at different temperatures. The level of control consists of *n* experimental parameters of the machine that can be driven (typically Hamiltonian parameters), together with the possibility to put the machine in contact with one of the thermal baths. The *n* control parameters are parametrised as *<sup>λ</sup>*(*s*) <sup>≡</sup> *λs<sup>τ</sup>* with *s* ∈ (0, 1)—note that this notation decouples the duration *τ* of each process from its shape *λ*(*s*). We assume in very general terms that the low-dissipation condition holds and it is described by an underlying thermodynamic metric, as presented in Section 2. That is, for an isothermal transformation at temperature *T* = *β*−1, we rewrite Equation (4) as

$$
\Delta Q = T \left( \Delta S - \frac{\sigma}{\tau} \right) \tag{40}
$$

$$
\sigma = \int\_0^1 ds \,\vec{\lambda}'^T(s) g\_{\vec{\lambda}} \vec{\lambda}'(s) \tag{41}
$$

which follows from identifying *<sup>w</sup>*diss <sup>=</sup> *<sup>w</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>F</sup>* <sup>=</sup> *<sup>T</sup>*Δ*<sup>S</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>Q</sup>* <sup>=</sup> *<sup>T</sup>σ*/*<sup>τ</sup>* and by recalling *<sup>λ</sup>*(*s*) <sup>≡</sup> *λsτ*, which has derivative *<sup>λ</sup>* <sup>≡</sup> *<sup>∂</sup> ∂s <sup>λ</sup>* <sup>=</sup> *<sup>τ</sup>* ˙ *λ*. Notice that in most of what follows, the exact form of *g <sup>λ</sup>* does not significantly change the results. In this sense, most of the derivations are common to any system that has first-order losses described by some quadratic form, as in linear response theory.

We consider a machine performing *M* transformations close to equilibrium (in general with different baths), each described by some heat exchange and some dissipation in the low-dissipation regime, with an output

$$
\Delta W\_{\text{out}} = \sum\_{i}^{M} \Delta Q\_{i} = \sum\_{i=1}^{M} T\_{i} \Delta S\_{i} - \frac{T\_{i} \sigma\_{i}}{\tau\_{i}} \,. \tag{42}
$$

The output being a sum of heat exchanges is guaranteed when considering cycling machines, or when the output of interest is the heat extraction from a subset of the sources. This framework thus includes a variety of tasks: cooling, work extraction, Landauer erasure, Carnot cycles, and generalised Carnot engines with multiple baths or finite size baths (see examples below). In any such a process, three main features can be optimised, corresponding to different levels of control over the machine:


In the following, we elaborate on the above features and show how to optimise them, which can be done independently or sequentially. In particular, following the above order in Section 4.1 we optimize the time duration of each transformation *τ<sup>i</sup>* and show a principle of constant dissipation rate optimality; in Section 4.2 we discuss consequences of the considerations presented in Section 3 when the experimental control is such to allow variations of the curve *γ* defined by *λ*(*s*); and in Section 4.3 we discuss the cases in which a full optimisation can be carried out, so that all the degrees of freedom listed above can be optimised.

## *4.1. Tuning the Speed: Optimality of Constant Dissipation Rate*

Here, we suppose initially that the only control available on the machine (42) is the time tuning of each step *τi*. We wish to maximise the power output *P* = Δ*Wout*/ ∑*<sup>j</sup> τ<sup>j</sup>* for a given loss, or equivalently we fix the (maximum) amount of dissipated work,

$$\sum\_{i} \frac{T\_i \sigma\_i}{\tau\_i} \equiv w\_{\text{diss}} \tag{43}$$

and maximize *P*. The power can be written as

$$P = \frac{\left(\sum\_{i} T\_i \Delta S\_i\right) - w\_{\text{diss}}}{\sum\_{j} \tau\_j},\tag{44}$$

hence, maximising it is equivalent to minimising ∑*<sup>j</sup> τ<sup>j</sup>* with the constraint (43). This can be stated as

**Principle 1.** *Maximising the power at fixed dissipation is equivalent to minimising the dissipation at given duration.*

This remark is important as the main result of this subsection (the optimality of constant thermodynamic speed, or dissipation rate) will thus be valid for all machines performing tasks that are limited by the above trade-off. Examples are: maximising the power, minimising the dissipation (or entropy production) with fixed total time, or hybrid figures of merit combinations, such as maximising the power with a fixed amount of total loss. For a discussion of what machines maximise their outputs when the irreversible entropy production is minimised see [71].

The maximisation of (44) can be done differentiating w.r.t *τ<sup>i</sup>* and using Lagrange multipliers, or directly with a Cauchy–Schwarz inequality

$$w\_{\rm diss} \sum\_{i} \tau\_i = \left(\sum\_{j} \frac{T\_j \sigma\_j}{\tau\_j}\right) \left(\sum\_{i} \tau\_i\right) \ge \left(\sum\_{j} \sqrt{T\_j \sigma\_j}\right)^2 \tag{45}$$

which is saturated when all *Tjσj*/*τ*<sup>2</sup> *<sup>j</sup>* are equal, that is

$$\tau\_{\dot{j}} = \frac{\sqrt{T\_{\dot{j}} \sigma\_{\dot{j}}} (\sum\_{i} \sqrt{T\_{i} \sigma\_{i}})}{w\_{\text{diss}}} \tag{46}$$

$$P\_{w\_{\rm diss}} = \frac{w\_{\rm diss}(\sum\_{i} T\_{i} \Delta S\_{i}) - w\_{\rm diss}^{2}}{(\sum\_{j} \sqrt{T\_{j} \sigma\_{j}^{2}})^{2}} \,. \tag{47}$$

Notice that the fact that *Tjσj*/*τ*<sup>2</sup> *<sup>j</sup>* is the same ∀*j* means that the rate of dissipation is constant for each of the *N* steps of the protocol. In particular, when the dissipation is described by an underlying thermodynamic metric (41), this implies the optimality of constant thermodynamic velocity *T λTg λ λ* = *const*., which can be seen by dividing each transformation into infinitesimal steps, i.e., expressing

$$T\_i \Delta S\_{\vec{i}} - \frac{T\_i \sigma\_{\vec{i}}}{\tau\_{\vec{i}}} = \int\_{\gamma^{(i)}} T d\mathcal{S} - \frac{T d \vec{\lambda}^T \mathcal{g}\_{\vec{\lambda}} d\vec{\lambda}}{d\tau} \tag{48}$$

and applying the above reasoning, which concludes that each of the infinitesimal *Td λ<sup>T</sup> g λd λ <sup>d</sup>τ*<sup>2</sup> must be equal. The "thermodynamic length inequality" inequality (6) ([6,72]) is indeed saturated when its integrand is constant, and coincides with the continuous version of (45). These considerations can be summed up saying that for the class of machines considered here

**Principle 2.** *In optimal protocols, the speed of the control variation is constant (as measured from the underlying thermodynamic metric), leading to a constant entropy production rate.*

The optimality of constant entropy production rate was noted already in the first seminal papers [73] in the context of endoreversible engines, and appeared in many works thereafter (for an historical perspective, see also [74,75]). The above formulation manifests the universality of this principle whenever a trade-off between output rate and losses is present in the regime where losses are linear in the average speed of the process.

The power (46) can be further maximised choosing *w*diss = <sup>1</sup> <sup>2</sup> ∑*<sup>i</sup> Ti*Δ*Si* to obtain the durations leading to the maximum power, in this case

$$P\_{\text{max}} = \frac{(\sum\_{i} T\_{i} \Delta S\_{i})^{2}}{4(\sum\_{j} \sqrt{\sigma\_{j}})^{2}} \,. \tag{49}$$

At maximum power the losses thus correspond to half of the quasistatic output: this corresponds to the "7th principle of control thermodynamics" pointed out by Salamon et al. in [74], whose general validity was unknown: we can state it holds (at least) for all machines described by (42).

We give here an example of application of the time tuning optimisation just described.

## Multi-Bath Carnot Engine

A generalised Carnot engine consists of a sequence of isotherms in contact with different thermal baths, alternated with adiabats as in the standard Carnot cycle. The total work output can be expressed as the sum of the heat exchanges due to cycling conditions, as in Equation (42), with ∑*<sup>i</sup>* Δ*Si* = 0. All the results described above apply and the maximum power obtainable by tuning the time durations of the isotherms is thus as in Equation (49). Moreover, in Appendix A we further analyze this result assuming that all the baths have the same spectral density ∝ *ωα*, described by the ohmicity *α*. Under this hypothesis and the assumption that all the isotherms are small enough (see details in Appendix A), we show how this can be translated in the maximum power being expressed by

$$P\_{\text{max}}^{\text{multi-}-\text{Carnot}} = \frac{(\sum\_{i} T\_{i} dS\_{i})^{2}}{4\kappa\_{0}T\_{0} \left(\sum\_{i} (\frac{T\_{i}}{T\_{0}})^{\frac{1}{2}\overline{\alpha}^{a}} |dS\_{i}|\right)^{2}} \tag{50}$$

where *κ*<sup>0</sup> represents the local ratio between *σ*<sup>0</sup> and (Δ*S*0)<sup>2</sup> at some reference temperature *T*0, and satisfies *κi*/*κ<sup>j</sup>* = (*Ti*/*Tj*)−*α*. In the Appendix A, we show how in this case, the power is upper bounded by the same power when it is obtained by the use of the highest and lowest temperature only, which leads to the maximum power of a standard Carnot Engine (cf. Section 4.3 or [31])

$$P\_{\text{max}}^{\text{multi-Carnot}} \le P\_{\text{max}}^{\text{Carnot}} = \frac{(\Delta S)^2}{\sigma\_h} \frac{(T\_h - T\_c)^2}{4T\_h \left(1 + (\frac{T\_c}{T\_h})^{\frac{1}{2}\overline{a}}\right)^2} \,. \tag{51}$$

#### *4.2. Path Optimisation: Geodesics and Coherences*

When the control over the working fluid allows not only to vary the speed of the transformation, but includes possible modifications of the path *γ* of the trajectory *λ*(*s*), the machine can be substantially improved. The optimisation over *γ* is independent from the time tuning considered in the previous section. It consists of finding the shortest path *σ* = *γ λTg λ λ* between two fixed points for each isotherm (41) considered in the cycle. Indeed, when the extremal points of a trajectory are fixed, the quasistatic output is fixed and minimizing *σ* always improves both power and the efficiency.

More precisely, with the tools described in Section 2, each of the *σ<sup>i</sup>* in Equation (42) will be described as in (5) by some metric *g*(*i*) and some trajectory *<sup>λ</sup>*(*i*), in the form *<sup>σ</sup><sup>i</sup>* <sup>=</sup> *<sup>γ</sup>*(*i*) *λ<sup>T</sup>* (*i*)*g* (*i*) *λ λ* (*i*) . As mentioned earlier (see Section 2 or Section 4.1), by choosing the speed to be constant the above expression can be minimised to the thermodynamic length of the path *γ*(*i*)

$$\sigma\_i = \left( \int\_{\gamma^{(i)}} d\mathbf{s} \, \sqrt{\vec{\lambda}\_{(i)}^T \mathbf{g}\_{\tilde{\lambda}}^{(i)} \vec{\lambda}\_{(i)}'} \right)^2 \equiv l\_{\gamma^{(i)}}^2 \,. \tag{52}$$

This quantity depends only on the path *γ*(*i*) of the trajectory and not on its parametrisation *λ*(*s*), but it can be further minimised by considering its minimum among all the possible paths linking the extremal points, which then defines the geodesics distance between the extremal points

$$d\_{\vec{\lambda}(0), \vec{\lambda}(1)} = \min\_{\substack{\gamma \text{ with externally}}} l\_{\gamma} \tag{53}$$

These considerations can be stated as follows:

**Principle 3.** *In optimal protocols, the driving minimises the entropy production, i.e., it follows a geodesic on the thermodynamic manifold.*

In the quantum case, as showed in Section 3, the irreversible entropy production can be split in two independent parts, one due to the variation of the spectrum *H*˙ (*d*) *<sup>t</sup>* and one due to the rotation of the eigenvectors *H*˙ (*c*) *<sup>t</sup>* of the Hamiltonian, i.e., *<sup>H</sup>*˙ *<sup>t</sup>* <sup>=</sup> *<sup>H</sup>*˙ (*d*) *<sup>t</sup>* <sup>+</sup> *<sup>H</sup>*˙ (*c*) *<sup>t</sup>* and

$$w\_{\rm diss} = w\_{\rm diss}^{(d)} + w\_{\rm diss}^{(c)} \tag{54}$$

where *w*(*X*) diss = −*β dt* Tr *H*˙ (*X*) *<sup>t</sup>* <sup>L</sup> <sup>+</sup> *<sup>t</sup>* <sup>J</sup>*π<sup>t</sup> <sup>H</sup>*˙ (*X*) *t* , with *X* = *d*, *c*. Now, notice that the quasistatic (lossless) output of a thermal machine is given by the integral of the heat exchange, or the work exchange, computed on the equilibrium state *πt*, for example

$$\log w\_{\text{eq}} = \int dt \,\text{Tr} \left[ \pi\_t \dot{H}\_t \right] = \int dt \,\text{Tr} \left[ \pi\_t \dot{H}\_t^{(d)} \right] \,\text{\,\,\,}\tag{55}$$

which shows how the work exchange only depends on the diagonal variation of *H*, that is the spectrum variation. This easily follows from the fact that for thermal states at temperature *T* one has Δ*U* = *w* + Δ*Q* = *w* + *T*Δ*S* , where all the quantities depend uniquely on the spectrum of the final and initial control *H*0, *H<sup>τ</sup>* (which define as well the spectrum of *π*0, *πτ*). This means that given the most general control *Ht* <sup>=</sup> *UtH*(*d*) *<sup>t</sup> <sup>U</sup>*† *<sup>t</sup>* , where *<sup>H</sup>*(*d*) *<sup>t</sup>* is diagonal in a time-independent basis, all the lossless heat and work exchanges are the same for the protocol in which only the spectrum is varied, *H*(*d*) *<sup>t</sup>* . At the same time given *<sup>w</sup>*(*c*) diss <sup>≥</sup> 0, losses are clearly reduced using *<sup>H</sup>*(*d*) *<sup>t</sup>* . From this we learn that, for standard Markovian dissipators,

**Principle 4.** *Quantum coherences are not created in optimal protocols, i.e., non-commutativity* [*Ht*, *Ht*] = 0 *is avoided.*

The effect of coherences inducing losses in the power was noted already in [67] in the context of linear response theory of slowly driven engines with slowly driven temperature, and more recently in [42]. A different approach to quantum dynamics, namely quantum jump trajectories, shows again the detrimental effects of coherence creation [70]. Moreover, notice that if the degree of control on the thermal machine allows to eliminate any coherence creation, using commutative controls all the metrics defined in Equation (26) collapse into the classical one and the geodesics distance between states is given by (31), and the bound (33) can be saturated.

We show here an example of application for a cooling process.

## Cooling/Work Extraction

Suppose we are interested only in a subset of the heat currents that are part protocol, meaning that relevant output is the heat extracted from one (or multiple) thermal sources, as in a generalised refrigerator model. To fix the ideas for a single bath to be cooled the cooling rate is

$$P^{\text{cooling}} = \frac{T\_{\text{c}} \Delta S\_{\text{c}} - \frac{T\_{\text{c}} \sigma\_{\text{c}}}{\tau\_{\text{ex}} + \tau\_{\text{c}}}}{\tau\_{\text{ex}} + \tau\_{\text{c}}} \equiv \frac{T\_{\text{c}} \Delta S\_{\text{c}} - w\_{\text{diss}}}{\tau\_{\text{ex}} + \tau\_{\text{c}}} \tag{56}$$

where now *τex* is additional time spent on parts of the cycle that do not contribute to the cooling output. The optimisation for fixed loss *w*diss applies as from (46) leading to *τ<sup>c</sup>* = *Tcσc*/*w*diss , and a power

$$P\_{\mathcal{w}\_{\text{diss}}}^{\text{cooling}} = \frac{T\_{\text{c}} \Delta S\_{\text{c}} - w\_{\text{diss}}}{\tau\_{\text{c}\text{x}} + T\_{\text{c}} \sigma\_{\text{c}} w\_{\text{diss}}^{-1}} \, \text{} \tag{57}$$

which clearly increases as *σ<sup>c</sup>* is minimised. The overall maximum of the cooling rate becomes for a suitable choice of *w*diss

$$P\_{\text{max}}^{\text{cooling}} = T\_{\text{c}} \sigma\_{\text{c}} \frac{\left(\sqrt{\Delta S\_{\text{c}} \pi\_{\text{c}x}/\sigma\_{\text{c}} + 1} - 1\right)^{2}}{\tau\_{\text{c}x}^{2}} = T\_{\text{c}} \frac{\Delta S\_{\text{c}}^{2}}{4\sigma\_{\text{c}}} - T\_{\text{c}} \frac{\Delta S\_{\text{c}}^{3}}{8\sigma\_{\text{c}}^{2}} \tau\_{\text{c}x} + \mathcal{O}\left(\tau\_{\text{c}x}^{2}\right) \,. \tag{58}$$

The above expressions are all decreasing in the value of *σc*, which is minimal when obtained on the geodesics of the transformation, as from Equations (52) and (53). For example, let us assume that the cooling consists of a single transformation from *π<sup>x</sup>* to *πy*, with no additional time *<sup>τ</sup>ex* <sup>=</sup> 0, and full control on the Hamiltonian defining *<sup>π</sup>x*,*<sup>y</sup>* <sup>=</sup> *<sup>e</sup>*−*Hx*,*y*/*Tc*/Tr *e*−*Hx*,*y*/*Tc* . Then, the maximum cooling power is obtained for a coherence-free protocol [*Hx*, *Hy*] = 0 that leads to *<sup>σ</sup>*min <sup>=</sup> <sup>2</sup>*τ*eq arccos(Tr √*π<sup>x</sup>* √*π<sup>y</sup>* ) from (30), whereas the maximum cooling rate is obtained by substituting it into (58). If the control does not allow for coherence-less transformations, or the Lindbladian has several time-scales, upper bounds on the cooling rate can be obtained by the use of (33).

## *4.3. Choosing the Location: Total Optimisation*

After optimizing the time duration and trajectory of the transformations, the resulting optimal output rates only depend on the end points of the transformations. The final maximisation of such expressions is in general non-trivial. However, we note how the maximum power obtained in (51) is proportional (Δ*S*) <sup>2</sup> /*σ*, which is maximal when *σ* takes the geodesics value described above (53). Thus, this last quantity

$$\frac{(\Delta S)^2}{\sigma} = \frac{\left(\mathcal{S}\_{\vec{\lambda}(0)} - \mathcal{S}\_{\vec{\lambda}(1)}\right)^2}{d\_{\vec{\lambda}(0),\vec{\lambda}(1)}^2} \tag{59}$$

can be maximised by changing the extremal of the transformation. The same quantity appears as the leading term for the cooling rate in (58). We find this to be a strikingly general feature of all thermal machines whose dynamical information ultimately consists of just one simple isothermal transformation close to equilibrium. This is clearly the case for a single heat extraction from a bath as in (58), but it happens also, e.g., for Carnot engines, which, due to the trivial dynamics at the quenches, have all relevant quantities which can be expressed solely in terms of the two isotherms. For example, power and efficiency of a Carnot engine read:

$$P^{\text{Carnot}} = \frac{\Delta S(T\_h - T\_c) - \left(\frac{T\_c \sigma\_c}{\tau\_c} + \frac{T\_h \sigma\_h}{\tau\_h}\right)}{\tau\_c + \tau\_h}, \qquad \eta = \frac{Q\_h + Q\_c}{Q\_h} = 1 - \frac{T\_c(\Delta S + \frac{\sigma\_c}{\tau\_c})}{T\_h(\Delta S - \frac{\sigma\_h}{\tau\_h})}, \tag{60}$$

where Δ*S* is the variation of entropy during the hot isotherm, and the irreversible entropy productions are proportional to each other on optimal protocols *σh*/*σ<sup>c</sup>* = (*Tc*/*Th*)−*α*, according to the spectral density of the baths [31,55] (cf. Appendix A). The two isotherms are thus *symmetric*, in the sense that by construction they have an opposite entropy variation Δ*Sh* = −Δ*Sc*, and the trajectories follow the same geodesics to link the endpoints [31,55]. After time optimisation on *τc*, *τ<sup>h</sup>* in such a case it is clear from dimensional analysis that the resulting power can only be proportional to (Δ*S*)2/*σ<sup>h</sup>* (or equivalently (Δ*S*)2/*σ<sup>c</sup>* due to proportionality) multiplied by a function with the dimension of temperature.

In more detail, it has been shown recently [31] that is possible to express the maximum power at any given efficiency *η* = (1 − *δ*)*η<sup>C</sup>* = (1 − *δ*)(1 − *Tc*/*Th*) for a Carnot engine (see also [51,52]). We report here for simplicity only on the case where *α* = 0, thus *σ<sup>c</sup>* = *σ<sup>h</sup>* = *σ*, as

$$P\_{\delta}^{\text{Carnot}} = \frac{\left(\Delta S\right)^{2}}{4\sigma} \frac{(T\_{\hbar} - T\_{c})^{2}\delta(1 - \delta)}{(1 - \delta)T\_{c} + \delta T\_{\hbar}} \tag{61}$$

The importance of the term (Δ*S*)2/*σ* was noted already in [49] as a natural unit of entropy over time, defining the performance of thermal machines in the low-dissipation regime for any trade-off between power and efficiency. The equivalent optimisation for a refrigerator has been conducted in [76], where one has a cooling power and COP coefficient (this time Δ*S* is defined to be positive on the cold isotherm)

$$P^{\text{Refrigerator}} = \frac{\Delta S T\_{\text{c}} - \frac{T\_{\text{c}} \sigma\_{\text{c}}}{\tau\_{\text{c}}}}{\tau\_{\text{c}} + \tau\_{\text{h}}}, \qquad \varepsilon = \frac{Q\_{\text{c}}}{|Q\_{\text{h}}| - Q\_{\text{c}}} = \frac{T\_{\text{c}} \left(\Delta S - \frac{\sigma\_{\text{c}}}{\tau\_{\text{c}}}\right)}{T\_{\text{h}} \left(\Delta S + \frac{\sigma\_{\text{h}}}{\tau\_{\text{h}}}\right) - T\_{\text{c}} \left(\Delta S - \frac{\sigma\_{\text{c}}}{\tau\_{\text{c}}}\right)}, \tag{62}$$

which leads to a maximum cooling power at given COP (again we report it for flat spectral density *σ<sup>c</sup>* = *σh*, see [76] for generalisations) *ε* = (1 − *δ*)*ε<sup>C</sup>* = (1 − *δ*)*Tc*/(*Th* − *Tc*)

$$P\_{\delta}^{\text{Refrigerator}} = \frac{(\Delta S)^2}{4\sigma} \frac{T\_c (T\_h - T\_c)\delta}{T\_h - \delta T\_c} \,. \tag{63}$$

Crucially, the maximisation of the (Δ*S*)2/*σ* term can always be obtained by the use of a Cauchy–Schwarz inequality [31], that is noticing that

$$\frac{(\int dS)^2}{\int ds \vec{\lambda}^T \mathcal{S}\_{\vec{\lambda}} \vec{\lambda}'} = \frac{\left(\int ds \,\vec{\partial} \mathcal{S}\_{\vec{\lambda}} \cdot \vec{\lambda}'\right)^2}{\int ds \vec{\lambda}'^T \mathcal{S}\_{\vec{\lambda}} \vec{\lambda}'} \le \int ds \,\vec{\partial} \mathcal{S}\_{\vec{\lambda}}^T \mathcal{g}\_{\vec{\lambda}}^{-1} \vec{\partial} \mathcal{S}\_{\vec{\lambda}} \le \max\_{\vec{\lambda}} \,\vec{\partial} \mathcal{S}\_{\vec{\lambda}}^T \mathcal{g}\_{\vec{\lambda}}^{-1} \vec{\partial} \mathcal{S}\_{\vec{\lambda}} \equiv \max\_{\vec{\lambda}} \,\mathcal{C}(\vec{\lambda}) \tag{64}$$

The upper bound in (64) can be saturated by performing an infinitesimal cycles around the point where *C*( *λ*) is maximised. In the meaningful case in which the observables *Xi* decay with a well defined timescale *τ*eq, the dissipation is described by the Kubo-Mori metric (see Section 3), and *C*( *λ*) is exactly the heat capacity of the system divided by the equilibration time, leading to [31]:

$$\frac{(\Delta S)^2}{\sigma} \le \max\_{\mathbf{G}} \frac{\mathcal{L}(\mathbf{G})}{\tau\_{\mathbf{eq}}}.\tag{65}$$

Here, *G* = *βH* is the adimensional Hamiltonian, and the thermal state and the heat capacity can be expressed as *π* = *e*−*G*/Tr *e*−*G* and <sup>C</sup>(*G*) = Tr *G*2*π* − Tr[*Gπ*] 2 . In other words,

**Principle 5.** *In order to optimise the power-efficiency trade-off, perform the finite-time Carnot cycle around the point where the ratio between heat capacity and relaxation time of the working medium is maximised.*

This general principle is illustrated in the next section for a two-level Carnot engine.

## **5. Case Study: Finite-Time Qubit Carnot Engine**

In what follows, we analyse the exactly solvable case of a heat engine where the engine consists of a driven two-level system:

$$H(t) = E(t)\sigma\_z. \tag{66}$$

We consider a finite-time Carnot cycle where the working substance is sequentially connected with two thermal baths at different temperatures (see details of the cycle in [31]), and focus on the low-dissipation regime where the results of Section 4 naturally apply. We model the relaxation with any of the two baths by an exponential decay to equilibrium with timescale *<sup>τ</sup>*eq, Tr[*Hρ*˙] <sup>=</sup> *<sup>τ</sup>*−<sup>1</sup> eq Tr[*H*(*<sup>π</sup>* <sup>−</sup> *<sup>ρ</sup>*)], which corresponds to the so-called reset master equation. In this case, the thermodynamic metric is given by the KMB metric.

Let us define *g* ≡ *βE* (with *β* being the inverse temperature of the bath the working substance is connected to), and let *gx* and *gy* be the two endpoints of the isotherms, with *gx* > *gy*. Let us also introduce the corresponding probabilities of the excited state:

$$p\_x = \frac{e^{-\mathcal{G}\_x}}{1 + e^{-\mathcal{G}\_x}},$$

$$p\_y = \frac{e^{-\mathcal{G}\_y}}{1 + e^{-\mathcal{G}\_y}},\tag{67}$$

with *px* < *py*. Then, we easily obtain:

$$
\Delta S = -p\_y \ln p\_y - (1 - p\_y) \ln(1 - p\_y) + p\_x \ln p\_x + (1 - p\_x) \ln(1 - p\_x). \tag{68}
$$

On the other hand, we can use (33) to lower bound the entropy production in the isothermal processes as:

$$
\sigma \ge \pi\_{\text{eq}} \left( 2 \arccos \left[ \sqrt{p\_{\text{x}} p\_{\text{y}}} + \sqrt{(1 - p\_{\text{x}})(1 - p\_{\text{y}})} \right] \right)^2. \tag{69}
$$

This bound can be saturated by following a geodesic, i.e., a protocol satisfying (7). Putting everything together, we can upper bound the relevant figure of merit (Δ*S*)2/*σ* for the power-efficiency optimisation as:

$$\frac{\left(\Delta S\right)^{2}}{\sigma} \le \frac{\left(-p\_{\mathcal{Y}}\ln p\_{\mathcal{Y}} - (1 - p\_{\mathcal{Y}})\ln(1 - p\_{\mathcal{Y}}) + p\_{\mathcal{X}}\ln p\_{\mathcal{X}} + (1 - p\_{\mathcal{X}})\ln(1 - p\_{\mathcal{X}})\right)^{2}}{\pi\_{\text{eq}}\left(2\arccos\left[\sqrt{p\_{\mathcal{X}}p\_{\mathcal{Y}}} + \sqrt{(1 - p\_{\mathcal{X}})(1 - p\_{\mathcal{Y}})}\right]\right)^{2}}.\tag{70}$$

Importantly, this expression is protocol-independent and can be saturated. Indeed, the maximal power of a finite-time Carnot engine (for a given efficiency *η* = (1 − *δ*)*ηC*) given a two-level system can then be written as (see (61)):

$$\max\_{\gamma} \, P\_{\delta}^{\text{Carnot}} = \frac{1}{4} \frac{(-p\_{\mathcal{Y}} \ln p\_{\mathcal{Y}} - (1 - p\_{\mathcal{Y}}) \ln(1 - p\_{\mathcal{Y}}) + p\_{x} \ln p\_{x} + (1 - p\_{x}) \ln(1 - p\_{x}))^{2}}{\pi\_{\text{eq}} \left( 2 \arccos \left[ \sqrt{p\_{\pi} p\_{\mathcal{Y}}} + \sqrt{(1 - p\_{x})(1 - p\_{y})} \right] \right)^{2}} \, \frac{(T\_{h} - T\_{c})^{2} \delta(1 - \delta)}{(1 - \delta)T\_{\delta} + \delta T\_{h}},\tag{71}$$

where the maximisation is meant over all possible protocols in the slow driving regime. We show the upper bound (70) as a function of *gx* in Figure 1 for various values of *gy*, including the optimal one, *gy* <sup>≈</sup> 2.4. It can be seen that the maximum of (Δ*S*)2/*<sup>σ</sup>* over {*gx*, *gy*} is bounded by the maximum of C/*τ*eq, where C is the heat capacity,

$$\mathcal{C} = \mathrm{g}^2 p (1 - p),$$

$$\mathcal{C} = \mathrm{g}^2 \mathfrak{p} (1 - p) \tag{72}$$

where *p* is the excited state probability *p* = *e*−*g*/(1 + *e*−*g*). This is in full agreement with (65) and [77], and is a particular illustration that the power of finite-time Carnot engines at any efficiency can be bounded by substituting the maximum value of <sup>C</sup>/*τ*eq to (Δ*S*)2/*<sup>σ</sup>* inside expression (61), as discussed in detail in Ref. [31].

**Figure 1.** We plot the upper bound of (Δ*S*)2/*σ*, given in (70), as a function of *gx* for different values of *gy* <sup>=</sup> {0.5, 1.5, 2.4}. The point where *gx* <sup>=</sup> *gy* <sup>≈</sup> 2.4 is the point where (Δ*S*)2/*<sup>σ</sup>* is maximised (this can be easily checked numerically), which is also the point of maximum heat capacity *C*. The heat capacity and its maximum are also plotted in dashed lines. We take *τ*eq = 1.

Summarising, here we have provided a tight upper bound on the relevant figure of merit (Δ*S*)2/*σ* for power (and efficiency) of a finite-time Carnot engine, for the particular case of a two-level driven system. We note that such optimisation for a low-dissipation Carnot cycle or an Otto cycle has been performed in [77], while exact total optimisation for a two-level system performing an arbitrary cycle was solved in Refs. [78,79], with both bosonic and fermionic baths. While our results apply in the high efficiency or low-dissipation regime, their strength lies in its simplicity: indeed, Equation (70) can be easily computed for larger working substances, and extensions to more complex relaxation processes with multiple timescales can also be relatively straightforwardly built (see Equation (64) and Ref. [31]). This contrasts with exact results in finite-time thermodynamics [78,80], which rely on non-trivial optimisation procedures that can become quickly unfeasible as the size of the working substance increases.

## **6. Conclusions and Outlook**

While originally developed for macroscopic systems, the geometric approach to finite-time thermodynamics is now finding renewed applications within the emerging fields of stochastic and quantum thermodynamics. In this paper, we have highlighted its utility for minimising dissipation in small scale systems operating close to equilibrium. We have derived lower bounds on thermodynamic length that provide a geometric refinement to the second law of thermodynamics and allow one to benchmark the attainable efficiency of quantum thermal machines. Alongside this, we summarised a set of key principles needed to optimise finite-time quantum low-dissipation engines in terms of efficiency and power, based on the computation of the thermodynamic metric tensor and length. Taken together, these principles provide a straightforward method for determining optimal thermodynamic processes. Indeed, we have seen that optimality is achieved by ensuring that the cycle follows a geodesic in the parameter space at constant velocity, while minimising the generation of quantum coherence and maximising the heat capacity relative to the relaxation time of the working system.

Interesting future directions for thermodynamic geometry in the quantum regime include the extension beyond the slow driving regime [81], the minimisation and characterisation of work and heat fluctuations [38–40,82], connections with strong coupling and speed-ups to isothermality [83], application to cooling processes and relations with the third law of thermodynamics [84–86], many-body systems and criticality [22,23,37].

**Author Contributions:** Conceptualization, P.A.; H.J.D.M.; M.P.-L. and M.S.; Investigation, P. A.; H.J.D.M.; M.P.-L. and M.S.; Writing—original draft, P.A.; H.J.D.M.; M.P.-L. and M.S.; Writing—review & editing, P.A.; H.J.D.M.; M.P.-L. and M.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** P.A. is supported by "la Caixa" Foundation (ID 100010434, fellowship code LCF/BQ/DI19/11730023). M.P.-L. acknowledges funding from Swiss National Science Foundation (Ambizione PZ00P2-186067). H.J.D.M. acknowledges support from the EPSRC through a Doctoral Prize. M.S. acknowledges funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 713729. Both P.A. and M.S. also acknowledge funding from Spanish MINECO (QIBEQI FIS2016-80773-P, Severo Ochoa SEV-2015-0522), Fundacio Cellex, Generalitat de Catalunya (SGR 1381 and CERCA Programme).

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Appendix A. Optimality of Lowest-Highest Temperature Use in Multi-Bath Carnot Engines**

A generalised, finite-time Carnot engine between multiple thermal sources can be described as in Equation (42) (where the adiabatic steps between the isotherms are assumed to happen on a much shorter timescale and thus neglected when compared to the *τi*s),

$$
\Delta W\_{\text{out}} = \sum\_{i}^{N} \Delta Q\_{i} = \sum\_{i=1}^{N} T\_{i} \Delta S\_{i} - \frac{T\_{i} \sigma\_{i}}{\tau\_{i}} \, , \tag{A1}
$$

with Δ*W* ≥ 0 and where the index *i* runs over multiple thermal baths, possibly with infinitesimal steps, including as a possibility the case in which the reservoirs have finite size and change temperature

during the process (notice that in the case of finite size baths the total dissipation ∑*<sup>i</sup> Tiσ<sup>i</sup> <sup>τ</sup><sup>i</sup>* is the natural measure of efficiency, as the total work extractable from the machine sources is finite and obtainable in the quasistatic regime). All the results of Section 4.1 apply, and the maximum power obtainable after tuning the *τi*s can be written

$$P\_{\text{max}} = \frac{(\sum\_{i} T\_{i} dS\_{i})^{2}}{4(\sum\_{j} \sqrt{T\_{j} \sigma\_{j}})^{2}}.\tag{A2}$$

To analyze further this result, we consider here the following property that holds for simple models where all the baths have the same spectral density

$$
\sigma\_i = \kappa\_0 \left(\frac{T\_i}{T\_0}\right)^{-\kappa} dS\_i^2 \tag{A3}
$$

where *α* represents the spectral density exponent of the baths (their ohmicity), *T*<sup>0</sup> is a reference temperature that can be chosen at will, and *κ*<sup>0</sup> a constant that depends on the local thermal state. This property holds if the steps of the transformation are performed "parallel" to each other and are small enough for the state to be almost always the same. More precisely, baths with the same spectral density satisfy the property

$$\mathbf{g}\_{H\_1}^{(T\_1)} = \left(\frac{T\_1}{T\_2}\right)^{-a} \mathbf{g}\_{H\_2}^{(T\_2)} \quad \text{when} \quad \frac{H\_1}{T\_1} = \frac{H\_2}{T\_2} \,. \tag{A4}$$

Here, *g* is the metric that defines the dissipation in terms of the variation of *dG* ≡ *dH*/*T* (cf. Equation (17)), and the property *H*1/*T*<sup>1</sup> = *H*2/*T*<sup>2</sup> means that the thermal state is the same *π*<sup>1</sup> = *π*2. The absolute value of the variation of entropy is instead the same if *dG*<sup>1</sup> = ±*dG*2, as in such a case

$$|d\mathbf{S}\_1| = |\text{Tr}\left[d\pi\_1\mathbf{G}\_1\right]| = |\text{Tr}\left[d\pi\_2\mathbf{G}\_2\right]| = |d\mathbf{S}\_2|\,. \tag{A5}$$

Combining the above two equations, we obtain (A3). For more details see [55] or the supplementary material of [31]. For such a case we obtain substituting (A3)

$$P = \frac{(\sum\_{i} T\_i dS\_i)^2}{4\kappa\_0 T\_0 \left(\sum\_{i} (\frac{T\_i}{T\_0})^{\frac{1}{2}\frac{d}{2}} |dS\_i|\right)^2} \,. \tag{A6}$$

Moreover, for a cycle we have ∑*<sup>i</sup> dSi* = 0 and we can divide the *N* steps into those having *dSk*<sup>+</sup> > 0 (which we will indicate with the index *k*<sup>+</sup> and those having *dSk*<sup>−</sup> < 0 (with index *k*−). We have thus ∑*k*<sup>+</sup> *dSk*<sup>+</sup> = − ∑*k*<sup>−</sup> *dSk*<sup>−</sup> ≡ S. The power (A12) can then be expressed in terms of the "weights" associated to each step for the positive and negative entropy variations. That is, we define

$$p\_{k^{+}} = \frac{dS\_{k^{+}}}{\mathcal{S}} \qquad p\_{k^{-}} = -\frac{dS\_{k^{-}}}{\mathcal{S}} \tag{A7}$$

The vectors *pk*<sup>+</sup> and *pk*<sup>−</sup> are normalised probability vectors and the power (A12) can be written as

$$4\kappa\_0 T\_0 \vec{P} = \frac{(\sum\_{k^+} T\_{k^+} p\_{k^+} - \sum\_{k^-} T\_{k^-} p\_{k^-})^2}{\left(\sum\_{k^+} (\frac{T\_{k^+}}{T\_0})^{\frac{1-\mathfrak{a}}{2}} p\_{k^+} + \sum\_{k^-} (\frac{T\_{k^-}}{T\_0})^{\frac{1-\mathfrak{a}}{2}} p\_{k^-}\right)^2} = \left(\frac{\vec{T}\_+ \cdot \vec{p}\_+ - \vec{T}\_- \cdot \vec{p}\_-}{\vec{T}\_+' \cdot \vec{p}\_+ + \vec{T}\_-' \cdot \vec{p}\_-}\right)^2 \tag{A8}$$

where we defined 4 positive vectors *T*+, *<sup>T</sup>*−, *T* <sup>+</sup>, *T* <sup>−</sup> <sup>&</sup>gt; 0. Being allowed to modify separately we positive and negative weight (essentially by tuning the size of the entropy variations (A7)) it is possible to maximize the above quantity by noting that for any probability vector *p*, positive vectors *B* > 0, vector *C* , positive constant *b* > 0, and constant *c*, it holds

$$\frac{c + \vec{C} \cdot \vec{p}}{b + \vec{B} \cdot \vec{p}} \le \max\_{i} \frac{c + C\_{i}}{b + B\_{i}} \tag{A9}$$

which is saturated by choosing *pi* = *δi*¯*<sup>i</sup>* , where ¯*i* is the index saturating the maximum of (A9). Applying twice the above inequality to 4*κ*0*T*0*P*¯ of Equation (A8) we obtain

$$\sqrt{4\kappa\_0 T\_0 P} \le \max\_{ij} \frac{T\_{+i} - T\_{-j}}{T'\_{+i} + T'\_{-j}}.\tag{A10}$$

Given that *T* <sup>±</sup>*<sup>i</sup>* <sup>=</sup> *<sup>T</sup>*<sup>±</sup> 1−*α* 2 *<sup>i</sup>* , we study the function

$$f(x,y) = \frac{x-y}{x^6+y^6} \qquad x \ge y \ge 0 \tag{A11}$$

and find that it is always decreasing in *y*. Also, it increases always in *x* provided that *β* ≤ 1. We thus conclude that for *α* ≥ −1 the maximisation on the right-hand side of (A10) is obtained by using the highest and lowest temperature available, that we will call *Th* and *Tc* respectively. We thus find that

$$\bar{P} \le \frac{(T\_h - T\_c)^2}{4\kappa\_0 T\_0 \left( (\frac{T\_h}{T\_0})^{\frac{1-\kappa}{2}} + (\frac{T\_c}{T\_0})^{\frac{1-\kappa}{2}} \right)^2} \tag{A12}$$

which is saturated when *dSc* = −*dSh* and all the rest are null. This shows that under the assumption of equal spectral density the power is bounded by the power obtainable by using only the extremal baths.

## **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
