# Pareto-optimal cycles for power, efficiency and fluctuations of quantum heat
engines using reinforcement learning

Paolo A. Erdman  Freie Universität Berlin, Department of Mathematics and
Computer Science, Arnimallee 6, 14195 Berlin, Germany  Alberto Rolandi
Department of Applied Physics, University of Geneva, 1211 Geneva, Switzerland
Paolo Abiuso  ICFO – Institut de Ciències Fotoǹiques, The Barcelona Institute
of Science and Technology, 08860 Castelldefels (Barcelona), Spain  Martí
Perarnau-Llobet<EMAIL_ADDRESS>Department of Applied Physics,
University of Geneva, 1211 Geneva, Switzerland  Frank Noé  frank.noe@fu-
berlin.de Microsoft Research AI4Science, Karl-Liebknecht Str. 32, 10178
Berlin, Germany  Freie Universität Berlin, Department of Mathematics and
Computer Science, Arnimallee 6, 14195 Berlin, Germany  Freie Universität
Berlin, Department of Physics, Arnimallee 6, 14195 Berlin, Germany  Rice
University, Department of Chemistry, Houston, TX 77005, USA

###### Abstract

The full optimization of a quantum heat engine requires operating at high
power, high efficiency, and high stability (i.e. low power fluctuations).
However, these three objectives cannot be simultaneously optimized - as
indicated by the so-called thermodynamic uncertainty relations - and a
systematic approach to finding optimal balances between them including power
fluctuations has, as yet, been elusive. Here we propose such a general
framework to identify Pareto-optimal cycles for driven quantum heat engines
that trade-off power, efficiency, and fluctuations. We then employ
reinforcement learning to identify the Pareto front of a quantum dot based
engine and find abrupt changes in the form of optimal cycles when switching
between optimizing two and three objectives. We further derive analytical
results in the fast and slow-driving regimes that accurately describe
different regions of the Pareto front.

##### Introduction.

Stochastic heat engines are devices that convert between heat and work at the
nanoscale [1, 2, 3]. Steady-state heat engines (SSHE) perform work against
external thermodynamic forces (e.g. a chemical potential difference) after
reaching a non-equilibrium steady state [4], while periodically driven heat
engines (PDHE) perform work against external driving fields through time-
dependent cycles [5]. The performance of heat engines is usually characterized
by the output power and efficiency, and their optimization has been thoroughly
addressed in literature [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. However, in contrast to their
macroscopic counterpart, the performance of quantum and microscopic engines is
strongly influenced by power fluctuations. While early works have started
optimizing power fluctuation [31, 32, 33, 34, 35], a framework to fully
optimize the performance of microscopic heat engines that accounts for power
fluctuations is currently lacking; this letter fills this void.

An ideal engine operates at high power, high efficiency, and low power
fluctuations; however, such quantities usually cannot be optimized
simultaneously, but one must seek trade-offs. In SSHEs, a rigorous
manifestation of this trade-off is given by the thermodynamic uncertainty
relations [36, 37, 38, 39, 40, 41, 42, 43, 44]. For “classical” stochastic
SSHE (i.e. in the absence of quantum coherence) operating between two
Markovian reservoirs at inverse temperatures $\beta_{\text{C}}$ (cold) and
$\beta_{\text{H}}$ (hot), they read [38]:

$\xi\equiv\frac{2}{\beta_{\text{C}}}\frac{\expectationvalue*{P}}{\expectationvalue*{\Delta
P}}\frac{\eta}{\eta_{\text{c}}-\eta}\leq 1,$ (1)

where $\expectationvalue*{P}$ and $\expectationvalue*{\Delta P}$ are
respectively the average power and power fluctuations, $\eta$ is the
efficiency, and $\eta_{\text{c}}\equiv 1-\beta_{\text{H}}/\beta_{\text{C}}$ is
the Carnot efficiency. Such thermodynamic uncertainty relations imply, for
example, that high efficiency can only be attained at the expense of low power
or high power fluctuations. The thermodynamic uncertainty relation inequality
(1) can be violated with quantum coherence [45, 46, 47, 48, 49, 50, 51, 52,
53] and in PDHE [54, 55, 56, 57, 58, 59]. This has motivated various
generalised thermodynamic uncertainty relations [60, 61, 62, 63, 64, 65], in
particular for time-symmetric driving [55, 40] and slowly driven stochastic
engines [66, 67].

Despite their importance, thermodynamic uncertainty relations provide an
incomplete picture of the trade-off: while high values of $\xi$ may appear
more favorable, this does not give us any information on the individual
objectives. Indeed, [56, 66] have shown that high values of $\xi$ can be
achieved in the limit of vanishing power, while often the goal is to operate
at high power or efficiency.

In this paper, we propose a framework to optimize arbitrary trade-offs between
power, efficiency, and power fluctuations in arbitrary PDHE described by
Lindblad dynamics [68, 69, 70, 71]; this framework enables the use of various
optimization techniques, such as Pontryagin Minimum Principle [72] or
Reinforcement Learning (RL) [73] to find Pareto-optimal cycles, i.e. those
cycles where no objective can be further improved without sacrificing another
one. We then employ RL to fully optimize a quantum dot (QD) based heat engine
[74]. We characterize the Pareto front, i.e. the set of values
$\\{\expectationvalue*{P},\expectationvalue*{\Delta P},\eta\\}$ corresponding
to Pareto-optimal cycles, and evaluate the thermodynamic uncertainty relation
ratio $\xi$ on such optimal cycles. Furthermore, we derive analytical results
for the Pareto front and $\xi$ in the fast [6, 8, 20, 23, 75, 76, 27] and slow
[77, 9, 78, 79, 80, 81, 82, 83, 84, 85, 26, 86, 87, 88] driving regime, i.e.
when the period of the cycle is respectively much shorter or much longer than
the thermalization timescale of the system.

##### Multi-objective optimization of quantum heat engines.

Figure 1: A quantum system (gray circle) is coupled to a thermal bath (left
box) characterized by a controllable inverse temperature $\beta(t)$. The
coupling produces a heat flux $J(t)$. Control parameters $\vec{u}(t)$ allow us
to control the state of the system and the power $P(t)$ extracted from the
system.

We describe a PDHE as a quantum system coupled to a heat bath whose inverse
temperature $\beta(t)$ can be tuned in time between two extremal values
$\beta_{\text{H}}$ and $\beta_{\text{C}}$ with
$\beta_{\text{H}}\leq\beta_{\text{C}}$ (Fig. 1). The coupling produces a heat
flux $J(t)$ from the bath to the quantum system. The PDHE is further
controlled by time-dependent control parameters $\vec{u}(t)$ that allow
exchanging work with the system, producing power $P(t)$. A thermodynamic cycle
is then described by periodic functions $\beta(t)$ and $\vec{u}(t)$ with
period $\tau$. This framework includes standard PDHE, in which the system is
sequentially put in contact with two baths (by abruptly changing the values of
$\beta(t)$) and cases where $\beta(t)$ varies smoothly in time [89, 90, 84,
66, 87, 91, 92]. We assume that the dynamics of the system is described by a
Markovian master equation, i.e that the reduced density matrix $\rho_{t}$ of
the quantum system satisfies

$\dot{\rho}_{t}=\mathcal{L}_{\vec{u}(t),\beta(t)}\rho_{t},$ (2)

where $\mathcal{L}_{\vec{u}(t),\beta(t)}$ is the Lindbladian describing the
evolution of the system [70].

We are interested in characterizing the performance of the PDHE computing the
average power $\expectationvalue*{P}$, average irreversible entropy production
$\expectationvalue*{\Sigma}$ and average power fluctuations
$\expectationvalue*{\Delta P}$ in the asymptotic limit cycle [89, 90, 66],
i.e., in the limit of infinite repetitions of the cycle. In such limit,
$\rho_{t}$ becomes periodic with the same periodicity $\tau$ of the control
(Sup. Mat. [93]).

Given the density matrix $\rho_{t}$, the calculation of
$\expectationvalue*{P}$ and $\expectationvalue*{\Sigma}$ can be performed
using the standard approach first put forward in Ref. [94] (see Sup. Mat. [93]
for details). Defining the time average $\expectationvalue{O}$ of an arbitrary
quantity $O(t)$ as
$\expectationvalue{O}\equiv\tau^{-1}\int_{0}^{\tau}O(t)\,dt,$ we can calculate
$\expectationvalue*{P}$ and $\expectationvalue*{\Sigma}$ by averaging

$\displaystyle P(t)$ $\displaystyle=-\text{Tr}[\rho_{t}\dot{H}_{\vec{u}(t)}],$
$\displaystyle\Sigma(t)$
$\displaystyle=-\text{Tr}[\dot{\rho}_{t}H_{\vec{u}(t)}]\beta(t).$ (3)

Notice that we compute the entropy production $\Sigma(t)\equiv-J(t)\beta(t)$
neglecting the entropy variation $\Delta S$ of the quantum system, since the
periodicity of the state in the limiting cycle implies that $\Delta S=0$ after
each repetition of the cycle.

The average fluctuations $\expectationvalue*{\Delta P}$, however, cannot be
expressed as a time average of a function of the state $\rho_{t}$, since they
involves a two-point correlation function. Indeed, from Ref. [82] we can
express them as

$\expectationvalue*{\Delta
P}=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}dt\Tr[s_{t}\dot{H}_{\vec{u}(t)}],$
(4)

where we define

$s_{t}\equiv\int_{0}^{t}dt^{\prime}P(t,t^{\prime})[\dot{H}_{\vec{u}(t^{\prime})}\rho_{t^{\prime}}]+\expectationvalue{w}_{t}\rho_{t}+h.c.\;.$
(5)

In Eq. (5)
$P(t,t^{\prime})\equiv\overleftarrow{T}\exp[\int_{t^{\prime}}^{t}dt^{\prime\prime}\mathcal{L}_{\vec{u}(t^{\prime\prime}),\beta(t^{\prime\prime})}]$
is the propagator,
$\expectationvalue{w}_{t}\equiv-\int_{0}^{t}dt^{\prime}\Tr[\rho_{t^{\prime}}\dot{H}_{\vec{u}(t^{\prime})}]$
is the total average work extracted between time $0$ and $t$, and $h.c.$
represents the complex conjugate of the right-hand-side.

Here, we overcome the difficulty of computing nested integrals and two-point
correlation function in Eqs. (4) and (5) by noticing that $s_{t}$ is a
traceless Hermitian operator that satisfies

$\dot{s}_{t}=\mathcal{L}_{\vec{u}(t),\beta(t)}s_{t}+\\{\rho_{t},\dot{H}_{\vec{u}(t)}\\}-2\text{Tr}[\rho_{t}\dot{H}_{\vec{u}(t)}]\rho_{t},$
(6)

and becomes periodic with period $\tau$ in the limiting cycle (Sup. Mat.
[93]).

Therefore, by considering $(\rho_{t},s_{t})$ as an “extended state” satisfying
the equations of motion in (2) and (6) in the time-interval $[0,\tau]$ with
periodic boundary conditions, we can compute $\expectationvalue{P}$,
$\expectationvalue{\Sigma}$ and $\expectationvalue*{\Delta P}$ as time-
averages of $P(t)$, $\Sigma(t)$, and $\Delta P(t)$, where $P(t)$ and
$\Sigma(t)$ are defined in Eq. (3), and where

$\Delta P(t)\equiv\Tr[s_{t}\dot{H}_{\vec{u}(t)}].$ (7)

Notice that these are now linear functionals of the extended state.

To identify Pareto-optimal cycles we introduce the dimensionless figure of
merit

$\expectationvalue*{F}=a\frac{\expectationvalue*{P}}{P_{\text{max}}}-b\frac{\expectationvalue*{\Delta
P}}{\Delta
P(P_{\text{max}})}-c\frac{\expectationvalue*{\Sigma}}{\Sigma(P_{\text{max}})},$
(8)

where $a,b,c\geq 0$ are three scalar weights, satisfying $a+b+c=1$, that
determine how much we are interested in each of the three objectives, and
$P_{\text{max}}$, $\Delta P(P_{\text{max}})$ and $\Sigma(P_{\text{max}})$ are
respectively the average power, fluctuations and entropy production of the
cycle that maximizes the power. Notice that, given the relation between
entropy production and efficiency, cycles that are Pareto-optimal for
$\\{\expectationvalue{P},\expectationvalue{\Delta P},\eta\\}$, are also
Pareto-optimal for $\\{\expectationvalue{P},\expectationvalue{\Delta
P},\expectationvalue*{\Sigma}\\}$ (Sup. Mat. [93]). The positive sign in front
of $\expectationvalue*{{P}}$ in Eq. (8) ensures that we are maximizing the
power, while the negative sign in front of $\expectationvalue*{\Delta{P}}$ and
$\expectationvalue*{\Sigma}$ ensures that we are minimizing power fluctuations
and the entropy production. Interestingly, if convex, it has been shown that
the full Pareto front can be identified repeating the optimization of
$\expectationvalue*{F}$ for many values of $a$, $b$, and $c$ [95, 96].

Since $\expectationvalue{F}$ is a linear combination of the average
thermodynamic quantities, using Eqs. (3) and (7) we can express
$\expectationvalue{F}$ as a time integral of a function of the extended state
$(\rho_{t},s_{t})$ and of the controls $\vec{u}(t)$ and $\beta(t)$:

$\expectationvalue{F}=\int_{0}^{\tau}G\left(\rho_{t},s_{t},\vec{u}(t),\beta(t)\right)dt,$
(9)

where $G(\rho_{t},s_{t},\vec{u}(t),\beta(t))$ is a suitable function. The
optimization of $\expectationvalue*{F}$ in this form is precisely the type of
problem that can be readily tackled using optimization techniques such as
Pontryagin Minimum Principle [72] or RL [73]. In this manuscript, we employ
the latter.

##### QD heat engine.

In the following, we compute Pareto-optimal cycles in a minimal heat engine
consisting of a two-level system coupled to a Fermionic bath with flat density
of states. This represents a model of a single-level QD [6, 10, 74]. The
Hamiltonian reads

$H_{u(t)}=u(t)\frac{E_{0}}{2}\sigma_{z},$ (10)

where $u(t)$ is our single control parameter, $E_{0}$ is a fixed energy scale
and $\sigma_{z}$ is a Pauli matrix. Denoting with $\ket{1}$ the excited state
of $H_{u(t)}$, and defining $p_{t}\equiv\matrixelement*{1}{\rho_{t}}{1}$ as
the probability of being in the excited state, the Lindblad equation (2)
becomes $\dot{p}_{t}=-\gamma(p_{t}-\pi_{\vec{u}(t),\beta(t)}),$ where
$\gamma^{-1}$ is the thermalization timescale arising from the coupling
between system and bath, and $\pi_{\vec{u}(t),\beta(t)}=f(\beta(t)E_{0}u(t))$
is the excited level population of the instantaneous Gibbs state, with
$f(x)\equiv(1+e^{x})^{-1}$ [23].

##### Optimal cycles with RL and analytical results.

We optimize $\expectationvalue*{F}$ of the QD heat engine using three
different tools: RL, analytics in the fast-driving regime, and analytics in
the slow-driving regime.

The RL-based method allows us to numerically optimize $\expectationvalue*{F}$
without making any approximations on the dynamics, exploring all possible
(time-discretized) time dependent controls $\beta(t)$ and $u(t)$ subject to
the constraints $\beta(t)\in[\beta_{\text{H}},\beta_{\text{C}}]$ and
$u(t)\in[u_{\text{min}},u_{\text{max}}]$ (thus beyond fixed structures such as
Otto cycles), and identifying automatically also the optimal period. The RL
method, based on the soft actor-critic algorithm [97, 98, 99] and generalized
from [100, 101], additionally includes the crucial impact of power
fluctuations, and identifies Pareto-optimal cycles (Sup. Mat. [93] for
technical details and for benefits of using RL). Machine learning methods have
been employed for other quantum thermodynamic [102, 103, 104] and quantum
control [105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117]
tasks.

The fast-driving regime assumes that $\tau\ll\gamma^{-1}$. Interestingly,
without any assumption on the driving speed, we show [93] that any trade-off
between power and entropy production ($b=0$ in Eq. (8)) in the QD engine is
maximized by Otto cycles in the fast-driving regime, i.e. switching between
two values of $\beta(t)$ and $u(t)$ “as fast as possible” [23, 27]. We thus
expect such “fast-Otto cycles” to be nearly optimal in the high power or
efficiency regime. Moreover, we derive analytical expressions to compute and
optimize $\\{\expectationvalue{P},\expectationvalue{\Delta
P},\expectationvalue*{\Sigma}\\}$ efficiently in arbitrary systems in the fast
driving regime [93].

The slow-driving regime corresponds to the opposite limit, i.e.
$\tau\gg\gamma^{-1}$. Since entropy production and power fluctuations can be
minimized by considering quasi-static cycles (see e.g. [56, 66]), we expect
this regime to be nearly optimal in the low power regime, i.e. for low values
of $a$ in Eq. (8). To make analytical progress in this regime, we maximize Eq.
(8) assuming a finite-time Carnot cycle (c.f. Sup. Mat. [93]). The obtained
results naturally generalize previous considerations for low-dissipation
engines [9, 10, 13, 14, 21, 22, 26] to account for the role of fluctuations
(see also [32]). The main technical tool is the geometric concept of
“thermodynamic length”[77, 86] which yields the-first order correction in
$(\gamma\tau)^{-1}$ from the quasi-static limit.

Figure 2: Optimization of $\expectationvalue*{F}$ at different values of $a$
and $c$, with $b=1-a-c$, for a QD-based PDHE. Each dot in panel (a) displays,
as a function of $c$ and $a$, whether
$\expectationvalue*{F}_{\text{RL}}>\expectationvalue{F}_{\text{FAST}}$ (blue
dots) or not (red dots). Points with $a\sim 0$ are not displayed since, in
such a regime, optimal cycles become infinitely long (to minimize entropy
production and fluctuations) and the RL method does not converge reliably
[93]. (b): contour plot of $\expectationvalue*{F}_{\text{RL}}$, as a function
of $c$ and $a$, using the data-points of (a). The black line represents the
curve below which $\expectationvalue*{F}_{\text{FAST}}=0$. (c,d,e): cycles,
described by piece-wise constant values of $\beta$ (black dots) and $u$ (blue
dots) as a function of $t$, identified at the three values of $a$ and $c$
highlighted in black in panel (a) (respectively from top to bottom). The inset
in panel (c) represents a zoom into the corresponding cycle, which is a fast-
Otto cycle. Parameters: $\beta_{\text{C}}=2$, $\beta_{\text{H}}=1$,
$u_{\text{min}}=0.2$, $u_{\text{max}}=1.1$ and $E_{0}=2.5$.

We now present the results. Each point in Fig. 2a corresponds to a separate
optimization of $\expectationvalue*{F}$ with weights $c$ and $a$ displayed on
the x-y axis. Since $b=1-a-c$, points lying on the sides of the triangle
(highlighted in yellow) correspond to optimizing the trade-off between 2
objectives, whereas points inside the triangle take all 3 objectives into
account. Denoting the figure of merit optimized with RL and with fast-Otto
cycles with $\expectationvalue*{F}_{\text{RL}}$ and
$\expectationvalue*{F}_{\text{FAST}}$, in Fig. 2a we show blue (red) dots when
$\expectationvalue{F}_{\text{RL}}>\expectationvalue{F}_{\text{FAST}}$
($\expectationvalue{F}_{\text{RL}}\leq\expectationvalue{F}_{\text{FAST}}$),
while Fig. 2b is a contour plot of $\expectationvalue{F}_{\text{RL}}$. As
expected, there are red dots when $b=0$ (along the hypotenuse), but it turns
out that fast-Otto cycles are optimal also when $c=0$. However, as soon as all
3 weights are finite, the optimal cycles identified with RL change abruptly
and outperform fast-Otto cycles. Furthermore, we notice that while
$\expectationvalue{F}_{\text{RL}}$ is positive for all values of the weights,
$\expectationvalue{F}_{\text{FAST}}=0$ below the black curve shown in Fig. 2b
(see Sup. Mat. [93] for its analytic expression).

To visualize the changes in protocol space, in Figs. 2(c,d,e) we show the
cycles identified with RL at the three different values of the weights
highlighted by a black circle in Fig. 2a (respectively from top to bottom).
Since RL identifies piece-wise constant controls, the cycle is displayed as
dots corresponding to the value of $\beta(t)$ (black dots) and $u(t)$ (blue
dots) at each small time-step. First, we notice that the inverse temperature
abruptly switches between $\beta_{\text{H}}$ and $\beta_{\text{C}}$ for all
values of the weights, so that in this engine no gain arises when smoothly
varying the temperature. As expected, the cycle identified by RL in Fig. 2c,
corresponding to the black point on the hypotenuse in Fig. 2a, is a fast-Otto
cycle (a “zoom” in a short time interval is shown in the inset). However,
moving down in weight space to the black dot at $a=0.6$ and $c=0.2$, we see
that the corresponding cycle (Fig. 2d) now displays a finite period, with
linear modulations of $u(t)$ at fixed temperatures, and a discontinuity of
$u(t)$ when switching between $\beta_{\text{H}}$ and $\beta_{\text{C}}$. The
cycle in Fig. 2e, corresponding to the lowest black dot at $a=0.2$ and
$c=0.4$, displays an extremely long period $\tau\approx 125\gamma^{-1}$, which
is far in the slow-driving regime. Optimal cycles, therefore, interpolate
between the fast and the slow-driving regimes as we move in weight space (Fig.
2a) from the sides to the lower and central region (i.e. switching from 2 to 3
objectives).

Figure 3: Pareto-front found optimizing $\expectationvalue*{F}$ with fast-Otto
cycles in the limit of small temperature differences [panel (a)], optimizing
$\expectationvalue*{F}$ in the slow-driving regime [panel (b)], and
numerically using RL [panel (c)]. The system parameters are as in Fig. 2. All
panel display $\expectationvalue*{\Delta{P}}/\Delta P(P_{\text{max}})$ as a
function of $\expectationvalue*{{P}}/P_{\text{max}}$ (x-axis) and of
$\eta/\eta_{\text{c}}$ (color). The black curve represents the outer border of
the Pareto-front derived analytically in Sup. Mat. [93].

In Fig. 3 we display the Pareto-front, i.e. we plot the value of
${P}/P_{\text{max}}$, $\Delta{P}/\Delta P(P_{\text{max}})$, and of the
efficiency $\eta/\eta_{\text{c}}$, found maximizing $\expectationvalue*{F}$
for various values of the weights. Figure 3a is derived in the fast-driving
regime assuming a small temperature difference, while Fig. 3b is derived in
the slow-driving regime. The RL results, shown in Fig. 3c, correspond to the
points in Fig. 2a. First, we notice that, by definition of the Pareto front,
the “outer border” corresponds to points where we only maximize the trade-off
between the two objectives $\expectationvalue*{P}$ and
$\expectationvalue*{\Delta{P}}$. Since these points are optimized by fast-Otto
cycles, the black border of Fig. 3a, also shown in Figs. 3b, 3c, is exact and
given by (see Sup. Mat. [93] for details)

$\frac{\expectationvalue*{P}}{P_{\text{max}}}=2\sqrt{\frac{\expectationvalue*{\Delta
P}}{\Delta P(P_{\text{max}})}}-\frac{\expectationvalue*{\Delta P}}{\Delta
P(P_{\text{max}})}.$ (11)

Moreover, in this setup, we can establish an exact mapping between the
performance of a SSHE and of our PDHE operated with fast-Otto cycles (cf. Sup.
Mat. [93]). Since SSHE satisfy Eq. (1), also fast-Otto cycles have $\xi\leq
1$. Furthermore, for small temperature differences, $\xi=1$. This allows us to
fully determine the internal part of the Pareto front in the fast-driving
regime using the thermodynamic uncertainty relations, i.e.
${P}/P_{\text{max}}=(\Delta{P}/\Delta
P(P_{\text{max}}))(\eta_{\text{c}}-\eta)/\eta$. Indeed, the linear contour
lines in Fig. 3a stem from the linearity between ${P}$ and $\Delta{P}$, the
angular coefficient being determined by the efficiency.

Comparing Figs. 3(a,b), we see where the fast and slow-driving regimes are
optimal. As expected, the slow-driving Pareto front cannot reach the black
border, especially in the high-power area, where fast-Otto cycles are optimal.
However, in the low power and low fluctuation regime, cycles in the slow-
driving substantially outperform fast-Otto cycles by delivering a higher
efficiency (purple and blue regions in Fig. 3b).

Figure 4: (a): contour plot of the SSHE thermodynamic uncertainty relationship
ratio $\xi$ as a function of $c$ and $a$. (b,c,d): log-log plot of $\xi$,
color mapped as in panel (a), as a function of
$\expectationvalue*{{P}}/P_{\text{max}}$,
$\expectationvalue*{\Delta{P}}/\Delta P(P_{\text{max}})$ and
$\expectationvalue*{\Sigma}/\Sigma(P_{\text{max}})$, respectively. Every point
corresponds to the same RL optimization performed in Fig. 2. The black line is
the behavior of $\xi$ derived analytically in Sup. Mat. [93] in the slow-
driving regime for small values of $\expectationvalue*{{P}}$,
$\expectationvalue*{\Delta{P}}$ and $\expectationvalue*{\Sigma}$.

Interestingly, the RL points in Fig. 3c capture the best features of both
regimes. RL can describe the high-power and low fluctuation regime displaying
both red and blue/green dots near the lower border. The red dots are fast-Otto
cycles that are optimal exactly along the border but deliver a low efficiency.
The blue/green dots instead are finite-time cycles that deliver a much higher
efficiency by sacrificing a very small amount of power and fluctuations. This
dramatic enhancement of the efficiency as we depart from the lower border is
another signature of the abrupt change in optimal cycles.

_Violation of thermodynamic uncertainty relation_. At last, we analyze the
behavior of the thermodynamic uncertainty relation ratio $\xi$, which
represents a relevant quantity combining the three objectives, computing it on
Pareto-optimal cycles (recall that $\xi\leq 1$ for classical stochastic SSHE
but PDHE can violate this bound [54, 55, 56, 57, 58, 59]). In Fig. 4a we show
a contour plot of $\xi$, computed with RL, as a function of $a$ and $c$.
Because of the mapping between SSHE and fast-Otto cycles [93], we have $\xi=1$
along the sides of the triangle, where only 2 objectives are optimized.
However, this mapping breaks down for finite-time cycles, allowing us to
observe a strong increase of $\xi$ in the green/purple region in Fig. 4a. As
shown in Fig. 2, this region corresponds to long cycles operated in the slow-
driving regime, where violations of thermodynamic uncertainty relations had
already been reported [56, 66]. In Figs. 4(b,c,d) we show a log-log plot of
$\xi$ respectively as a function of ${P}/P_{\text{max}}$, $\Delta{P}/\Delta
P(P_{\text{max}})$, and ${\Sigma}/\Sigma(P_{\text{max}})$ with the same color-
map as in Fig. 4a. We see that $\xi$ diverges in the limit of low power, low
fluctuations, and low entropy production as a power law. Indeed, using the
slow-driving approximation, we analytically prove that $\xi$ diverges as
$\expectationvalue*{{P}}^{-2}$, $\expectationvalue*{\Delta{P}}^{-1}$, and
$\expectationvalue*{\Sigma}^{-1}$. Such relations, plotted as black lines,
nicely agree with our RL results.

##### Conclusions.

We introduced a general framework to identify Pareto-optimal cycles between
power, efficiency, and power fluctuations in quantum or classical stochastic
heat engines, paving the way for their systematic optimization using optimal
control techniques such as the Pontryagin Minimum Principle [72] or
Reinforcement Learning [73]. As opposed to previous literature reviewed above,
we account for the crucial impact of power fluctuations. We then employed RL
to optimize a quantum dot-based heat engine, solving its exact finite-time and
out-of-equilibrium dynamics, providing us with new physical insights. We
observe an abrupt change in Pareto-optimal cycles when switching from the
optimization of $2$ objectives, where Otto cycles in the fast-driving regime
are optimal, to 3 objectives, where the optimal cycles have a finite period.
This feature, which shares analogies with the phase transition in protocol
space observed in Ref. [96], corresponds to a large enhancement of one of the
objectives while only slightly decreasing the other ones. Furthermore, we find
an exact mapping between Otto cycles in the fast-driving regime and SSHEs,
implying that a violation of the thermodynamic uncertainty relation ratio
$\xi$ in Eq. (1) requires the optimization of all $3$ objectives. We then find
that $\xi$ becomes arbitrarily large in the slow-driving regime. Cycles found
with RL display the best features analytically identified in the fast and slow
driving regimes.

##### Acknowledgments.

FN gratefully acknowledges funding by the BMBF (Berlin Institute for the
Foundations of Learning and Data – BIFOLD), the European Research Commission
(ERC CoG 772230) and the Berlin Mathematics Center MATH+ (AA1-6, AA2-8). PAE
gratefully acknowledges funding by the Berlin Mathematics Center MATH+ (AA1-6,
AA2-18). AR and MPL acknowledge funding from Swiss National Science Foundation
through an Ambizione Grant No. PZ00P2-186067. PA is supported by “la Caixa”
Foundation (ID 100010434, Grant No. LCF/BQ/DI19/11730023), and by the
Government of Spain (FIS2020-TRANQI and Severo Ochoa CEX2019-000910-S),
Fundacio Cellex, Fundacio Mir-Puig, Generalitat de Catalunya (CERCA, AGAUR SGR
1381).

## References

  * Giazotto  _et al._ [2006] F. Giazotto, T. T. Heikkilä, A. Luukanen, A. M. Savin, and J. P. Pekola, Opportunities for mesoscopics in thermometry and refrigeration: Physics and applications, Rev. Mod. Phys. 78, 217 (2006).
  * Pekola [2015] J. P. Pekola, Towards quantum thermodynamics in electronic circuits, Nat. Phys. 11, 118 (2015).
  * Vinjanampathy and Anders [2016] S. Vinjanampathy and J. Anders, Quantum thermodynamics, Contemp. Phys. 57, 545 (2016).
  * Benenti  _et al._ [2017] G. Benenti, G. Casati, K. Saito, and R. S.Whitney, Fundamental aspects of steady-state conversion of heat to work at the nanoscale, Phys. Rep. 694, 1 (2017).
  * Binder  _et al._ [2019] F. Binder, L. Correa, C. Gogolin, J. Anders, and G. Adesso,  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_ , (Springer International Publishing, Cham, Switzerland, 2019).
  * Geva and Kosloff [1992] E. Geva and R. Kosloff, A quantum-mechanical heat engine operating in finite time. a model consisting of spin-1/2 systems as the working fluid, J. Chem. Phys. 96, 3054 (1992).
  * Feldmann and Kosloff [2000] T. Feldmann and R. Kosloff, Performance of discrete heat engines and heat pumps in finite time, Phys. Rev. E 61, 4774 (2000).
  * Schmiedl and Seifert [2008] T. Schmiedl and U. Seifert, Efficiency at maximum power: An analytically solvable model for stochastic heat engines, Europhys. Lett. 81, 20003 (2008).
  * Esposito  _et al._ [2010a] M. Esposito, R. Kawai, K. Lindenberg, and C. V. den Broeck, Efficiency at maximum power of low-dissipation carnot engines, Phys. Rev. Lett. 105, 150603 (2010a).
  * Esposito  _et al._ [2010b] M. Esposito, R. Kawai, K. Lindenberg, and C. Van den Broeck, Quantum-dot carnot engine at maximum power, Phys. Rev. E 81, 041106 (2010b).
  * Allahverdyan  _et al._ [2013] A. E. Allahverdyan, K. V. Hovhannisyan, A. V. Melkikh, and S. G. Gevorkian, Carnot cycle at finite power: Attainability of maximal efficiency, Phys. Rev. Lett. 111, 050601 (2013).
  * Juergens  _et al._ [2013] S. Juergens, F. Haupt, M. Moskalets, and J. Splettstoesser, Thermoelectric performance of a driven double quantum dot, Phys. Rev. B 87, 245423 (2013).
  * Holubec and Ryabov [2015] V. Holubec and A. Ryabov, Efficiency at and near maximum power of low-dissipation heat engines, Phys. Rev. E 92, 052125 (2015).
  * Holubec and Ryabov [2016] V. Holubec and A. Ryabov, Maximum efficiency of low-dissipation heat engines at arbitrary power,  J. Stat. Mech. Theory Exp. 2016, 073204 (2016).
  * Karimi and Pekola [2016] B. Karimi and J. P. Pekola, Otto refrigerator based on a superconducting qubit: Classical and quantum performance, Phys. Rev. B 94, 184503 (2016).
  * Campisi and Fazio [2016] M. Campisi and R. Fazio, The power of a critical heat engine, Nat. Commun. 7, 11895 (2016).
  * [17] K. Yamamoto, O. Entin-Wohlman, A. Aharony, and N. Hatano, Efficiency bounds on thermoelectric transport in magnetic fields: The role of inelastic processes, Phys. Rev. B 94, 121402(R) (2016).
  * Kosloff and Rezek [2017] R. Kosloff and Y. Rezek, The quantum harmonic otto cycle, Entropy 19, 136 (2017).
  * [19] V. Patel, V. Savsani, and A. Mudgal, Many-objective thermodynamic optimization of Stirling heat engine, Energy 125, 629 (2017).
  * Cavina  _et al._ [2018] V. Cavina, A. Mari, A. Carlini, and V. Giovannetti, Optimal thermodynamic control in open quantum systems, Phys. Rev. A 98, 012139 (2018).
  * Abiuso and Giovannetti [2019] P. Abiuso and V. Giovannetti, Non-markov enhancement of maximum power for quantum thermal machines, Phys. Rev. A 99, 052106 (2019).
  * Ma  _et al._ [2018] Y.-H. Ma, D. Xu, H. Dong, and C.-P. Sun, Universal constraint for efficiency and power of a low-dissipation heat engine, Phys. Rev. E 98, 042112 (2018).
  * Erdman  _et al._ [2019] P. A. Erdman, V. Cavina, R. Fazio, F. Taddei, and V. Giovannetti, Maximum power and corresponding efficiency for two-level heat engines and refrigerators: optimality of fast cycles, New J. Phys. 21, 103049 (2019).
  * Menczel  _et al._ [2019] P. Menczel, T. Pyhäranta, C. Flindt, and K. Brandner, Two-stroke optimization scheme for mesoscopic refrigerators, Phys. Rev. B 99, 224306 (2019).
  * Chen  _et al._ [2019] J. Chen, C. Sun, and H. Dong, Boosting the performance of quantum otto heat engines, Phys. Rev. E 100, 032144 (2019).
  * Abiuso and Perarnau-Llobet [2020] P. Abiuso and M. Perarnau-Llobet, Optimal cycles for low-dissipation heat engines, Phys. Rev. Lett. 124, 110606 (2020).
  * Cavina  _et al._ [2021] V. Cavina, P. A. Erdman, P. Abiuso, L. Tolomeo, and V. Giovannetti, Maximum-power heat engines and refrigerators in the fast-driving regime, Phys. Rev. A 104, 032226 (2021).
  * [28] R. A. Vargas-Hernández, R. T. Q. Chen, K. A. Jung, and P. Brumer, Fully differentiable optimization protocols for non-equilibrium steady states, New J. Phys. 23, 123006 (2021).
  * Terrén Alonso  _et al._ [2022] P. Terrén Alonso, P. Abiuso, M. Perarnau-Llobet, and L. Arrachea, Geometric optimization of nonequilibrium adiabatic thermal machines and implementation in a qubit system, PRX Quantum 3, 010326 (2022).
  * Ye  _et al._ [2022] Z. Ye, F. Cerisola, P. Abiuso, J. Anders, M. Perarnau-Llobet, and V. Holubec, Optimal finite-time heat engines under constrained control Phys. Rev. Research 4, 043130 (2022).
  * Holubec [2014] V. Holubec, An exactly solvable model of a stochastic heat engine: optimization of power, power fluctuations and efficiency, J. Stat. Mech. Theory Exp. 2014, P05022 (2014).
  * Miller and Mehboudi [2020] H. J. Miller and M. Mehboudi, Geometry of work fluctuations versus efficiency in microscopic thermal machines, Phys. Rev. Lett. 125, 260602 (2020).
  * Denzler and Lutz [2021] T. Denzler and E. Lutz, Power fluctuations in a finite-time quantum carnot engine, Phys. Rev. Research 3, L032041 (2021).
  * Saryal and Agarwalla [2021] S. Saryal and B. K. Agarwalla, Bounds on fluctuations for finite-time quantum otto cycle, Phys. Rev. E 103, L060103 (2021).
  * Mehboudi and Miller [2021] M. Mehboudi and H. J. D. Miller, Thermodynamic length and work optimisation for gaussian quantum states, Phys. Rev. A 105, 062434 (2022).
  * Barato and Seifert [2015] A. C. Barato and U. Seifert, Thermodynamic uncertainty relation for biomolecular processes, Phys. Rev. Lett. 114, 158101 (2015).
  * Gingrich  _et al._ [2016] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, Dissipation bounds all steady-state current fluctuations, Phys. Rev. Lett. 116, 120601 (2016).
  * Pietzonka and Seifert [2018] P. Pietzonka and U. Seifert, Universal trade-off between power, efficiency, and constancy in steady-state heat engines, Phys. Rev. Lett. 120, 190602 (2018).
  * Guarnieri  _et al._ [2019] G. Guarnieri, G. T. Landi, S. R. Clark, and J. Goold, Thermodynamics of precision in quantum nonequilibrium steady states, Phys. Rev. Research 1, 033021 (2019).
  * Timpanaro  _et al._ [2019] A. M. Timpanaro, G. Guarnieri, J. Goold, and G. T. Landi, Thermodynamic uncertainty relations from exchange fluctuation theorems, Phys. Rev. Lett. 123, 090604 (2019).
  * Falasco  _et al._ [2020] G. Falasco, M. Esposito, and J.-C. Delvenne, Unifying thermodynamic uncertainty relations, New J. Phys. 22, 053046 (2020).
  * Friedman  _et al._ [2020] H. M. Friedman, B. K. Agarwalla, O. Shein-Lumbroso, O. Tal, and D. Segal, Thermodynamic uncertainty relation in atomic-scale quantum conductors, Phys. Rev. B 101, 195423 (2020).
  * Horowitz and Gingrich [2020] J. M. Horowitz and T. R. Gingrich, Thermodynamic uncertainty relations constrain non-equilibrium fluctuations, Nat. Phys. 16, 15 (2020).
  * Saryal  _et al._ [2021] S. Saryal, M. Gerry, I. Khait, D. Segal, and B. K. Agarwalla, Universal bounds on fluctuations in continuous thermal machines, Phys. Rev. Lett. 127, 190603 (2021).
  * Agarwalla and Segal [2018] B. K. Agarwalla and D. Segal, Assessing the validity of the thermodynamic uncertainty relation in quantum systems, Phys. Rev. B 98, 155438 (2018).
  * Brandner  _et al._ [2018] K. Brandner, T. Hanazato, and K. Saito, Thermodynamic bounds on precision in ballistic multiterminal transport, Phys. Rev. Lett. 120, 090601 (2018).
  * Ptaszyński [2018] K. Ptaszyński, Coherence-enhanced constancy of a quantum thermoelectric generator, Phys. Rev. B 98, 085425 (2018).
  * Carollo  _et al._ [2019] F. Carollo, R. L. Jack, and J. P. Garrahan, Unraveling the large deviation statistics of markovian open quantum systems, Phys. Rev. Lett. 122, 130605 (2019).
  * Liu and Segal [2019] J. Liu and D. Segal, Thermodynamic uncertainty relation in quantum thermoelectric junctions, Phys. Rev. E 99, 062141 (2019).
  * Kalaee  _et al._ [2021] A. A. S. Kalaee, A. Wacker, and P. P. Potts, Violating the thermodynamic uncertainty relation in the three-level maser, Phys. Rev. E 104, L012103 (2021).
  * Rignon-Bret  _et al._ [2021] A. Rignon-Bret, G. Guarnieri, J. Goold, and M. T. Mitchison, Thermodynamics of precision in quantum nanomachines, Phys. Rev. E 103, 012133 (2021).
  * Ehrlich and Schaller [2021] T. Ehrlich and G. Schaller, Broadband frequency filters with quantum dot chains, Phys. Rev. B 104, 045424 (2021).
  * Gerry and Segal [2022] M. Gerry and D. Segal, Absence and recovery of cost-precision tradeoff relations in quantum transport, Phys. Rev. B 105, 155401 (2022).
  * Barato and Seifert [2016] A. C. Barato and U. Seifert, Cost and precision of brownian clocks, Phys. Rev. X 6, 041053 (2016).
  * Proesmans and den Broeck [2017] K. Proesmans and C. V. den Broeck, Discrete-time thermodynamic uncertainty relation, EPL 119, 20001 (2017).
  * Holubec and Ryabov [2018] V. Holubec and A. Ryabov, Cycling tames power fluctuations near optimum efficiency, Phys. Rev. Lett. 121, 120601 (2018).
  * Cangemi  _et al._ [2020] L. M. Cangemi, V. Cataudella, G. Benenti, M. Sassetti, and G. De Filippis, Violation of thermodynamics uncertainty relations in a periodically driven work-to-work converter from weak to strong dissipation, Phys. Rev. B 102, 165418 (2020).
  * Menczel  _et al._ [2021] P. Menczel, E. Loisa, K. Brandner, and C. Flindt, Thermodynamic uncertainty relations for coherently driven open quantum systems, J. Phys. A: Math. Theor. 54, 314002 (2021).
  * Lu  _et al._ [2022] J. Lu, Z. Wang, J. Peng, C. Wang, J.-H. Jiang, and J. Ren, Geometric thermodynamic uncertainty relation in a periodically driven thermoelectric heat engine, Phys. Rev. B 105, 115428 (2022).
  * Koyuk  _et al._ [2019] T. Koyuk, U. Seifert, and P. Pietzonka, A generalization of the thermodynamic uncertainty relation to periodically driven systems, J. Phys. A: Math. Theor. 52, 02LT02 (2019).
  * Koyuk and Seifert [2019] T. Koyuk and U. Seifert, Operationally accessible bounds on fluctuations and entropy production in periodically driven systems, Phys. Rev. Lett. 122, 230601 (2019).
  * Koyuk and Seifert [2020] T. Koyuk and U. Seifert, Thermodynamic uncertainty relation for time-dependent driving, Phys. Rev. Lett. 125, 260604 (2020).
  * Van Vu and Hasegawa [2020] T. Van Vu and Y. Hasegawa, Thermodynamic uncertainty relations under arbitrary control protocols, Phys. Rev. Research 2, 013060 (2020).
  * Potanina  _et al._ [2021] E. Potanina, C. Flindt, M. Moskalets, and K. Brandner, Thermodynamic bounds on coherent transport in periodically driven conductors, Phys. Rev. X 11, 021013 (2021).
  * Hasegawa [2021] Y. Hasegawa, Irreversibility, loschmidt echo, and thermodynamic uncertainty relation, Phys. Rev. Lett. 127, 240602 (2021).
  * Miller  _et al._ [2021a] H. J. D. Miller, M. H. Mohammady, M. Perarnau-Llobet, and G. Guarnieri, Thermodynamic uncertainty relation in slowly driven quantum heat engines, Phys. Rev. Lett. 126, 210603 (2021a).
  * Miller  _et al._ [2021b] H. J. D. Miller, M. H. Mohammady, M. Perarnau-Llobet, and G. Guarnieri, Joint statistics of work and entropy production along quantum trajectories, Phys. Rev. E 103, 052138 (2021b).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of N‐level systems, J. Math. Phys. 17, 821 (1976).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys 48, 119 (1976).
  * Breuer and Petruccione [2002] H. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, New York, 2002).
  * Yamaguchi  _et al._ [2017] M. Yamaguchi, T. Yuge, and T. Ogawa, Markovian quantum master equation beyond adiabatic regime, Phys. Rev. E 95, 012136 (2017).
  * Kirk [2004] D. E. Kirk,  _Optimal Control Theory: an Introduction_ (Dover, 2004)..
  * Sutton and Barto [2018] R. S. Sutton and A. G. Barto,  _Reinforcement learning: An introduction_ (MIT press, 2018).
  * Sothmann  _et al._ [2015] B. Sothmann, R. Sánchez, and A. N. Jordan, Thermoelectric energy harvesting with quantum dots, Nanotechnology 26, 032001 (2015).
  * Pekola  _et al._ [2019] J. P. Pekola, B. Karimi, G. Thomas, and D. V. Averin, Supremacy of incoherent sudden cycles, Phys. Rev. B 100, 085405 (2019).
  * Cangemi  _et al._ [2021] L. M. Cangemi, M. Carrega, A. De Candia, V. Cataudella, G. De Filippis, M. Sassetti, and G. Benenti, Optimal energy conversion through antiadiabatic driving breaking time-reversal symmetry, Phys. Rev. Research 3, 013237 (2021).
  * Salamon and Berry [1983] P. Salamon and R. S. Berry, Thermodynamic length and dissipated availability, Phys. Rev. Lett. 51, 1127 (1983).
  * Wang  _et al._ [2011] J. Wang, J. He, and X. He, Performance analysis of a two-state quantum heat engine working with a single-mode radiation field in a cavity, Phys. Rev. E 84, 041127 (2011).
  * Avron  _et al._ [2012] J. E. Avron, M. Fraas, G. M. Graf, and P. Grech, Adiabatic theorems for generators of contracting evolutions, Commun. Math. Phys. 314, 163 (2012).
  * Ludovico  _et al._ [2016] M. F. Ludovico, F. Battista, F. von Oppen, and L. Arrachea, Adiabatic response and quantum thermoelectrics for ac-driven quantum systems, Phys. Rev. B 93, 075136 (2016).
  * Cavina  _et al._ [2017] V. Cavina, A. Mari, and V. Giovannetti, Slow dynamics and thermodynamics of open quantum systems, Phys. Rev. Lett. 119, 050601 (2017).
  * Miller  _et al._ [2019] H. J. D. Miller, M. Scandi, J. Anders, and M. Perarnau-Llobet, Work fluctuations in slow processes: Quantum signatures and optimal control, Phys. Rev. Lett. 123, 230603 (2019).
  * Scandi and Perarnau-Llobet [2019] M. Scandi and M. Perarnau-Llobet, Thermodynamic length in open quantum systems, Quantum 3, 197 (2019).
  * Brandner and Saito [2020] K. Brandner and K. Saito, Thermodynamic geometry of microscopic heat engines, Phys. Rev. Lett. 124, 040602 (2020).
  * Bhandari  _et al._ [2020] B. Bhandari, P. T. Alonso, F. Taddei, F. von Oppen, R. Fazio, and L. Arrachea, Geometric properties of adiabatic quantum thermal machines, Phys. Rev. B 102, 155407 (2020).
  * Abiuso  _et al._ [2020] P. Abiuso, H. J. D. Miller, M. Perarnau-Llobet, and M. Scandi, Geometric optimisation of quantum thermodynamic processes, Entropy 22 (2020).
  * Frim and DeWeese [2021] A. G. Frim and M. R. DeWeese, Geometric bound on the efficiency of irreversible thermodynamic cycles, Phys. Rev. Lett. 128, 230601 (2022).
  * Eglinton and Brandner [2022] J. Eglinton and K. Brandner, Geometric bounds on the power of adiabatic thermal machines, Phys. Rev. E 105, L052102 (2022).
  * Brandner  _et al._ [2015] K. Brandner, K. Saito, and U. Seifert, Thermodynamics of micro- and nano-systems driven by periodic temperature variations, Phys. Rev. X 5, 031019 (2015).
  * Brandner and Seifert [2016] K. Brandner and U. Seifert, Periodic thermodynamics of open quantum systems, Phys. Rev. E 93, 062134 (2016).
  * Movilla Miangolarra  _et al._ [2021a] O. Movilla Miangolarra, R. Fu, A. Taghvaei, Y. Chen, and T. T. Georgiou, Underdamped stochastic thermodynamic engines in contact with a heat bath with arbitrary temperature profile, Phys. Rev. E 103, 062103 (2021a).
  * Movilla Miangolarra  _et al._ [2021b] O. Movilla Miangolarra, A. Taghvaei, R. Fu, Y. Chen, and T. T. Georgiou, Energy harvesting from anisotropic fluctuations, Phys. Rev. E 104, 044101 (2021b).
  * [93] See Supplemental Material for details on the the general framework, the reinforcement learning method, and the fast and slow driving regimes.
  * Alicki [1979] R. Alicki, The quantum open system as a model of the heat engine, J. Phys. A: Math. Gen. 12, L103 (1979).
  * [95] L. F. Seoane and R. Solé, Multiobjective optimiza- tion and phase transitions, in Proceedings of ECCS 2014, edited by S. Battiston, F. De Pellegrini, G. Caldarelli, and E. Merelli (Springer International Publishing, Cham, 2016) pp. 259–270.
  * Solon and Horowitz [2018] A. P. Solon and J. M. Horowitz, Phase transition in protocols minimizing work fluctuations, Phys. Rev. Lett. 120, 180605 (2018).
  * Haarnoja  _et al._ [2018a] T. Haarnoja, A. Zhou, K. Hartikainen, G. Tucker, S. Ha, J. Tan, V. Kumar, H. Zhu, A. Gupta, P. Abbeel, _et al._ , Soft actor-critic algorithms and applications, arXiv:1812.05905 (2018a).
  * Haarnoja  _et al._ [2018b] T. Haarnoja, S. Ha, A. Zhou, J. Tan, G. Tucker, and S. Levine, Learning to walk via deep reinforcement learning, arXiv:1912.11077 (2018b).
  * Haarnoja  _et al._ [2018c] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor, in  _International Conference on Machine Learning_, Vol. 80 (PMLR, 2018) p. 1861.
  * Erdman and Noé [2022a] P. A. Erdman and F. Noé, Identifying optimal cycles in quantum thermal machines with reinforcement-learning, NPJ Quantum Inf. 8, 1 (2022a).
  * Erdman and Noé [2022b] P. A. Erdman and F. Noé, Driving black-box quantum thermal machines with optimal power/efficiency trade-offs using reinforcement learning, arXiv:2204.04785 (2022b).
  * Sgroi  _et al._ [2021] P. Sgroi, G. M. Palma, and M. Paternostro, Reinforcement learning approach to nonequilibrium quantum thermodynamics,  Phys. Rev. Lett. 126, 020601 (2021).
  * Khait  _et al._ [2022] I. Khait, J. Carrasquilla, and D. Segal, Optimal control of quantum thermal machines using machine learning, Phys. Rev. Research 4, L012029 (2022).
  * [104] Y. Ashida and T. Sagawa, Learning the best nanoscale heat engines through evolving network topology, Commun. Phys. 4, 45 (2021).
  * Bukov  _et al._ [2018] M. Bukov, A. G. R. Day, D. Sels, P. Weinberg, A. Polkovnikov, and P. Mehta, Reinforcement learning in different phases of quantum control, Phys. Rev. X 8, 031086 (2018).
  * Fösel  _et al._ [2018] T. Fösel, P. Tighineanu, T. Weiss, and F. Marquardt, Reinforcement learning with neural networks for quantum feedback, Phys. Rev. X 8, 031084 (2018).
  * An and Zhou [2019] Z. An and D. Zhou, Deep reinforcement learning for quantum gate control, EPL 126, 60002 (2019).
  * Niu  _et al._ [2019] M. Y. Niu, S. Boixo, V. N. Smelyanskiy, and H. Neven, Universal quantum control through deep reinforcement learning, NPJ Quantum Inf. 5, 33 (2019).
  * Zhang  _et al._ [2019] X.-M. Zhang, Z. Wei, R. Asad, X.-C. Yang, and X. Wang, When does reinforcement learning stand out in quantum control? a comparative study on state preparation, NPJ Quantum Inf. 5, 85 (2019).
  * Dalgaard  _et al._ [2020] M. Dalgaard, F. Motzoi, J. J. Sørensen, and J. Sherson, Global optimization of quantum dynamics with alphazero deep exploration, NPJ Quantum Inf. 6, 6 (2020).
  * Mackeprang  _et al._ [2020] J. Mackeprang, D. B. R. Dasari, and J. Wrachtrup, A reinforcement learning approach for quantum state engineering, Quantum Mach. Intell. 2, 5 (2020).
  * Schäfer  _et al._ [2020] F. Schäfer, M. Kloc, C. Bruder, and N. Lörch, A differentiable programming method for quantum control, Mach. Learn.: Sci. Technol. 1, 035009 (2020).
  * Sweke  _et al._ [2020] R. Sweke, M. S. Kesselring, E. P. L. van Nieuwenburg, and J. Eisert, Reinforcement learning decoders for fault-tolerant quantum computation, Mach. Learn.: Sci. Technol. 2, 025005 (2020).
  * Schäfer  _et al._ [2021] F. Schäfer, P. Sekatski, M. Koppenhöfer, C. Bruder, and M. Kloc, Control of stochastic quantum dynamics by differentiable programming, Mach. Learn.: Sci. Technol. 2, 035004 (2021).
  * Porotti  _et al._ [2021] R. Porotti, A. Essig, B. Huard, and F. Marquardt, Deep reinforcement learning for quantum state preparation with weak nonlinear measurements, Quantum 6, 747 (2022).
  * Metz and Bukov [2022] F. Metz and M. Bukov, Self-correcting quantum many-body control using reinforcement learning with tensor networks, arXiv:2201.11790 (2022).
  * Coopmans  _et al._ [2022] L. Coopmans, S. Campbell, G. D. Chiara, and A. Kiely, Optimal control in disordered quantum systems, Phys. Rev. Research 4, 043138 (2022).
  * Menczel and Brandner [2019] P. Menczel and K. Brandner, Limit cycles in periodically driven open quantum systems, J. Phys. A Math. Theor. 52, 43LT01 (2019).
  * Christodoulou [2019] P. Christodoulou, Soft actor-critic for discrete action settings, arXiv:1910.07207 (2019).
  * Delalleau  _et al._ [2019] O. Delalleau, M. Peter, E. Alonso, and A. Logut, Discrete and continuous action representation for practical rl in video games, arXiv:1912.11077 (2019).
  * Mnih  _et al._ [2015] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, _et al._ , Human-level control through deep reinforcement learning, Nature 518, 529 (2015).
  * Silver  _et al._ [2017] D. Silver, J. Schrittwieser, K. Simonyan, I. Antonoglou, A. Huang, A. Guez, T. Hubert, L. Baker, M. Lai, A. Bolton, _et al._ , Mastering the game of go without human knowledge, Nature 550, 354 (2017).
  * Vinyals  _et al._ [2019] O. Vinyals, I. Babuschkin, W. M. Czarnecki, M. Mathieu, A. Dudzik, J. Chung, D. H. Choi, R. Powell, T. Ewalds, P. Georgiev, _et al._ , Grandmaster level in starcraft ii using multi-agent reinforcement learning, Nature 575, 350 (2019).
  * [124] J. Degrave, F. Felici, J. Buchli, M. Neunert, B. Tracey, F. Carpanese, T. Ewalds, R. Hafner, A. Abdolmaleki, D. de las Casas, _et al._ , Magnetic control of tokamak plasmas through deep reinforcement learning, Nature 602, 414 (2022)..
  * Schlögl [1985] F. Schlögl, Thermodynamic metric and stochastic measures, Zeitschrift für Physik B Condensed Matter 59, 449 (1985).
  * Nulton  _et al._ [1985] J. Nulton, P. Salamon, B. Andresen, and Q. Anmin, Quasistatic processes as step equilibrations, J. Chem. Phys. 83, 334 (1985) .
  * Crooks [2007] G. E. Crooks, Measuring thermodynamic length, Phys. Rev. Lett. 99, 100602 (2007).
  * Sivak and Crooks [2012] D. A. Sivak and G. E. Crooks, Thermodynamic metrics and optimal paths, Phys. Rev. Lett. 108, 190602 (2012).
  * Bonança and Deffner [2014] M. V. S. Bonança and S. Deffner, Optimal driving of isothermal processes close to equilibrium, J. Chem. Phys. 140, 244119 (2014) .
  * Deffner and Bonança [2020] S. Deffner and M. V. S. Bonança, Thermodynamic control —an old paradigm with new applications, EPL 131, 20001 (2020).
  * Weber [1956] J. Weber, Fluctuation dissipation theorem, Phys. Rev. 101, 1620 (1956).
  * Kubo [1966] R. Kubo, The fluctuation-dissipation theorem, Rep. Prog. Phys. 29, 255 (1966).
  * Hermans [1991] J. Hermans, Simple analysis of noise and hysteresis in (slow-growth) free energy simulations, J. Phys. Chem. 95, 9029 (1991).
  * Jarzynksi [1997] C. Jarzynksi, A nonequilibrium equality for free energy differences, Phys. Rev. Lett. 78, 2690 (1997).
  * Hendrix and Jarzynski [2001] D. A. Hendrix and C. Jarzynski, A “fast growth” method of computing free energy differences, J. Chem. Phys. 114, 5974 (2001) .
  * Jenčová [2004] A. Jenčová, Geodesic distances on density matrices, J. Math. Phys. 45, 1787 (2004).

Supplemental Material for “Pareto-optimal cycles for power, efficiency and
fluctuations of quantum heat engines using reinforcement learning”

##  I General Framework and Reinforcement Learning approach

In this section we describe a general framework to compute the power, entropy
production and power fluctuations of a quantum thermal machine. We then show
how it allows us to use RL to optimize the trade-off between these three
objectives. Assuming that the dynamics is Markovian, the reduced state evolves
according to Eq. (2) of the main text, i.e.

$\dot{\rho}_{t}=\mathcal{L}_{\vec{u}(t),\beta(t)}\rho_{t}.$ (S1)

As usual, we employ the first law to split the internal energy
$U(t)=\Tr[\rho_{t}H_{\vec{u}(t)}]$ of the system into absorbed heat and output
work [94]:

$\displaystyle dU=dQ-dW\;,\quad dQ=\Tr[d\rho_{t}H_{\vec{u}(t)}]\;,\quad
dW=-\Tr[\rho_{t}dH_{\vec{u}(t)}]\;.$ (S2)

The instantaneous power _output_ is thus given by Eq. (3) of the main text,
i.e.

$\displaystyle P(t)=-\Tr[\rho_{t}\dot{H}_{\vec{u}(t)}]$ (S3)

and, according to Clausius theorem for cycles (in which we set $\Delta S=0$ of
the quantum system since we consider a periodic driving in the limiting
cycle), the instantaneous irreversible entropy production is given by Eq. (3)
of the main text, i.e.

$\displaystyle\Sigma(t)=-\beta(t)\Tr[\dot{\rho}_{t}H_{\vec{u}(t)}]\;.$ (S4)

For what concerns power fluctuations, let $\langle\sigma^{2}\rangle_{t}$ be
the work fluctuations between time $0$ at $t$. Using Eq. (2) of Ref. [82], we
can write the fluctuations as

$\langle\sigma^{2}\rangle_{t}=\int_{0}^{t}dt_{1}\Tr[\dot{H}_{\vec{u}(t_{1})}s_{t_{1}}],$
(S5)

where we introduce

$s_{t_{1}}\equiv\int_{0}^{t_{1}}dt_{2}P(t_{1},t_{2})[\dot{H}_{\vec{u}(t_{2})}\rho_{t_{2}}]+\expectationvalue{w}_{t_{1}}\rho_{t_{1}}+h.c.\;.$
(S6)

Here
$P(t_{1},t_{2})\equiv\overleftarrow{T}\exp[\int_{t_{2}}^{t_{1}}d\tau\mathcal{L}_{\vec{u}(\tau),\beta(\tau)}]$
is the propagator,
$\expectationvalue{w}_{t_{1}}\equiv-\int_{0}^{t_{1}}dt_{2}\Tr[\rho_{t_{2}}\dot{H}_{\vec{u}(t_{2})}]$
is the total average work extracted between time $0$ and $t_{1}$, $h.c.$
represents the complex conjugate of the whole right-hand-side, and in deriving
Eq. (S5) we used that $[P(t_{1},t_{2})O]^{\dagger}=P(t_{1},t_{2})O^{\dagger}$,
which stems from
$[\mathcal{L}_{\vec{u}(\tau),\beta(\tau)}O]^{\dagger}=\mathcal{L}_{\vec{u}(\tau),\beta(\tau)}O^{\dagger}$.
Taking the time-derivative of $s_{t}$, we find that it satisfies the equation
of motion given in Eq. (6) of the main text, i.e.

$\dot{s}_{t}=\mathcal{L}_{\vec{u}(t),\beta(t)}s_{t}+\\{\rho,\dot{H}_{\vec{u}(t)}\\}-2\Tr[\rho_{t}\dot{H}_{\vec{u}(t)}]\rho_{t},$
(S7)

where $\\{\cdot,\cdot\\}$ denotes the anti-commutator. Using the definition of
$s_{t}$ given in Eq. (S6), we see that $s_{t}$ is an Hermitian operator, and
that $\Tr[s_{t}]=0$. This last property stems from the fact that
$P(t_{1},t_{2})[\cdot]$ is trace-preserving.

###  I.1 “Extended state” and formulation of the optimization problem

We now define $(\rho_{t},s_{t})$ as an “extended state”. We have just shown
that the equation of motion of the extended state is a set of first-order
differential equations [see Eqs. (S1) and (S7)]. Furthermore, the average
power $\expectationvalue{P}\equiv\expectationvalue{w}_{\tau}/\tau$, the
average fluctuations $\expectationvalue{\Delta
P}\equiv\langle\sigma^{2}\rangle_{\tau}/\tau$ and the average entropy
production $\expectationvalue{\Sigma}\equiv(1/\tau)\int_{0}^{\tau}\Sigma(t)dt$
can be expressed as a time-average of the extended state and of the controls
$\vec{u}(t)$ and $\beta(t)$, i.e.

$\displaystyle\expectationvalue{P}$
$\displaystyle=-\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[\rho_{t}\dot{H}_{\vec{u}(t)}]=\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[\dot{\rho}_{t}H_{\vec{u}(t)}],$
(S8) $\displaystyle\expectationvalue{\Delta P}$
$\displaystyle=\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[s_{t}\dot{H}_{\vec{u}(t)}],$
$\displaystyle\expectationvalue{\Sigma}$
$\displaystyle=-\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[\dot{\rho}_{t}H_{\vec{u}(t)}]\beta(t).$

Therefore, also the figure of merit $\expectationvalue{F}$ in Eq. (8) of the
main text can be expressed as a time integral of a function of the extended
state $(\rho_{t},s_{t})$ and of the controls $\vec{u}(t)$ and $\beta(t)$, i.e.

$\expectationvalue{F}=\int_{0}^{\tau}G\left(\rho_{t},s_{t},\vec{u}(t),\beta(t)\right)dt,$
(S9)

where $G(\rho_{t},s_{t},\vec{u}(t),\beta(t))$ is a suitable function. This is
precisely the type of problem that can readily tackled using optimal control
techniques such as Pontryagin Minimum Principle [72] or Reinforcement Learning
[73].

As a final remark, one can study heat engines in the limiting cycle imposing
periodic boundary conditions for the extended state, i.e.
$(\rho_{0},s_{0})=(\rho_{\tau},s_{\tau})$. Indeed, we are interested in
studying thermodynamic cycles, i.e. when both the control and the quantum
state $\rho_{t}$ are periodic with period $\tau$. Assuming that the
Lindbladian $\mathcal{L}_{\vec{u}(t),\beta(t)}$ has a single fixed point, it
can be shown that a periodic control eventually drives the state $\rho_{t}$
towards the limiting cycle solution, where $\rho_{t}$ becomes periodic with
the same period $\tau$ of the driving [118]. Furthermore, the equation of
motion of $s_{t}$ has the same homogeneous part
($\mathcal{L}_{\vec{u}(t),\beta(t)}[s_{t}]$) as $\rho_{t}$, and the non-
homogeneous term becomes periodic once $\rho_{t}$ reaches the limiting cycle.
Therefore, also $s_{t}$ naturally reaches a limiting cycle solution where it
is periodic with the same period $\tau$.

###  I.2 Optimizing the entropy production instead of the efficiency

Here we discuss the relation between minimizing the entropy production or
maximizing the efficiency. We generalize the discussion of Ref. [101], where
power fluctuations are not considered. We start by noticing that we can
express the efficiency of a heat engine in terms of the average power and
entropy production, i.e.

$\eta=\eta_{\text{c}}\,[1+\expectationvalue{\Sigma}/(\beta_{\text{C}}\expectationvalue*{P})]^{-1}.$
(S10)

We show that, thanks to this dependence of $\eta$ on $\expectationvalue*{P}$
and $\expectationvalue*{\Sigma}$, optimizing a trade-off between high power,
low fluctuations and high efficiency yields all the Pareto optimal trade-offs
between high power, low fluctuations, and low entropy-production up to a
change of the weights $(a,b,c)$.

First, we provide a non-rigorous argument as follows: notice that a point
$\\{\expectationvalue{P},\expectationvalue{\Delta
P},\expectationvalue{\Sigma}\\}$ belongs to the Pareto front iff there exists
no cycle outperforming it in all 3 quantities, i.e. larger power, smaller
fluctuations and smaller entropy production. Due to the functional dependence
of $\eta(\expectationvalue{\Sigma},\expectationvalue{P})$ (S10), the
corresponding point $\\{\expectationvalue{P},\expectationvalue{\Delta
P},\eta(\expectationvalue{\Sigma},\expectationvalue{P})\\}$ will belong to the
power-fluctuations-efficiency Pareto front. In fact, if the point
$\\{\expectationvalue{P},\expectationvalue{\Delta
P},\eta(\expectationvalue{\Sigma},\expectationvalue{P})\\}$ is not on the
Pareto front, then a cycle
$\\{\expectationvalue{P}^{\prime},\expectationvalue{\Delta
P}^{\prime},\eta^{\prime}\\}$ exists such that one of the following statements
holds

$\displaystyle\text{(Case
1)}\quad\expectationvalue{P}^{\prime}>\expectationvalue{P}\;,\quad\expectationvalue{\Delta
P}^{\prime}=\expectationvalue{\Delta
P}\;,\quad\eta^{\prime}=\eta(\expectationvalue{\Sigma},\expectationvalue{P})\;;$
(S11) $\displaystyle\text{(Case
2)}\quad\expectationvalue{P}^{\prime}=\expectationvalue{P}\;,\quad\expectationvalue{\Delta
P}^{\prime}<\expectationvalue{\Delta
P}\;,\quad\eta^{\prime}=\eta(\expectationvalue{\Sigma},\expectationvalue{P})\;;$
(S12) $\displaystyle\text{(Case
3)}\quad\expectationvalue{P}^{\prime}=\expectationvalue{P}\;,\quad\expectationvalue{\Delta
P}^{\prime}=\expectationvalue{\Delta
P}\;,\quad\eta^{\prime}>\eta(\expectationvalue{\Sigma},\expectationvalue{P})\;.$
(S13)

Notice that in (Case 2) the corresponding
$\\{\expectationvalue{P}^{\prime},\expectationvalue{\Delta
P}^{\prime},\expectationvalue{\Sigma}(\expectationvalue{P}^{\prime},\eta^{\prime})\\}\equiv\\{\expectationvalue{P},\expectationvalue{\Delta
P}^{\prime},\expectationvalue{\Sigma}\\}$ clearly violates the power-
fluctuations-entropy production Pareto front, i.e. outperforms the point
$\\{\expectationvalue{P},\expectationvalue{\Delta
P},\expectationvalue{\Sigma}\\}$ because of smaller fluctuations. Similarly,
by inverting (S10)

$\Sigma(P,\eta)=\frac{\eta_{\text{c}}-\eta}{\eta}\beta_{\text{C}}P$ (S14)

we see that in (Case 3) the corresponding point violates the Pareto front via
a reduction of $\expectationvalue{\Sigma}$. Finally (Case 1) is non trivial,
but it is intuitive (and empirically verified) that in power-fluctuations-
efficiency trade-offs it is always possible to reduce the power in favor of
efficiency and fluctuations, (typically by slowing down the whole cycle, i.e.
increasing $\tau$). This means that the occurrence of (Case 1) induces the
existence of (Case 2) and (Case 3), closing the argument.

Besides the above intuitive argument, to prove such equivalence
mathematically, we can prove that the cycles that maximize

$\expectationvalue*{G}\equiv
a\frac{\expectationvalue*{P}}{P_{\text{max}}}-b\frac{\expectationvalue*{\Delta
P}}{\Delta P(P_{\text{max}})}+c\,\frac{\eta}{\eta_{\text{c}}}$ (S15)

for some values of the weights $a\geq 0$, $b\geq 0$, $c\geq 0$ such that
$a+b+c=1$, also maximize the figure of merit in Eq. (8) of the main text for
some (possibly different) non-negative values of the weights summing to one.
To simplify the proof and the notation, we consider the following two
functions

$\displaystyle F(a,b,c;\theta)$ $\displaystyle=aP(\theta)-b\Delta
P(\theta)-c\Sigma(P(\theta),\eta(\theta)),$ (S16) $\displaystyle
G(a,b,c;\theta)$ $\displaystyle=aP(\theta)-b\Delta P(\theta)+c\eta(\theta),$

where $P(\theta)$, $\Delta P(\theta)$, and $\eta(\theta)$ represent the
average power, fluctuations and efficiency of a cycle parameterized by a set
of parameters $\theta$, and $\Sigma(P,\eta)$ is given by Eq. (S14).

We wish to prove the following. Given three fixed scalars
$a_{1},b_{1},c_{1}>0$, that do not necessarily sum to 1, let $\theta_{1}$ be
the value of $\theta$ that locally maximizes $G(a_{1},b_{1},c_{1};\theta)$.
Then, it is always possible to identify three positive scalars
$a_{2},b_{2},c_{2}>0$, such that the same parameters $\theta_{1}$ (i.e. the
same cycle) is a local maximum for $F(a_{2},b_{2},c_{2};\theta)$. In the
following, we will use that

$\displaystyle\partial_{P}\Sigma$ $\displaystyle\geq 0,$
$\displaystyle\partial_{\eta}\Sigma$ $\displaystyle<0,$ (S17)

and that the Hessian $H^{(\Sigma)}$ of $\Sigma(P,\eta)$ is given by

$H^{(\Sigma)}=\begin{pmatrix}0&-\beta_{\text{C}}\frac{\eta_{\text{c}}}{\eta^{2}}\\\
-\beta_{\text{C}}\frac{\eta_{\text{c}}}{\eta^{2}}&2\beta_{\text{C}}P\frac{\eta_{\text{c}}}{\eta^{3}}\end{pmatrix}.$
(S18)

Proof: by assumption, $\theta_{1}$ is a local maximum for
$G(a_{1},b_{1},c_{1};\theta)$. Denoting with $\partial_{i}$ the partial
derivative in $(\theta)_{i}$, we thus have

$0=\partial_{i}G(a_{1},b_{1},c_{1};\theta_{1})=a_{1}\partial_{i}P(\theta_{1})-b_{1}\partial_{i}\Delta
P(\theta)+c_{1}\partial_{i}\eta(\theta_{1}).$ (S19)

Now we compute the derivative in $\theta$ of $F(a_{2},b_{2},c_{2};\theta)$,
where $a_{2},b_{2},c_{2}>0$ are three arbitrary scalars, and we evaluate it in
$\theta_{1}$. We have

$\partial_{i}F(a_{2},b_{2},c_{2};\theta_{1})=(a_{2}-c_{2}\partial_{P}\Sigma)\partial_{i}P(\theta_{1})-b_{2}\partial_{i}\Delta
P(\theta)-(c_{2}\partial_{\eta}\Sigma)\partial_{i}\eta(\theta_{1}).$ (S20)

Therefore, if we choose $a_{2},b_{2},c_{2}$ such that

$\begin{pmatrix}a_{1}\\\ b_{1}\\\
c_{1}\end{pmatrix}=\begin{pmatrix}1&0&-\partial_{P}\Sigma\\\ 0&1&0\\\
0&0&-\partial_{\eta}\Sigma\end{pmatrix}\begin{pmatrix}a_{2}\\\ b_{2}\\\
c_{2}\end{pmatrix},$ (S21)

thanks to Eq. (S19) we have that

$0=\partial_{i}F(a_{2},b_{2},c_{2};\theta_{1}),$ (S22)

meaning that the same parameters $\theta_{1}$ that nullifies the gradient of
$G$, nullifies also the gradient of $F$ at a different choice of the weights,
given by Eq. (S21). The invertibility of Eq. (S21) (i.e. a non-null
determinant of the matrix) is guaranteed by Eq. (S17). We also have to make
sure that if $a_{1},b_{1},c_{1}>0$, then also $a_{2},b_{2},c_{2}>0$. To do
this, we invert Eq. (S21), finding

$\begin{pmatrix}a_{2}\\\ b_{2}\\\
c_{2}\end{pmatrix}=\begin{pmatrix}1&0&-\partial_{P}\Sigma/(\partial_{\eta}\Sigma)\\\
0&1&0\\\ 0&0&-1/(\partial_{\eta}\Sigma)\end{pmatrix}\begin{pmatrix}a_{1}\\\
b_{1}\\\ c_{1}\end{pmatrix}.$ (S23)

It is now easy to see that also the weights $a_{2},b_{2},c_{2}$ are positive
using Eq. (S17).

To conclude the proof, we show that $\theta_{1}$ is a local maximum for
$F(a_{2},b_{2},c_{2};\theta)$ by showing that its Hessian is negative semi-
definite. Since, by hypothesis, $\theta_{1}$ is a local maximum for
$G(a_{1},b_{1},c_{1};\theta)$, we have that the Hessian matrix

$H^{(G)}_{ij}\equiv\partial_{ij}G(a_{1},b_{1},c_{1};\theta_{1})=a_{1}\partial_{ij}P-b_{1}\partial_{ij}\Delta
P+c_{1}\partial_{ij}\eta$ (S24)

is negative semi-definite. We now compute the Hessian $H^{(F)}$ of
$F(a_{2},b_{2},c_{2};\theta)$ in $\theta=\theta_{1}$:

$H^{(F)}_{ij}=a_{2}\partial_{ij}P-b_{2}\partial_{ij}\Delta
P-c_{2}\left[\partial_{P}\Sigma\,\partial_{ij}P+\partial_{\eta}\Sigma\,\partial_{ij}\eta+Q_{ij}\right],$
(S25)

where

$Q_{ij}=\begin{pmatrix}\partial_{i}P&\partial_{i}\eta\end{pmatrix}H^{(\Sigma)}\begin{pmatrix}\partial_{j}P\\\
\partial_{j}\eta\end{pmatrix},$ (S26)

and $H^{(\Sigma)}$ is the Hessian of $\Sigma(P,\eta)$ computed in
$P(\theta_{1})$ and $\eta(\theta_{1})$. Since we are interested in studying
the Hessian of $F(a_{2},b_{2},c_{2};\theta_{1})$ in the special point
$(a_{2},b_{2},c_{2})$ previously identified, we substitute Eq. (S23) into Eq.
(S25), yielding

$H^{(F)}_{ij}=H^{(G)}_{ij}+\frac{b_{1}}{\partial_{\eta}\Sigma}Q_{ij}.$ (S27)

We now prove that $H^{(F)}_{ij}$ is negative semi-definite since it is the sum
of negative semi-definite matrices. By hypothesis $H^{(G)}_{ij}$ is negative
semi-definite. Recalling Eq. (S17) and that $b_{1}>0$, we now need to show
that $Q_{ij}$ is positive semi-definite. Plugging Eq. (S18) into Eq. (S26)
yields

$Q_{ij}=\beta_{\text{C}}\frac{\eta_{\text{c}}}{\eta^{2}}\partial_{i}\eta\,\partial_{j}\eta\,R_{ij},$
(S28)

where

$\displaystyle R_{ij}$ $\displaystyle\equiv 2P+S_{ij}+S^{T}_{ij},$
$\displaystyle S_{ij}$
$\displaystyle=-\frac{\partial_{i}P}{\partial_{i}\eta}.$ (S29)

We now show that if $R_{ij}$ is positive semi-definite, then also $Q_{ij}$ is
positive semi-definite. By definition, $Q_{ij}$ is positive semi-definite if,
for any set of coefficient $a_{i}$, we have that
$\sum_{ij}a_{i}Q_{ij}a_{j}\geq 0$. Assuming $R_{ij}$ to be positive semi-
definite, and using that $\beta_{\text{C}},\eta_{\text{c}},\eta>0$, we have

$\sum_{ij}a_{i}Q_{ij}a_{j}=\beta_{\text{C}}\frac{\eta_{\text{c}}}{\eta^{2}}\sum_{ij}x_{i}R_{ij}x_{j}\geq
0,$ (S30)

where we define $x_{i}\equiv\partial_{i}\eta\,a_{i}$. We thus have to prove
the positivity of $R_{ij}$. We prove this showing that it is the sum of $3$
positive semi-definite matrices. Indeed, the first term in Eq. (S29), $2P$, is
proportional to a matrix with $1$ in all entries. Trivially, this matrix has
$1$ positive eigenvalue, and all other ones are null, so it is positive semi-
definite. At last, $S_{ij}$ and its transpose have the same positivity, so we
focus only on $S_{ij}$. $S_{ij}$ is a matrix with all equal columns. This
means that it has all null eigenvalues, except for a single one that we denote
with $\lambda$. Since the trace of a matrix is equal to the sum of the
eigenvalues, we have $\lambda=\mathrm{Tr}[S]=\sum_{i}S_{ii}$. Using the
optimality condition in Eq. (S19), we see that each entry of $S$ is positive,
i.e. $S_{ij}>0$. Therefore $\lambda>0$, thus $S$ is positive semi-definite,
concluding the proof that $H^{(F)}_{ij}$ is negative semi-definite.

To conclude, we notice that we can always renormalize $a_{2},b_{2},c_{2}$,
such that they sum to $1$, preserving the same exact optimization problem.

###  I.3 Identifying optimal cycles with reinforcement learning

In this subsection we show how the formulation of the optimization problem
given in Sec. I.1 in terms of an extended state allows us to use RL to
optimize a figure of merit containing power fluctuations, and we provide
details on the RL implementation. The RL method that we use is based on the
Soft Actor-Critic algorithm [99], introduced in the context of robotics and
video-games [97, 98, 119, 120], generalized to optimize multiple objectives.
RL has received great attention for its success at mastering tasks beyond
human-level such as playing games [121, 122, 123], for robotic applications
[97, 98], and for controlling plasmas in a tokamak [124]. RL has been recently
used for quantum control [105, 107, 110, 111, 112, 114, 115, 116],
outperforming previous state-of-the-art methods [108, 109], for fault-tolerant
quantum computation [106, 113], and in the field of quantum thermodynamics
[104, 100, 102, 101]. RL also allows to optimize “blackbox systems”, i.e. to
perform an optimization without requiring any explicit knowledge of the
dynamics of the system being optimized, thus potentially applicable to
experimental devices [102, 101]. In light of this, we expect RL to be a
powerful tool also in this context. Furthermore, since RL is based on the
Markovian Decision Process framework (see below and Ref. [73] for details), it
is a natural choice for systems described by Markovian dynamics. In principle,
also other optimization methods, such as Pontryagin’s Minimum Principle, could
be used [72]. However, this method only provides, in general, necessary
conditions for an optimal control. Furthermore, it could in practice get stuck
in local maxima, and it requires some hand-tuning in the case of non-analytic
controls. As we see from our results, optimal cycles are often given by of
piece-wise continuous function, which are not analytic. This would complicate
the use of Pontryagin’s Minimum Principle.

The method here described is an extension of the methods described in Refs.
[100, 101]. We therefore refer to the Ref. [100] for an explanation of the RL
method, and we adopt its notation to explain all the differences and
generalizations that we put forward in this paper.

First, we define the optimization problem detailing the choice of the state
space, the action space, and of the reward. Let us discretize time in time
steps $t_{i}=i\Delta t$. In Ref. [100] just the power is optimized, so the
state of the environment is described by the density matrix $\rho$. Indeed, RL
is based on the Markov Decision Process assumption, meaning that the reward at
a given time-step must be a function of the last state and action. In Ref.
[100], the reward is the power averaged over the last time-step, which can be
computed from the density matrix and from the value of the control at the last
time-step. However, power fluctuations cannot be computed solely from the
knowledge of $\rho$. Here, to optimize also power fluctuations, we employ the
extended state as state of the environment, i.e. we choose as state space
$\mathcal{S}=\\{(\rho,s,\vec{u})|\rho\in\mathcal{D},s\in\mathcal{E},\vec{u}\in\mathcal{U}\\}$,
where $\mathcal{D}=\\{\rho|\rho\geq 0,\Tr[\rho]=1\\}$ is the space of density
matrices, $\mathcal{E}=\\{s|s=s^{\dagger},\Tr[s]=0\\}$ is the space of
traceless Hermitian operators, and $\mathcal{U}$ is the continuous set of
allowed control parameters. The state at time-step $i$ is then given by
$s_{i}=(\rho_{t_{i}},s_{t_{i}},\vec{u}(t_{i-1}))$. The action space is given
by
$\mathcal{A}=\\{(\vec{u},\beta)|\vec{u}\in\mathcal{U},\beta\in[\beta_{\text{H}},\beta_{\text{C}}]\\}$,
such that the action at time-step $t_{i}$ is $a_{i}=(\vec{u}_{i},\beta_{i})$.
The controls $\vec{u}(t)$ and $\beta(t)$ are then kept constant at
$\vec{u}_{i}$ $\beta_{i}$ for the current time-step. The reward is given by

$r_{i}=\frac{1}{\Delta
t}\int_{t_{i-1}}^{t_{i}}\left[a\frac{{P}(t)}{P_{\text{max}}}-b\frac{\Delta{P}(t)}{\Delta
P(P_{\text{max}})}-c\frac{\Sigma(t)}{\Sigma(P_{\text{max}})}\right]dt,$ (S31)

which corresponds to the figure of merit averaged only over the last time-
step. Together with Eqs. (3) and (7) of the main text, we see that this choice
satisfies the Markov Decision Process assumption.

As in Ref. [100], we formulate the optimization problem as a discounted,
continuing RL task, where the aim is to learn a policy $\pi(a|s)$ that, at
each time step $t_{i}$, maximizes the return, i.e. the long-term average of
the rewards:

$r_{i+1}+\gamma
r_{i+2}+\gamma^{2}r_{i+3}+\dots=\sum_{k=0}^{\infty}\gamma^{k}r_{i+1+k},$ (S32)

where $\gamma\in[0,1)$ is the discount factor which determines how much we are
interested in future rewards, as opposed to immediate rewards. By choosing
$\gamma$ close enough to $1$, we are optimizing the figure of merit averaged
over a long time-scale, such that the method should automatically discover to
perform cycles (see Refs. [100, 101] for additional details).

In this work, we employ the soft actor-critic method [97, 98] as implemented
in Ref. [100] with the following changes.

  * •

Policy parameterization. Here we do not have a discrete action, but we have
multiple continuous actions given by $\vec{\xi}\equiv(\vec{u},\beta)$. We
therefore need to parameterize the policy $\pi(\vec{\xi}|s)$ as a multi-
variate distribution function. We use the same approach as in Ref. [100], but
we replace the normal distribution with a multivariate normal distribution. To
this end, we employ a Neural Network (NN) that takes the state $s$ as input,
and outputs a vector $\vec{\mu}$ and a matrix $M$. We use a multilayer
perceptron NN with $2$ hidden layer. We then produce samples from the
multivariate normal distribution with mean $\vec{\mu}$ and covariance matrix
$M^{T}M+\lambda\mathbb{1}$, where $\lambda>0$ is a small real number added for
numerical stability. Then, as in Ref. [100], we determine a sample of
$\pi(\vec{\xi}|s)$ applying a “squash function”, in this case a hyperbolic
tangent, to each variable of the multivariate distribution. This ensures that
each action lies in the specified continuous interval $\mathcal{U}$.

  * •

Automatic “temperature” tuning. In the soft-actor critic method, the
exploration-exploitation balance is governed by the hyper-parameter
$\varepsilon$, known in the RL literature as the “temperature” parameter. In
Ref. [100], $\varepsilon$ was scheduled during training. However, its value
depends on the magnitude of the reward, making it quite model-dependent. Here,
instead, we follow Ref. [97] to automatically tune $\varepsilon$ during
training. The idea is to change $\varepsilon$ during training as to set the
average entropy of the policy to a target value $\bar{H}$. Large values of
$\bar{H}$ produce an exploratory policy, while smaller values of $\bar{H}$
produce a more deterministic policy. As detailed in Ref. [97] and implemented
in Ref. [101], the tuning of $\varepsilon$ is done by minimizing the following
loss function

$L_{H}(\varepsilon)\equiv\varepsilon\underset{\begin{subarray}{c}s\sim\mathcal{B}\end{subarray}}{\text{E}}\left[H(\pi(\cdot|s))-\bar{H}\right],$
(S33)

where $H(P)\equiv\underset{\begin{subarray}{c}x\sim
P\end{subarray}}{\text{E}}[-\ln P(x)]$ is the entropy of the distribution
$P(x)$. For numerical stability, if $\varepsilon$ becomes negative during
training, we reset it to a small positive number. To favour exploration in the
early phases of the training, while still obtaining a final policy that is
nearly deterministic, we exponentially schedule the target entropy $\bar{H}$
during training as follows

$\bar{H}(n_{\text{steps}})=\bar{H}_{\text{end}}+(\bar{H}_{\text{start}}-\bar{H}_{\text{end}})\exp(-n_{\text{steps}}/\bar{H}_{\text{decay}}),$
(S34)

where $\bar{H}_{\text{start}}$, $\bar{H}_{\text{end}}$ and
$\bar{H}_{\text{decay}}$ are training hyperparameters.

  * •

Training steps. The training steps are performed as in Ref. [100], i.e.
repeatedly using the ADAM algorithm to minimize the loss functions $L_{\pi}$
and $L_{Q}$ computed over a batch of experience drawn from the replay buffer
$\mathcal{B}$. Here, in addition, every time an optimization step of $L_{\pi}$
and $L_{Q}$ is performed, we also update $\varepsilon$ performing one
optimization step of $L_{H}$ computed over the same batch of experience.

  * •

Optimal cycle evaluation. Once the training is complete, all cycles and values
of the power, power fluctuations, entropy production and efficiency reported
in the main text are computed evaluating the deterministic policy. More
specifically, we determine the optimal cycle using the deterministic policy,
i.e. choosing actions according to the mean of the multivariate normal
distribution, instead of sampling from it. This turns the stochastic policy
into a deterministic one that typically performs better.

##  II Fast-Driving for generic systems

In this section we derive analytic expressions for the main quantities (S8) in
the regime of fast-driving, i.e. when the period $\tau$ of the cycle is small
compared to the relaxation timescales of the system. As this regime is
rigorously justified for stochastic dynamics, and for easiness of
presentation, we will assume diagonal operators that commute at all times (cf.
Ref. [27]).

First, for what concerns power and entropy production, notice, from (S8) that
both can be computed from the integral of

$\displaystyle\Tr[\dot{\rho}_{t}H_{\vec{u}(t)}]=\Tr[\mathcal{L}_{\vec{u}(t),\beta(t)}[\rho_{t}]H_{\vec{u}(t)}]\;.$
(S35)

In the fast driving approximation, the cycle is much quicker than the
relaxation time of the system. The consequence is that
$\rho_{t}\simeq\rho^{(0)}$ tends to a steady-state solution, which can be
computed as the leading order of a perturbative expansion [27]. It follows
that power and entropy production can be computed as

$\displaystyle\expectationvalue{P}$
$\displaystyle=\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[\mathcal{L}_{\vec{u}(t),\beta(t)}[\rho^{(0)}]H_{\vec{u}(t)}],$
(S36) $\displaystyle\expectationvalue{\Sigma}$
$\displaystyle=-\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[\mathcal{L}_{\vec{u}(t),\beta(t)}[\rho^{(0)}]H_{\vec{u}(t)}]\beta(t).$
(S37)

The generic expression for $\rho^{(0)}$ is given in Ref. [27]. At the same
time, for the purposes of the present work, it is enough to consider the
simple case in which the dynamics is described by a scalar master equation
$\mathcal{L}[...]=\gamma\pi_{t}\Tr[...]-\gamma\mathbb{1}$, where

$\displaystyle\pi_{t}=\frac{e^{-\beta(t)H_{\vec{u}(t)}}}{\Tr[e^{-\beta(t)H_{\vec{u}(t)}}]}$
(S38)

is the instantaneous steady state and $\gamma^{-1}$ the relaxation timescale.
In such case it is clear that, as $\rho^{(0)}$ is the steady state, it
satisfies

$\displaystyle\int_{0}^{\tau}dt\gamma\left(\pi_{t}-\rho^{(0)}\right)=0\;,$
(S39)

leading to

$\displaystyle\rho^{(0)}=\frac{\int_{0}^{\tau}dt\
\pi_{t}}{\tau}=\expectationvalue{\pi}\;,$ (S40)

and therefore

$\displaystyle\expectationvalue{P}$
$\displaystyle=\frac{1}{\tau}\int_{0}^{\tau}dt\gamma\Tr[(\pi_{t}-\expectationvalue{\pi})H_{\vec{u}(t)}],$
(S41) $\displaystyle\expectationvalue{\Sigma}$
$\displaystyle=-\frac{1}{\tau}\int_{0}^{\tau}dt\gamma\Tr[(\pi_{t}-\expectationvalue{\pi})H_{\vec{u}(t)}]\beta(t).$
(S42)

The most general expression for the heat current in the fast-driving regime is
given in Ref. [27]. In the following, we derive new analytical expressions
that allow to efficiently compute also the power fluctuations in the same
regime.

###  II.1 Fluctuations in the Fast driving regime for arbitrary stochastic
(classical) systems

As mentioned above, we will simplify the discussion by considering semi-
classical systems represented by commuting operators at all time.

####  II.1.1 Simple Lindbladian

First, we start again by considering the simple scalar master equation
$\mathcal{L}[...]=\gamma\pi_{t}\Tr[...]-\gamma\mathbb{1}$, with
$\pi_{t}=e^{-\beta(t)H_{\vec{u}(t)}}/\Tr[e^{-\beta(t)H_{\vec{u}(t)}}]$. In
such case, the equations for $\rho_{t}$ (S1) and $s_{t}$ (S7) become (from now
on we simplify the notation as $H_{\vec{u}(t)}\equiv H_{t}$)

$\displaystyle\dot{\rho}_{t}=$ $\displaystyle\gamma(\pi_{t}-\rho_{t})\;,$
(S43) $\displaystyle\dot{s}_{t}=$ $\displaystyle-\gamma
s_{t}+2\rho_{t}\dot{H}_{t}-2\Tr[\rho_{t}\dot{H}_{t}]\rho_{t}\equiv-\gamma
s_{t}+2\rho_{t}{\tilde{\dot{H}}_{\rho}}_{t}\;,$ (S44)

where we defined

$\displaystyle\tilde{A}_{\rho}=A-\rho\Tr[\rho A]\;.$ (S45)

The asymptotic solutions for $\pi_{t}$ and $s_{t}$, in the case of periodic
driving with period $\tau$, are

$\displaystyle\rho_{t}=$
$\displaystyle\frac{\int_{t-\tau}^{t}dt^{\prime}\;\gamma\pi_{t^{\prime}}e^{\gamma(t^{\prime}-t)}}{1-e^{-\gamma\tau}}\;,$
(S46) $\displaystyle s_{t}=$
$\displaystyle\frac{\int_{t-\tau}^{t}dt^{\prime}\;2\rho_{t^{\prime}}{\tilde{\dot{H}}_{\rho}}_{t^{\prime}}e^{\gamma(t^{\prime}-t)}}{1-e^{-\gamma\tau}}\;.$
(S47)

Notice that until now, no approximation has been used. To compute all relevant
quantities, we will use the fast-driving condition $\gamma\tau\ll 1$ and
expand in powers of $\gamma\tau$. That is,

$\displaystyle\rho_{t}=\rho^{(0)}_{t}+\rho^{(1)}_{t}+\dots\quad\rho^{(i)}_{t}\sim\mathcal{O}(\gamma^{i}\tau^{i})\;,$
(S48) $\displaystyle s_{t}=s^{(0)}_{t}+s^{(1)}_{t}+\dots\quad
s^{(i)}_{t}\sim\mathcal{O}(\gamma^{i}\tau^{i})\;.$ (S49)

It is easy to verify, from (S46),

$\displaystyle\rho^{(0)}_{t}$
$\displaystyle\equiv\rho^{(0)}=\expectationvalue{\pi}\;,$ (S50)
$\displaystyle\dot{\rho}^{(0)}_{t}$ $\displaystyle=0\;,$ (S51)
$\displaystyle\dot{\rho}^{(1)}_{t}$
$\displaystyle=\gamma(\pi_{t}-\rho^{(0)})\;.$ (S52)

The expansion of $s_{t}$ is slightly more complex, but for our purposes we
only need $s^{(0)}$, which can be found integrating by parts (S47),

$\displaystyle
s_{t}=\frac{\int_{t-\tau}^{t}ds\;2\rho_{s}{\tilde{\dot{H}}_{\rho}}_{s}e^{\gamma(s-t)}}{1-e^{-\gamma\tau}}=2\rho_{t}{{\tilde{H}}_{\rho_{t}}}-\frac{2}{1-e^{-\gamma\tau}}\int_{t-\tau}^{t}dt^{\prime}\left(\dot{\rho}_{t^{\prime}}{\tilde{H}_{\rho_{t^{\prime}}}}e^{\gamma(t^{\prime}-t)}+\gamma\rho_{t^{\prime}}{\tilde{H}_{\rho_{t^{\prime}}}}e^{\gamma(t^{\prime}-t)}-\rho_{t^{\prime}}\Tr[H_{t^{\prime}}\dot{\rho}_{t^{\prime}}]e^{\gamma(t^{\prime}-t)}\right)\;.$
(S53)

In this expression we can substitute all the quantities at the leading order,
obtaining, using the fluctuations notation $\delta A:=A-\langle A\rangle$,

$\displaystyle
s^{(0)}_{t}=2\rho^{(0)}\tilde{H}_{\rho^{(0)}}-\frac{2}{\gamma\tau}\int_{t-\tau}^{t}dt^{\prime}\left(\gamma\delta\pi_{t^{\prime}}\tilde{H}_{\rho^{(0)}}+\gamma\rho^{(0)}\tilde{H}_{\rho^{(0)}}-\gamma\rho^{(0)}\Tr[H_{t^{\prime}}\delta\pi_{t^{\prime}}]\right)$
(S54)

Notice that although $\rho^{(0)}$ is constant, the quantity
$\tilde{H}_{\rho^{(0)}}$ (cf. Eq. (S45)) is time-dependent as $H_{t}$ is time
dependent. The above expression can be further simplified by dropping all the
time-dependence in the notation and expanding the time-average at leading
order

$\displaystyle
s^{(0)}=2\left(\rho^{(0)}\tilde{H}_{\rho^{(0)}}-\langle\rho^{(0)}\tilde{H}_{\rho^{(0)}}\rangle\right)+2\bigg{\langle}\rho^{(0)}\Tr[\delta
H\delta\pi]+\delta\pi\Tr[\delta H\rho^{(0)}]-\delta\pi\delta
H\bigg{\rangle}\;.$ (S55)

##### Fluctuations.

The fluctuations become, integrating by parts (S8),

$\displaystyle\Delta
P=\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[\dot{H}_{t}s_{t}]=-\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[H_{t}\dot{s}_{t}]=\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[H_{t}(\gamma
s_{t}-2\tilde{\dot{H}}_{\rho_{t}})]\;.$ (S56)

Therefore we want to compute

$\displaystyle I_{2}=$
$\displaystyle-\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[2H_{t}\tilde{\dot{H}}_{\rho_{t}}]\;,$
(S57) $\displaystyle I_{1}=$
$\displaystyle\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[\gamma H_{t}s_{t}]\;,$ (S58)

at leading order. Let’s compute the two terms

$I_{2}=-\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[2H_{t}\tilde{\dot{H}}_{\rho_{t}}]=-\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[2H_{t}(\dot{H}_{t}\rho_{t}-\rho_{t}\Tr[\dot{H}_{t}\rho_{t}])]=\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[H_{t}^{2}\dot{\rho}_{t}-2\Tr[\rho_{t}H_{t}]\Tr[\dot{\rho}_{t}H_{t}]]=\\\
=\gamma\left(\langle\Tr[H^{2}\delta\pi]\rangle-2\langle\Tr[H\delta\pi]\Tr[H\rho^{(0)}]\rangle\right)\;,$
(S59)

where for the last equality we used the equation of motions
$\dot{\rho}_{t}=\gamma(\pi_{t}-\rho_{t})$ a the leading order, where
$\rho_{t}\sim\rho^{(0)}=\langle\pi\rangle$. We remind that we use the notation
$\delta A:=A-\langle A\rangle$. Notice that for $I_{1}$ we only need $s$ at
order $\mathcal{O}(\gamma^{0}\tau^{0})$. $I_{1}$ can therefore be computed as

$\displaystyle I_{1}=\gamma\langle\Tr[Hs^{(0)}]\rangle\;,$ (S60)

and the total fluctuations are, after substituting $s^{(0)}$ (S55) and some
tedious algebra,

$\frac{\Delta P}{2\gamma}=\frac{I_{2}+I_{1}}{2\gamma}=\\\
=\expectationvalue{\Tr[\rho^{(0)}H^{2}]-\Tr[\rho^{(0)}H]^{2}}-\left(\Tr[\rho^{(0)}\expectationvalue{H}^{2}]-\Tr[\rho^{(0)}\expectationvalue{H}]^{2}\right)+\frac{1}{2}\langle\Tr[\delta
H^{2}\delta\pi]\rangle-\langle\Tr[\delta H\delta\pi]\Tr[\delta
H\rho^{(0)}]\rangle\;.$ (S61)

Under the same assumptions, the power (S41) can be expressed as

$P=\gamma\langle\Tr[(\pi-\langle\pi\rangle)(H-\langle
H\rangle)]\rangle\equiv\gamma\Tr[\text{Cov}[\pi,H]],$ (S62)

and similarly the entropy production (S42)

$\Sigma=\gamma\langle\Tr[(\pi-\langle\pi\rangle)(\beta H-\langle\beta
H\rangle)]\rangle\equiv\gamma\Tr[\text{Cov}[\pi,\beta H]],$ (S63)

where the covariance is with respect to time.

####  II.1.2 General Stochastic Classical Case

Moving to the case of a general Lindbladian, the expression for power and
entropy production can be found as from Reference [27]. For what concerns the
fluctuations, these can be again expressed as (see (S56))

$\displaystyle\Delta
P=-\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[H_{t}\dot{s}_{t}]=-\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[H_{t}(\mathcal{L}_{t}[s_{t}]+2\tilde{\dot{H}}_{\rho_{t}})]=I_{1}+I_{2}\;.$
(S64)

In the next expressions we drop all the time-dependencies to simplify the
notation. The computation of $I_{1}$ and $I_{2}$ is in this case

$\displaystyle I_{2}=$
$\displaystyle-\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[2H\tilde{\dot{H}}_{\rho}]=\langle\Tr[H^{2}\mathcal{L}[\rho^{(0)}]]\rangle-2\langle\Tr[H\mathcal{L}[\rho^{(0)}]]\Tr[H\rho^{(0)}]\rangle\;,$
(S65) $\displaystyle I_{1}=$
$\displaystyle-\frac{1}{\tau}\int_{0}^{\tau}dt\Tr[H\mathcal{L}[s_{0}]]\;,$
(S66)

where the first expression can be obtained integrating by parts, and
everything is expanded at the leading order in $|\mathcal{L}|\tau$,
$\rho^{(0)}$ is the steady state of the fast driving (see [27] for an
analytical expression) and $s_{0}$ is the leading order term of $s$, which can
be computed from the analytical solution

$\displaystyle
s_{t}=\int_{-\infty}^{t}dt^{\prime}\;U_{t,t^{\prime}}[2\rho_{t^{\prime}}\tilde{\dot{H}}_{\rho_{t^{\prime}}}]\;,$
(S67)

$U_{t,t^{\prime}}$ being the propagator, i.e. time-ordered exponential of
$\mathcal{L}$ between $t^{\prime}$ and $t\geq t^{\prime}$. Expanding the
leading term one obtains

$\displaystyle
s_{0}=2\rho^{(0)}\tilde{H}_{\rho^{(0)}}+2\langle\mathcal{L}\rangle^{-1}\left(\langle\mathcal{L}[\rho^{(0)}]\tilde{H}_{\rho^{(0)}}\rangle-\rho^{(0)}\Tr[\langle
H\mathcal{L}[\rho^{(0)}]\rangle]-\langle\mathcal{L}[\rho^{(0)}\tilde{H}_{\rho^{(0)}}]\rangle\right)\;.$
(S68)

##  III Quantum Dot Based Heat Engine

###  III.1 Optimality of the Otto cycle in the fast driving for power and
entropy production trade-offs

Here we prove that, among all possible cycles described by $\beta(t)$ and
$u(t)$, Otto cycles in the fast driving regime maximize an arbitrary trade-off
between power and entropy production. The proof stems from a simple
generalization of Appendix A of Ref. [23], where it is shown that, among all
possible Otto cycles, Otto cycles in the fast driving regime maximize the
power.

To this end, we consider the figure of merit $\expectationvalue{F}$ with
$b=0$:

$\expectationvalue{F}=a\frac{\expectationvalue{P}}{P_{\text{max}}}-c\frac{\expectationvalue{{\Sigma}}}{\Sigma(P_{\text{max}})}.$
(S69)

By neglecting fluctuations, we no not need the extended state, so we can write
$\expectationvalue{F}$ as the time average of a function of the state $p_{t}$
(defined in the main text) and of the controls $\beta(t)$ and $u(t)$, i.e. as

$\expectationvalue{F}=\frac{1}{\tau}\int_{0}^{\tau}f(p_{t},u(t),\beta(t))dt,$
(S70)

where $f(p,u,\beta)$ is a suitable function. This assumption is sufficient to
use the argument of [23], with $\expectationvalue{F}$ replacing the average
power, to prove that optimal cycles are Otto cycles in the fast-driving
regime. Here we outline the general idea, referring to [23] for additional
details. As argued here and in [23], given a cycle with period $\tau$, also
the state $p_{t}$ will become periodic with the same period $\tau$. We can
thus represent an arbitrary cycle cycle as a closed curve in the $u-p$ plane
and in the $\beta-p$ plane. It can be shown that, given any fixed cycle $u(t)$
and $\beta(t)$, we can always define two sub-cycles such that one of the two
has a larger or equal $\expectationvalue{F}$ than the original one. The two
sub-cycles are defined by introducing quenches in the control, i.e. “cutting
vertically” the cycle represented in the $u-p$ and $\beta-p$ plane. By
reiterating this process over and over, we end up with an Otto cycle in the
fast driving regime, thus concluding the proof.

###  III.2 Fast driving regime

Using Eqs. (S61,S62,S63), we can easily compute the power, entropy production
and fluctuations of a generic Otto cycle in the fast driving regime. We fix an
Otto cycle as in Fig. S1, where $\epsilon_{\text{H}}$ and
$\epsilon_{\text{C}}$ are the values of $u(t)E_{0}$ while in contact
respectively with the hot (H) and cold (C) bath, and $\theta_{\text{H}}$ and
$\theta_{\text{C}}$ represent the time fraction of the cycle spent in contact
with the respective bath ($\theta_{\text{C}}+\theta_{\text{H}}=1$).

Figure S1: Schematic representation of the Otto cycle. The value of the gap
$\epsilon(t)\equiv u(t)E_{0}$ is plotted as a function of $t/\tau$. The red
segment corresponds to setting $\beta(t)=\beta_{\text{H}}$, and the blue
segment to $\beta(t)=\beta_{\text{C}}$.

This yields

$\displaystyle\expectationvalue*{P}$
$\displaystyle=\gamma\theta_{\text{H}}\theta_{\text{C}}\,(\epsilon_{\text{H}}-\epsilon_{\text{C}})(f_{\text{H}}-f_{\text{C}}),$
(S71) $\displaystyle\expectationvalue*{\Delta P}$
$\displaystyle=2\gamma\,\theta_{\text{H}}\theta_{\text{C}}\,(\epsilon_{\text{H}}-\epsilon_{\text{C}})^{2}\left[\bar{f}(1-\bar{f})+(\frac{f_{\text{C}}+f_{\text{H}}}{2}-\bar{f})(1-2\bar{f})\right],$
$\displaystyle\expectationvalue*{\Sigma}$
$\displaystyle=-\gamma\theta_{\text{H}}\theta_{\text{C}}(\beta_{\text{H}}\epsilon_{\text{H}}-\beta_{\text{C}}\epsilon_{\text{C}})(f_{\text{H}}-f_{\text{C}}),$

where we define $f_{\alpha}=f(\beta_{\alpha}\epsilon_{\alpha})$, for
$\alpha=\text{H},\text{C}$ and the average excited state occupation
$\bar{f}=\theta_{\text{H}}f_{\text{H}}+\theta_{\text{C}}f_{\text{C}}$. Without
loss of generality, in the heat engine regime we can assume
$\epsilon_{\alpha}\geq 0$. This condition, together with the heat engine
condition $P\geq 0$, can be expressed as

$0\leq x_{\text{H}}\leq x_{\text{C}}\leq x_{\text{H}}(1+dT),$ (S72)

where we introduce the dimensionless quantities
$x_{\alpha}=\beta_{\alpha}\epsilon_{\alpha}$, and where we define the
dimensionless temperature difference

$dT\equiv\beta_{\text{C}}(\beta_{\text{H}}^{-1}-\beta_{\text{C}}^{-1}).$ (S73)

For simplicity, in this appendix we define a figure of merit $G$ with a
different normalization, i.e.

$\expectationvalue*{G}=a\left(\frac{\expectationvalue*{P}}{\gamma
T}\right)-b\left(\frac{\expectationvalue*{\Delta P}}{\gamma
T^{2}}\right)-c\left(\frac{\expectationvalue*{\Sigma}}{\gamma}\right),$ (S74)

where $T=\beta_{\text{C}}^{-1}$ and $k_{B}=1$.

We now assume that the temperature difference is small, i.e. that $dT\ll 1$.
We therefore study $\expectationvalue*{G}$ to leading order in $dT$. We carry
out the optimization with respect to the time fractions $\theta_{\alpha}$ and
with respect to the dimensionless energy gaps $x_{\alpha}$. We expand
$x_{\text{C}}$ as

$x_{\text{C}}=x_{\text{H}}(1+\delta x_{\text{C}}).$ (S75)

Imposing the heat engine condition in Eq. (S72), we find that $\delta
x_{\text{C}}$ must satisfy

$0\leq\delta x_{\text{C}}\leq dT;$ (S76)

therefore, $\delta x_{\text{C}}$ is a first order quantity in $dT$. Expanding
$\expectationvalue*{P}$, $\expectationvalue*{\Delta P}$ and
$\expectationvalue*{\Sigma}$ to leading order in $dT$, we find

$\displaystyle\expectationvalue*{P}/(\gamma T)$
$\displaystyle=\theta_{\text{C}}\theta_{\text{H}}g(x_{\text{H}})\delta
x_{\text{C}}(dT-\delta x_{\text{C}}),$ (S77)
$\displaystyle\expectationvalue*{\Delta P}/(\gamma T^{2})$
$\displaystyle=2\theta_{\text{C}}\theta_{\text{H}}g(x_{\text{H}})(dT-\delta
x_{\text{C}})^{2},$ $\displaystyle\expectationvalue*{\Sigma}/\gamma$
$\displaystyle=\theta_{\text{C}}\theta_{\text{H}}g(x_{\text{H}})\delta
x_{\text{C}}^{2},$

where we define

$g(x)\equiv\frac{x^{2}}{2(1+\cosh{x})}.$ (S78)

Plugging this expansion into Eq. (S74) yields

$\expectationvalue*{G}=\theta_{\text{H}}\theta_{\text{C}}g(x_{\text{H}})\left[\delta
x_{\text{C}}\left(adT-(a+c)\delta x_{\text{C}}\right)-2b\left(dT-\delta
x_{\text{C}}\right)^{2}\right].$ (S79)

First, we notice that $\expectationvalue*{P}$, $\expectationvalue*{\Delta P}$
and $\expectationvalue*{\Sigma}$ in Eq. (S77), valid in the small temperature
difference regime, exactly saturate the steady-state thermodynamic uncertainty
relation even before performing any optimization, i.e. $\xi=1$, where,
combining Eqs. (1) of the main text and (S10)

$\xi=2\frac{\expectationvalue*{P}^{2}}{\expectationvalue*{\Sigma}\expectationvalue*{\Delta
P}}.$ (S80)

We now maximize $\expectationvalue*{G}$ as written in Eq. (S79). The
optimization over $\theta_{\alpha}$ and $x_{\text{H}}$ depends on the term in
the square parenthesis. We therefore distinguish two cases:

####  III.2.1 Case 1

If the term in square parenthesis is non-null, then we trivially have that
$\theta_{\alpha}=1/2$ and $x_{\text{H}}=x_{\text{max}}$, where
$x_{\text{max}}$ is the value that maximizes $g(x)$. It is implicitly defined
by $x_{\text{max}}>0$ and by solving the transcendental equation

$x_{\text{max}}\tanh(x_{\text{max}}/2)=2.$ (S81)

We can then determine $\delta x_{\text{C}}$ by taking the derivative of $G$
and setting it to zero. This yields

$\delta x_{\text{C}}^{*1}=\frac{a+4b}{2a+4b+2c}dT.$ (S82)

It is easy to see that this solution always satisfies the heat engine
condition in Eq. (S76). In this point we have that

$\expectationvalue*{G}^{*1}=g_{\text{max}}\frac{a^{2}-8bc}{16(a+2b+c)}dT^{2},$
(S83)

where we define $g_{\text{max}}\equiv g(x_{\text{max}})\approx 2.40$, and

$\displaystyle\expectationvalue*{P}^{*1}/(\gamma T)$
$\displaystyle=g_{\text{max}}\frac{(a+4b)(a+2c)}{16(a+2b+c)^{2}}dT^{2},$
$\displaystyle~{}\quad~{}\expectationvalue*{\Delta P}^{*1}/(\gamma T^{2})$
$\displaystyle=g_{\text{max}}\frac{(a+2c)^{2}}{8(a+2b+c)^{2}}dT^{2},$ (S84)
$\displaystyle\expectationvalue*{\Sigma}^{*1}/\gamma$
$\displaystyle=g_{\text{max}}\frac{(a+4b)^{2}}{16(a+2b+c)^{2}}dT^{2},$
$\displaystyle\frac{\eta^{*1}}{\eta_{\text{c}}}$
$\displaystyle=\frac{a+2c}{2a+4b+2c}.$

We thus see that if $a^{2}-8bc>0$, then the figure of merit is positive. In
the opposite case, the figure of merit is negative, so doing nothing becomes
more convenient (since doing nothing, or doing Carnot cycles, gives
$\expectationvalue*{G}=0$). Notice that if $a^{2}=8bc$, then any value of
$\theta_{\alpha}$ and $x_{\text{H}}$ gives the same figure of merit, so all
these solutions must be included in the Pareto front.

Setting $a=1$ and $b=c=0$, we find the maximum power solution:

$\displaystyle\frac{P_{\text{max}}}{\gamma T}$
$\displaystyle=\frac{g_{\text{max}}}{16}dT^{2},$ $\displaystyle\frac{\Delta
P(P_{\text{max}})}{\gamma T^{2}}$
$\displaystyle=\frac{g_{\text{max}}}{8}dT^{2},$
$\displaystyle\frac{\Sigma(P_{\text{max}})}{\gamma}$
$\displaystyle=\frac{g_{\text{max}}}{16}dT^{2},$
$\displaystyle\frac{\eta}{\eta_{\text{c}}}$ $\displaystyle=\frac{1}{2}.$ (S85)

These relations imply

$\Delta P(P_{\text{max}})=2TP_{\text{max}}.$ (S86)

Setting $c=0$ and $b=1-a$ in Eq. (S84), we can find the Pareto front between
power and power fluctuations, i.e. the outer border of Fig. 3 of the main
text. After some algebra, it can be shown that

$\frac{\expectationvalue*{P}}{P_{\text{max}}}=2\sqrt{\frac{\expectationvalue*{\Delta
P}}{\Delta P(P_{\text{max}})}}-\frac{\expectationvalue*{\Delta P}}{\Delta
P(P_{\text{max}})}.$ (S87)

This corresponds to the black border shown in Fig. 3 of the main text.

####  III.2.2 Case 2

We now consider the case where the term in the square parenthesis of Eq. (S79)
is zero. This happens in two cases, i.e. when

$\displaystyle\delta x_{\text{C}}^{*2}$
$\displaystyle=\frac{a+4b-\sqrt{a^{2}-8bc}}{2a+4b+2c}dT,$ (S88)
$\displaystyle\delta x_{\text{C}}^{*3}$
$\displaystyle=\frac{a+4b+\sqrt{a^{2}-8bc}}{2a+4b+2c}dT.$

These solutions hold when $a^{2}-8bc\geq 0$. In the opposite case, there is no
solution. Both solutions can be shown to always satisfy the heat engine
condition in Eq. (S76), and in this point we have

$\expectationvalue*{G}^{*2,3}=0.$ (S89)

However, this zero value of the figure of merit occurs delivering a positive
power that is compensated by the negative sign in front of the finite
fluctuations and of the finite entropy production. It can be shown that such
solutions either lie on the boundary Eq. (S87), or inside. We further notice
that, looking at the figure of merit, both these solutions are suboptimal with
respect to case $1$, except for being equivalent in the special case when
$a^{2}=8bc$.

####  III.2.3 Fast driving Pareto front

As we have seen, Otto cycles in the fast driving, expanded to leading order in
the temperature difference, exactly satisfy the steady-state thermodynamic
uncertainty relation $\xi=1$. Furthermore, the trade-off between power and
power fluctuations is given by Eq. (S87). Therefore, the entire Pareto front,
plotted in Fig. 3a of the main text, is simply obtained by imposing the
thermodynamic uncertainty relation, together with the boundary of Eq. (S87).

However, it could in principle be possible that points inside the boundary of
Eq. (S87) may not be reached, i.e. there may not be any Otto cycle in the fast
driving regime reaching some of these points. Using the results of the
previous subsections, it can be verified that this is not the case, i.e. that
there is an Otto cycle in the fast driving regime corresponding to all points
shown in Fig. 3a of the main text.

####  III.2.4 Positive figure of merit boundary for Otto cycles in the fast-
driving regime

In this section we derive the expression for the black boundary shown in Fig.
2b of the main text. As we have seen in the previous two subsections, the
condition

$a^{2}=8bc,$ (S90)

with $b=1-a-c$, determines the boundary between positive and zero value of the
optimized figure of merit in the fast-driving regime. However, this condition
was found with the normalization of $\expectationvalue{P}$,
$\expectationvalue{\Delta P}$ and $\expectationvalue{\Sigma}$ defined as in
Eq. (S74), while the boundary shown in Fig. 2 of the main text is relative to
the normalization chosen in Eq. (8) of the main text. To account for this, we
need to “change the coordinate system” from $(a_{0},b_{0},c_{0})$ to
$(a_{1},b_{1},c_{1})$. Let us define

$\expectationvalue*{G_{0}}\equiv
a_{0}\expectationvalue*{P}-b_{0}\expectationvalue*{\Delta
P}-c_{0}\expectationvalue*{\Sigma},$ (S91)

and let us consider a given cycle that maximize $G_{0}$ at given
$(a_{0},b_{0},c_{0})$. We want to see how the coefficients
$(a_{1},b_{1},c_{1})$ must be chosen in order to find the same cycle when
maximizing

$\expectationvalue*{G_{1}}=a_{1}\left(\frac{\expectationvalue*{P}}{\lambda_{a}}\right)-b_{1}\left(\frac{\expectationvalue*{\Delta
P}}{\lambda_{b}}\right)-c_{1}\left(\frac{\expectationvalue*{\Sigma}}{\lambda_{c}}\right),$
(S92)

where $\lambda_{a}$, $\lambda_{b}$ and $\lambda_{c}$ are given coefficients.

It can be shown that the following identity holds:

$\expectationvalue*{G_{0}}=N\expectationvalue*{G1},$ (S93)

where we choose

$(a_{1},b_{1},c_{1})=(\lambda_{a}a_{0},\lambda_{b}b_{0},\lambda_{c}c_{0})/N,$
(S94)

with

$N=\lambda_{a}a_{0}+\lambda_{b}b_{0}+\lambda_{c}c_{0}.$ (S95)

Since $N$ is just a proportionality factor, they will share the same maximums.
Therefore Eq. (S94) defines the “change of coordinates”. In order to transform
the boundary in Eq. (S90) to the “new coordinate system”, we need to invert
Eq. (S94). Using that $b=1-a-c$ in both coordinate systems, we find

$\begin{pmatrix}a_{0}\\\
c_{0}\end{pmatrix}=\frac{\lambda_{b}}{a_{1}\lambda_{c}(\lambda_{b}-\lambda_{a})+c_{1}\lambda_{a}(\lambda_{b}-\lambda_{c})+\lambda_{a}\lambda_{c}}\begin{pmatrix}\lambda_{c}a_{1}\\\
\lambda_{a}c_{1}\end{pmatrix}.$ (S96)

Plugging this into Eq. (S90) yields

$a_{1}^{2}=\frac{\lambda_{a}^{2}}{\lambda_{b}\lambda_{c}}8c_{1}(1-a_{1}-c_{1}).$
(S97)

Using Eq. (S85), we choose

$\displaystyle\lambda_{a}$ $\displaystyle=\frac{g_{\text{max}}}{16}dT^{2},$
$\displaystyle\lambda_{b}$ $\displaystyle=\frac{g_{\text{max}}}{8}dT^{2},$
$\displaystyle\lambda_{c}$ $\displaystyle=\frac{g_{\text{max}}}{16}dT^{2},$
(S98)

and solve for $a$. Retaining the correct solution, we find

$a=2(c-\sqrt{c}).$ (S99)

This corresponds to the black curve shown in Fig. 2b of the main text.

###  III.3 Slow driving regime

Here we optimize the trade-off between power, entropy production and power
fluctuations for the quantum dot engine in the slow-driving regime. For this
regime we split the protocol of the cycle into four steps: (i) we fix the bath
at a cold temperature $T_{\text{C}}$ and we slowly vary the dot energy
$\epsilon(t)\equiv u(t)E_{0}$ for a time $\tau_{C}$ from
$\epsilon(0)=\epsilon_{A}$ to $\epsilon(\tau_{C}^{-})=\epsilon_{B}$. (ii) We
now change the bath temperature to $T_{\text{H}}$ ($>T_{\text{C}}$), all while
we are changing the energy to $\epsilon_{B}T_{\text{H}}/T_{\text{C}}$ in such
a way that the populations of the dot levels remain constant all along the
process (because $[H_{u(t)},H_{u(t^{\prime})}]=0\;\forall\,t,t^{\prime}$).
This property allows us to perform this step arbitrarily fast without
affecting the objectives. (iii) while keeping the bath temperature at
$T_{\text{H}}$ we slowly vary the dot energy from
$\epsilon(\tau_{C}^{+})=\epsilon_{B}T_{\text{H}}/T_{\text{C}}$ to
$\epsilon(\tau^{-})=\epsilon_{A}T_{\text{H}}/T_{\text{C}}$ in a time
$\tau_{H}=\tau-\tau_{C}$. (iv) A second quench is performed (in the same
manner as step (ii)) to bring the temperature to $T_{\text{C}}$ and the dot
energy to $\epsilon_{A}$, which closes the cycle.

During each isotherm we can divide the heat exchange in the following manner:

$Q_{x}=T_{x}\Delta S_{x}-W_{x}^{\textrm{diss}},\quad x=(H,C)$ (S100)

where $T_{x}\Delta S_{x}$ is the reversible contribution and corresponds to
the quasistatic limit $\tau_{x}\rightarrow\infty$. The reversible term is
given by $\Delta S\equiv\Delta S_{H}=-\Delta
S_{C}=S(\pi_{0})-S(\pi_{\tau_{C}})$ with $S(\pi)$ the Von Neumann entropy.
Crucially, we will from now on assume the _slow driving_ regime (aka. _low
dissipation_). That is, we assume that the relaxation time scale $\gamma^{-1}$
of the system is much smaller than the time it takes to complete the protocol.
This allows us to expand the relevant quantities in orders of
$1/(\gamma\tau_{x})$ and keep only the leading terms. In this regime the state
can be written as

$\rho_{t}=\pi_{t}+\frac{1}{\tau\gamma}\rho^{(1)}_{t}+\mathcal{O}(\frac{1}{\gamma^{2}\tau^{2}}).$
(S101)

It is important to note that when the quenches of steps (ii) and (iv) are
performed the adiabatic term $\pi_{t}$ only contributes reversible work, but
the corrective term $\frac{1}{\tau\gamma}\rho^{(1)}_{t}$ will cause dissipated
work of order $\frac{1}{\tau\gamma}$ . For simplicity we want to neglect those
terms, in order to do so we can add a waiting time $\tau_{w}=\tau\Delta s$ at
the end of each isotherm. By using the same Lindbladian model as in the fast
driving regime (S43), it is clear that during this waiting time the corrective
term is exponentially suppressed. Therefore the total protocol time is now
$\tau=\tau_{C}+\tau_{H}+2\tau_{w}=\frac{\tau_{C}+\tau_{H}}{1-2\Delta s}$.

As a result of the slow driving expansion, for each isotherm the irreversible
terms can be written as [77, 26, 83, 125, 126, 127, 128, 129, 81, 130]:

$W_{x}^{\textrm{diss}}=\frac{T_{x}}{\tau_{x}}\int
ds~{}m_{\epsilon\epsilon}\dot{\epsilon}(s)^{2}\equiv
T_{x}\frac{\sigma_{x}}{\tau_{x}}+\mathcal{O}(\frac{1}{\gamma^{2}\tau_{x}^{2}}),$
(S102)

where $m_{\epsilon\epsilon}=\gamma^{-1}\partial^{2}\ln\Tr[e^{-\beta
H_{u}}]/\partial\epsilon^{2}$ is the thermodynamic metric. Furthermore, since
we only change the temperature of the bath without affecting its spectral
density we have that the optimal protocol verifies $\sigma_{H}=\sigma_{C}$
[26, 81]. We will therefore drop this index going forward. In a cycle the work
extracted is $W=Q_{C}+Q_{H}$, therefore the power of the engine reads

$\langle P\rangle=\frac{1}{\tau}\Delta T\Delta
S-\frac{\sigma}{\tau}\left(\frac{T_{\text{H}}}{\tau_{H}}+\frac{T_{\text{C}}}{\tau_{C}}\right),$
(S103)

with $\Delta=T_{\text{H}}-T_{\text{C}}$. The slow driving regime allows us to
use fluctuation dissipation relations [131, 132, 133, 134, 135, 82] to compute
the fluctuations of the work: in particular we have
$\frac{1}{2}\textrm{Var}(W_{x})=T_{x}W_{x}^{\textrm{diss}}$. Instead over the
quenches we get $\textrm{Var}(W_{y})=\Delta T^{2}C_{y}$, where
$C_{y}=-\beta^{2}\partial^{2}\ln\Tr[e^{-\beta H_{u_{y}}}]/\partial\beta^{2}$
is the heat capacity, for $y=A,B$. Therefore the power fluctuations of the
full cycle are

$\langle\Delta P\rangle=\frac{\textrm{Var}(W)}{\tau}=\frac{1}{\tau}\Delta
T^{2}\left(C_{A}+C_{B}\right)+\frac{2\sigma}{\tau}\left(\frac{T_{\text{H}}^{2}}{\tau_{H}}+\frac{T_{\text{C}}^{2}}{\tau_{C}}\right).$
(S104)

The entropy production rate is, by definition,

$\langle\Sigma\rangle=\frac{\sigma}{\tau}\left(\frac{1}{\tau_{H}}+\frac{1}{\tau_{C}}\right).$
(S105)

In this appendix we define a figure of merit where the normalization has been
absorbed in the weights

$G=a\langle P\rangle-\frac{b}{2}\langle\Delta
P\rangle-c\langle\Sigma\rangle=\frac{1}{\tau}\left[\Delta
TA-\sigma\left(\frac{\delta_{H}^{2}}{\tau_{H}}+\frac{\delta_{C}^{2}}{\tau_{C}}\right)\right]$
(S106)

with $A=a\Delta S-\frac{1}{2}b\Delta T(C_{A}+C_{B})$ and
$\delta_{x}^{2}=aT_{x}+bT_{x}^{2}+c$.

####  III.3.1 Optimization of the cycle

From (S106) we can see that, for given bath temperatures, $A$ depends only on
the boundary values $\epsilon_{A}$ and $\epsilon_{B}$ and $\sigma$ only
depends on the protocol’s geometric shape between those boundaries. Therefore
all the dependence on protocol time has been rendered explicit and we can
optimize over it by setting the derivatives with respect to $\tau_{C}$ and
$\tau_{C}$ to zero: we find

$\displaystyle\tau_{H}$
$\displaystyle=2\left(\delta_{H}+\delta_{C}\right)\frac{\delta_{H}\sigma}{\Delta
TA},$ (S107) $\displaystyle\tau_{C}$
$\displaystyle=2\left(\delta_{H}+\delta_{C}\right)\frac{\delta_{C}\sigma}{\Delta
TA},$ (S108) $\displaystyle\tau$
$\displaystyle=2\frac{\left(\delta_{H}+\delta_{C}\right)^{2}}{1-2\Delta
s}\frac{\sigma}{\Delta TA}.$ (S109)

It is important to note that, for the time optimization to be valid, we
require that $A$ is strictly positive. If $A$ is negative we can see that $G$
is also negative. But a cycle in which we do nothing results in $G=0$ which is
”better” than a cycle where $G$ is negative, therefore we can always assume
$A$ to be positive. We must now check the positivity of $A$. It is quite easy
to see that if $a=0$ then $A$ is negative and therefore the optimal cycle is
to do nothing. But if $a>0$ we then have that $A$ is strictly larger than zero
if and only if

$\frac{\Delta S}{C_{A}+C_{B}}>\frac{b}{2a}\Delta T.$ (S110)

It is clear that this equation can always be satisfied by choosing endpoints
of the protocol ($\epsilon_{A}$ and $\epsilon_{B}$) that have vanishing heat
capacity.
When we insert the optimal times into the figure of merit we find

$G=\frac{A^{2}\Delta T^{2}}{4\sigma}\frac{1-2\Delta
s}{(\delta_{H}+\delta_{C})^{2}}.$ (S111)

During the waiting time at the end of each isotherm, the decay towards the
thermal state goes as $e^{-\tau\gamma\Delta s}$, to be able to neglect the
correction we have to set $\tau\gamma\Delta s=r$, where $r$ is some arbitrary
$\mathcal{O}(1)$ constant that is sufficiently large. Therefore we find
$\Delta s=\frac{1}{2}\frac{rA}{rA+(\delta_{H}+\delta_{C})^{2}\sigma}$, which
leads to a new form of $G$:

$G=\frac{1}{4}\frac{A^{2}\Delta T^{2}}{rA+(\delta_{H}+\delta_{C})^{2}\sigma}.$
(S112)

At this point, the only term that depends on the protocol shape (given its
boundaries) is $\sigma$. As one would expect, $G$ is maximized when $\sigma$
is minimized. It has been shown [86, 136] that $\sigma$ is minimized by the
thermodynamic length
$\sigma_{min}=L^{2}=4\arccos\\!\left(\sqrt{p_{A}p_{B}}+\sqrt{(1-p_{A})(1-p_{B})}\right)^{2}$,
where $p_{A,B}=(1+e^{\beta_{C}\epsilon_{A,B}})^{-1}$ is the occupation
probability of the excited state at the endpoints of the isotherms.

Now we only have to optimize over the endpoints of the protocol the figure of
merit $G=\frac{1}{4}\frac{A^{2}\Delta
T^{2}}{rA+(\delta_{H}+\delta_{C})^{2}L^{2}}$. By taking derivatives with
respect to $\epsilon_{A}$ and $\epsilon_{B}$ we find two transcendental
equations that cannot be solved analytically. Equivalently, we can optimize
with respect to $p_{A}$ and $p_{B}$, though the same problem persists. But
these variables are more practical for numerical optimization as they are
defined on the finite range $[0,0.5]$. Therefore we do a numerical
maximization of $G$ for this last step, from the form of $G$ it can be seen
that the couple $(p_{A},p_{B})$ that maximizes $G$ depends only on two
parameters (assuming the temperatures are given): $\alpha:=b\Delta T/a$ and
$\alpha^{\prime}:=c/a$. We can note that a choice of these parameters defines
which trade-off of the objectives we want. Furthermore we have to satisfy the
constraint $\frac{\Delta S}{C_{A}+C_{B}}>\alpha/2$ (see Eq. (S110)). From a
numerical point of view, in order to plot the Pareto front, we are interested
in a range of values of $\alpha$ and $\alpha^{\prime}$ that range from
$10^{-2}$ to $10^{2}$. The optimization was done using a _Mathematica_ script
and is shown in Fig. S2, then we can use this to compute the power,
efficiency, and power fluctuations for different trade-off ratios $\alpha$ and
$\alpha^{\prime}$ with which we obtain the Pareto front shown in panel b of
Fig. 3 of the main text.

Figure S2: Probabilities of occupation of the excited state at the endpoints
of the isotherms ($A$ and $B$) that maximize $G$ in the slow-driving regime.
These are displayed as a function of $\alpha=b\Delta T/a$ and
$\alpha^{\prime}=c/a$ (parameters: $\beta_{C}=2$, $\beta_{H}=1$ and $r=2$).
These values are used to produce the Pareto front in panel b of Fig. 3 of the
main text.

####  III.3.2 Violation of thermodynamic uncertainty relation.

In the low power limit it is known that PDHE can violate thermodynamic
uncertainty relations [54, 55, 56, 57, 58, 59]. We reach this limit by not
optimizing for power but instead prioritizing the other two objectives. We
take the limit of low power by taking $a\ll b,c$; this implies that
$\alpha,\alpha^{\prime}\gg 1$, by Eq. (S110) this implies that we need
vanishing heat capacity at $A$ and $B$ (because $\Delta S$ is upper-bounded by
$\ln 2$). This is seen in Fig. S2: $p_{A}$ approaches the largest values it
can while $p_{B}$ approaches the smallest values it can, which both minimize
heat capacity. Furthermore, this also implies that the thermodynamic distance
is as large as it can be (since it is a measure of the distance between the
states at $A$ and $B$), but crucially it remains finite. By applying this to
Eq. (S109) it is clear that the magnitude of the protocol time is determined
by the ratio $\frac{L^{2}}{A}=\frac{1}{a}\frac{L^{2}}{\Delta
S-\alpha(C_{A}+C_{B})/2}$, it is thus clear that the protocol time diverges in
this limit.

The fact that the protocol time diverges in this limit also explains the fact
that for large $\alpha$ we see a loss of the dependence on $\alpha^{\prime}$
in Fig. S2: in this limit we don’t need to fix $\tau\gamma\Delta s$ to a
constant, we only need to fix $\Delta s$ to an arbitrarily low constant since
$\tau$ goes to infinity. This allows us to keep a simpler form of the figure
of merit: $G=\frac{A^{2}\Delta T^{2}}{4L^{2}}\frac{1-2\Delta
s}{(\delta_{H}+\delta_{C})^{2}}$, in which we can note that the optimization
over the protocol endpoints only has to be done over the ratio $A^{2}/L^{2}$;
that does not depend on $\alpha^{\prime}$. To maximize this ratio we can take
the derivatives with respect to $p_{A}$ and $p_{B}$. As before we end up with
transcendental equations:

$\displaystyle\frac{\alpha^{-1}\Delta S-(C_{A}+C_{B})/2}{\Delta TL}$
$\displaystyle=\beta_{C}\epsilon_{A}\sqrt{p_{A}(1-p_{A})}\left[\alpha^{-1}+1-\frac{1}{2}\beta_{C}\epsilon_{A}(1-2p_{A})\right],$
(S113) $\displaystyle\frac{\alpha^{-1}\Delta S-(C_{A}+C_{B})/2}{\Delta TL}$
$\displaystyle=\beta_{C}\epsilon_{B}\sqrt{p_{B}(1-p_{B})}\left[\alpha^{-1}-1+\frac{1}{2}\beta_{C}\epsilon_{B}(1-2p_{B})\right].$
(S114)

For $\alpha^{-1}=0$ we can see that $p_{A}=\frac{1}{2}$ and $p_{B}=0$ solves
the equations, which is coherent with the numerical results of the previous
section. We can now expand around $\alpha^{-1}=0$ by setting
$p_{A}=\frac{1}{2}+p_{A}^{(1)}\alpha^{-1}+\mathcal{O}(\alpha^{-2})$ and
$p_{B}=0+p_{B}^{(1)}\alpha^{-1}+\mathcal{O}(\alpha^{-2})$. We therefore find

$\displaystyle p_{A}$ $\displaystyle=\frac{1}{2}-\frac{\ln
2}{\pi}\alpha^{-1}+\mathcal{O}(\alpha^{-2}),$ (S115) $\displaystyle p_{B}$
$\displaystyle=\mathcal{O}(\alpha^{-2}).$ (S116)

With this result we can compute the values of the objectives analytically to
leading order in the low power limit:

$\displaystyle\langle P\rangle$ $\displaystyle=2\alpha^{-1}\left(\frac{\ln
2}{\pi}\right)^{2}\frac{\Delta
T^{2}}{\left(\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}+\sqrt{\frac{c}{b}+T_{\text{C}}^{2}}\right)^{2}}+\mathcal{O}(\alpha^{-2}),$
(S117) $\displaystyle\langle\Delta P\rangle$
$\displaystyle=2\alpha^{-2}\left(\frac{\ln 2}{\pi}\right)^{2}\frac{\Delta
T^{2}}{\left(\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}+\sqrt{\frac{c}{b}+T_{\text{C}}^{2}}\right)^{3}}\left(\frac{T_{\text{H}}^{2}}{\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}}+\frac{T_{\text{C}}^{2}}{\sqrt{\frac{c}{b}+T_{\text{C}}^{2}}}\right)+\mathcal{O}(\alpha^{-3}),$
(S118) $\displaystyle\langle\Sigma\rangle$
$\displaystyle=\alpha^{-2}\left(\frac{\ln
2}{\pi}\right)^{2}\frac{1}{\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}\sqrt{\frac{c}{b}+T_{\text{C}}^{2}}}~{}\frac{\Delta
T^{2}}{\left(\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}+\sqrt{\frac{c}{b}+T_{\text{C}}^{2}}\right)^{2}}+\mathcal{O}(\alpha^{-3}).$
(S119)

In particular we notice that $P\propto\alpha^{-1}$, $\Delta
P\propto\alpha^{-2}$, $\langle\Sigma\rangle\propto\alpha^{-2}$. Which results
in $\xi\propto\alpha^{2}$, therefore we can violate the thermodynamic
uncertainty relations arbitrarily. Therefore we can write, in this limit,
$\xi$ in terms of $\langle P\rangle$, $\langle\Delta P\rangle$ and
$\langle\Sigma\rangle$

$\displaystyle\xi$ $\displaystyle=\frac{16\left(\frac{\ln
2}{\pi}\right)^{4}P^{-2}\Delta
T^{4}}{\left(\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}+\sqrt{\frac{c}{b}+T_{c\vphantom{h}}^{2}}\right)^{3}}\frac{\left(\frac{c}{b}+T_{\text{H}}^{2}\right)\left(\frac{c}{b}+T_{c\vphantom{h}}^{2}\right)}{T_{\text{C}}^{2}\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}+T_{\text{H}}^{2}\sqrt{\frac{c}{b}+T_{c\vphantom{h}}^{2}}},$
(S120) $\displaystyle\xi$ $\displaystyle=\frac{4\left(\frac{\ln
2}{\pi}\right)^{2}\langle\dot{\Sigma}\rangle^{-1}\Delta
T^{2}}{\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}+\sqrt{\frac{c}{b}+T_{c\vphantom{h}}^{2}}}~{}\frac{\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}\sqrt{\frac{c}{b}+T_{c\vphantom{h}}^{2}}}{T_{\text{C}}^{2}\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}+T_{\text{H}}^{2}\sqrt{\frac{c}{b}+T_{c\vphantom{h}}^{2}}},$
(S121) $\displaystyle\xi$ $\displaystyle=8\left(\frac{\ln
2}{\pi}\right)^{2}\Delta P^{-1}\Delta
T^{2}\frac{\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}\sqrt{\frac{c}{b}+T_{c\vphantom{h}}^{2}}}{\left(\sqrt{\frac{c}{b}+T_{\text{H}}^{2}}+\sqrt{\frac{c}{b}+T_{c\vphantom{h}}^{2}}\right)^{2}}.$
(S122)

By maximising over the choice of $b$ and $c$ for a given $\alpha$ we can then
obtain the black lines in Fig. 4 of the main text.

###  III.4 Numerical Optimization using RL

Here we provide the training details and hyperparameters used to produce all
RL results presented in the main text in Figs. 2, 3, and 4 of the main text.
The method is described in Sec. I.3. A separate training was performed for all
values of the weights $(a,c)$ reported in Fig. 2 of the main text. Since the
cycles vary dramatically based on the choice of $(a,c)$, we mainly employed
two sets of hyperparameters, denoted with v1 and v2, shown in Table 1. The
hyperparameter names in Table 1 not defined in this paper are defined as in
Ref. [100], except for $n_{\text{train-steps}}$ that represents the number of
steps that were performed during training. For values of $a$ in $[0.4,1]$, we
used v1 and performed a single training. The only two exceptions are
$(a,c)=(1,0)$ and $(0.55,0)$, where we obtained better cycles training a
second time respectively with $n_{\text{train-steps}}=340k$ and $280k$.
Instead, training for values of $a$ in $[0.2,0.35]$ was less stable, since the
low power regime produces much longer cycles that require more “long-term
planning”. Furthermore, the maximum value of the figure of merit becomes
closer to $0$ as $a$ decreases. Since also “doing nothing” (i.e. setting
$\vec{u}(t)$ and $\beta(t)$ to a time-independent constant) produces a null
figure of merit, it becomes harder and harder for the RL method to distinguish
optimal cycles from these trivial solutions. To overcome these difficulties,
we mainly used v2, and we repeated the training up to $3$ times for some
values of $(a,c)$ \- choosing then the cycle with the largest figure of merit.
In total, we repeated $54$ trainings, some of which trained with v1. At last,
in some cases we trained for more steps, up to $n_{\text{train-steps}}=471k$
steps, and we set $\Delta t=0.5$ in some other trains, including the cases
with $c=0$ or $c=1-a$.

Hyperparameter name | v1 | v2
---|---|---
Hidden layers | 2 | ”
Hidden layers units | 256 | ”
Initial random steps | 6k | ”
First update at step | 1000 | ”
Batch size | 256 | ”
Learning rate | 0.001 | ”
$\bar{H}_{\text{start}}$ | 0.4 | 0.28
$\bar{H}_{\text{end}}$ | -7 | ”
$\bar{H}_{\text{decay}}$ | 108k | 162k
$n_{\text{updates}}$ | 50 | ”
$\rho_{\text{polyak}}$ | 0.995 | ”
$\mathcal{B}_{\text{size}}$ | 200k | ”
$\Delta t$ | 0.5 | 2
Discount factor $\gamma$ | 0.9997 | 0.9998
$n_{\text{train-steps}}$ | 240k | 360k
Table 1: Hyperparameters used to produce the results shown in the main text.
The columns v1 and v2 correspond to two different sets of hyperparameters used
in different regimes. The values not specified in v2 are the same as in v1.

###  III.5 Equivalence between SSHE and Fast-Otto PDHE for a two-level system

In this section we show that there is mapping between power, fluctuations, and
entropy production produced by a two-level system in a steady-state heat
engine (SSHE) and the same system periodically driven (PDHE) in a fast-Otto
cycle. This affinity can be explained, intuitively speaking, as in both cases
the state of the system relaxes to a fixed point (cf. [27]), and the work
exchanges are defined by a single variation of energy level (chemical
potential for the SSHE, energy gap quench for the PDHE, see details below).

This mapping explains why, among other things, the Fast-Otto cycle satisfies
the SSHE thermodynamic uncertainty relations, as shown in the main text.

Specifically, consider a two-level system connected to two thermal baths,
modeled via a simple rate master equation, that is

$\displaystyle\dot{\vec{p}}_{t}=\gamma_{i}(\vec{\pi}_{i}-\vec{p}_{t})$ (S123)

where $\gamma_{i}$ is the rate of bath $i$, $\vec{p}$ the populations vector,
and $\vec{\pi}_{i}:=(f_{i},1-f_{i})$, with
$f_{i}=\frac{1}{1+e^{-\beta_{i}\epsilon_{i}}}$. Notice that due to
normalization $\vec{p}=(p,1-p)$ and the dynamical equation can be equivalently
written as $\dot{p}_{t}=\gamma_{i}(f_{i}-p_{t})$. Both in the steady state
case and the fast driving cycle, the state tends to a fixed value

$\displaystyle p_{t}\rightarrow p\;.$ (S124)

It is possible to compute the steady state $p$, and the resulting power,
efficiency and fluctuations for such system when driven on a Fast-Otto cycle,
made by two quenches rapidly alternating. The analytic expressions can be
obtained as from Appendix II (cf. Eq. (S71), but we also allow different rates
for the two baths now).

The resulting expressions are given in Table 2 (right column), in terms of the
rates $\gamma_{i}$, the chosen gaps $\epsilon_{i}$, the fraction of time
$\theta_{i}$ spent on each bath $i$ ($\theta_{1}+\theta_{2}=1$). It turns out
that such expressions become formally equivalent to those of a SSHE based on a
two-level system in a chemical potential gradient between the two baths
$\mu_{i}$ (the gap $\epsilon$ of the qubit is in this case fixed), for which
the standard expressions in the left column of Table 2 can be obtained from
standard references. In particular, for $\Delta P$ it is possible to consider
Eq. (14) of [49] and substitute Eq. (24) of [45] for the values of the
$\mathcal{T}$-coefficients in the weak coupling limit.

By looking at Table 2, it is possible to see how each quantity is formally
equivalent between the two columns, via the following mapping

$\displaystyle\gamma_{i}$ $\displaystyle\rightarrow\theta_{i}\gamma_{i}\;,$
(S125) $\displaystyle\epsilon-\mu_{i}$
$\displaystyle\rightarrow\epsilon_{i}\;.$ (S126)

Finally, let us notice that the expression for the fluctuations $\Delta P$ can
equivalently be expressed in other forms, such as (we express it for the SSHE
case, although it can be expressed similarly for the Fast-Otto case, via the
above mapping)

$\displaystyle\Delta
P=\frac{2\gamma_{1}\gamma_{2}}{\gamma_{1}+\gamma_{2}}(\mu_{1}-\mu_{2})^{2}\left(f_{+}(1-f_{+})+(\bar{f}-f_{+})^{2}\right)\;,$
(S127)

where $f_{+}:=(f_{1}+f_{2})/2$. Such expression shows immediately the
positiveness of the fluctuations. Another equivalent expression, which is
typically obtained in SSHE calculations is

$\displaystyle\Delta
P=\frac{\gamma_{1}\gamma_{2}}{\gamma_{1}+\gamma_{2}}(\mu_{1}-\mu_{2})^{2}\left(f_{1}(1-f_{2})+f_{2}(1-f_{1})\right)-\frac{2\gamma_{1}^{2}\gamma_{2}^{2}}{(\gamma_{1}+\gamma_{2})^{3}}(\mu_{1}-\mu_{2})^{2}\left(f_{1}-f_{2}\right)^{2}\;.$
(S128)

| SSHE | Fast-Otto PDHE
---|---|---
$f_{i}$ | $(1+e^{-\beta_{i}(\epsilon-\mu_{i})})^{-1}$ | $(1+e^{-\beta_{i}\epsilon_{i}})^{-1}$
$p$ | $\bar{f}:=\dfrac{\gamma_{1}f_{1}+\gamma_{2}f_{2}}{\gamma_{1}+\gamma_{2}}$ | $\bar{f}:=\dfrac{\theta_{1}\gamma_{1}f_{1}+\theta_{2}\gamma_{2}f_{2}}{\theta_{1}\gamma_{1}+\theta_{1}\gamma_{2}}$
$P$ | $\dfrac{\gamma_{1}\gamma_{2}}{\gamma_{1}+\gamma_{2}}(f_{1}-f_{2})(\mu_{1}-\mu_{2})$ | $\dfrac{\theta_{1}\theta_{2}\gamma_{1}\gamma_{2}}{\theta_{1}\gamma_{1}+\theta_{2}\gamma_{2}}(f_{1}-f_{2})(\epsilon_{2}-\epsilon_{1})$
$\eta$ | $1-\dfrac{\epsilon-\mu_{1}}{\epsilon-\mu_{2}}$ | $1-\dfrac{\epsilon_{1}}{\epsilon_{2}}$
$\Delta P$ | $2\frac{\gamma_{1}\gamma_{2}}{\gamma_{1}+\gamma_{2}}(\epsilon_{1}-\epsilon_{2})^{2}\left[\bar{f}(1-\bar{f})+\left(\frac{f_{1}+f_{2}}{2}-\bar{f}\right)(1-2\bar{f})\right]$ | $2\frac{\gamma_{1}\theta_{1}\gamma_{2}\theta_{2}}{\gamma_{1}\theta_{1}+\gamma_{2}\theta_{2}}(\epsilon_{1}-\epsilon_{2})^{2}\left[\bar{f}(1-\bar{f})+\left(\frac{f_{1}+f_{2}}{2}-\bar{f}\right)(1-2\bar{f})\right]$
Table 2: Mapping between a steady-state engine and a fast-Otto cycle for the
case of a two-level system.