# Model-free optimization of power/efficiency tradeoffs in quantum thermal
machines using reinforcement learning

Paolo A. Erdman<EMAIL_ADDRESS>Freie Universität Berlin, Department of
Mathematics and Computer Science, Arnimallee 6, 14195 Berlin, Germany  Frank
Noé<EMAIL_ADDRESS>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

A quantum thermal machine is an open quantum system that enables the
conversion between heat and work at the micro or nano-scale. Optimally
controlling such out-of-equilibrium systems is a crucial yet challenging task
with applications to quantum technologies and devices. We introduce a general
model-free framework based on Reinforcement Learning to identify out-of-
equilibrium thermodynamic cycles that are Pareto optimal trade-offs between
power and efficiency for quantum heat engines and refrigerators. The method
does not require any knowledge of the quantum thermal machine, nor of the
system model, nor of the quantum state. Instead, it only observes the heat
fluxes, so it is both applicable to simulations and experimental devices. We
test our method on a model of an experimentally realistic refrigerator based
on a superconducting qubit, and on a heat engine based on a quantum harmonic
oscillator. In both cases, we identify the Pareto-front representing optimal
power-efficiency tradeoffs, and the corresponding cycles. Such solutions
outperform previous proposals made in the literature, such as optimized Otto
cycles, reducing quantum friction.

## Introduction

A driving force of the research field of quantum thermodynamic is the quest of
understanding and designing quantum thermal machines (QTMs), i.e. devices that
convert between heat and work at the micro or nanoscale exploiting quantum
effects [1, 2, 3, 4, 5]. Such devices could be operated as heat engines, which
convert heat into work, or refrigerators, that extract heat from a cold bath.
Recent experiments have measured the heat flowing across these devices [6, 7,
8, 9], and early experimental realizations of QTMs have been reported [10, 11,
12, 13, 14, 15, 16, 17].

However, the optimal control of such devices, necessary to reveal their
maximum performance, is an extremely challenging task that could find
application in the control of quantum technologies and devices beyond QTMs.
The difficulties include: (i) having to operate in finite time, the state can
be driven far from equilibrium, where the thermal properties of the system are
model-specific; (ii) the optimization is a search over the space of all
possible time-dependent controls, which increases exponentially with the
number of time points describing the cycle; (iii) in experimental devices,
often subject to undesired effects such as noise and decoherence [18], we
could have a limited knowledge of the actual model describing the dynamics of
the QTM.

A further difficulty (iv) arises in QTMs, since the maximization of their
performance requires a multi-objective optimization. Indeed, the two main
quantities that describe the performance of a heat engine (refrigerator) are
the extracted power (cooling power) and the efficiency (coefficient of
performance). The optimal strategy to maximize the efficiency consists of
performing reversible transformations [19] which are, however, infinitely
slow, and thus deliver vanishing power. Conversely, maximum power is typically
reached at the expense of reduced efficiency. Therefore, one must seek optimal
trade-offs between the two.

The theoretical optimization of QTMs is typically carried out making
restrictive assumptions on the cycle. For example, optimal strategies have
been derived assuming the driving speed of the control to be slow [20, 21, 22,
23, 24, 25, 26, 27, 28, 29] or fast [30, 31, 32] compared to the
thermalization time. Other approaches consists of assuming a-priori a specific
shape of the cycle structure [33, 34, 35, 36, 37, 38, 39, 40], such as the
Otto cycle [41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56].
Shortcuts to adiabaticity [57, 58, 59, 60, 61, 62, 63, 64, 65] and variational
strategies [66, 67, 68] have also been employed.

In general, aside from variational approaches, there is no guarantee that
these regimes and cycles are optimal. Recently, reinforcement-learning (RL)
has been used to find cycles that maximize the power of QTMs without making
assumptions on the cycle structure [69]. However, this approach requires a
model of the system and the knowledge of the quantum state of the system,
which restricts its practical applicability. This calls for the development of
robust and general strategies that overcome all above-mentioned difficulties
(i-iv).

We propose a RL-based method with the following properties: (i) it finds
cycles yielding near Pareto-optimal trade-offs between power and efficiency,
i.e. the collection of cycles such that it is not possible to further improve
either power or efficiency, without decreasing the other one. (ii) It only
requires the heat currents as input, and not the quantum state of the system.
(iii) It is completely model-free. (iv) It does not make any assumption on the
cycle structure, nor on the driving speed. The RL method is based on the Soft
Actor-Critic algorithm [70, 71], introduced in the context of robotics and
video-games [72, 73], generalized to combined discrete and continuous actions
and to optimize multiple objectives. RL has received great attention for its
success at mastering tasks beyond human-level such as playing games [74, 75,
76], and for robotic applications [77]. RL has been recently used for quantum
control [78, 79, 80, 81, 82, 83, 84, 85, 86, 87], outperforming previous
state-of-the-art methods [88, 89], for fault-tolerant quantum computation [90,
91], and to minimize entropy production in closed quantum systems [92].

We prove the validity of our approach finding the full Pareto-front, i.e. the
collection of all Pareto-optimal cycles describing optimal power-efficency
tradeoffs, in two paradigmatic systems that have been well studied in
literature: a refrigerator based on an experimentally realistic
superconducting qubit [49, 6], and a heat engine based on a quantum harmonic
oscillator [43]. In both cases we find elaborate cycles that outperform
previous proposal mitigating quantum friction [93, 43, 49, 94, 95, 66, 55],
i.e. the detrimental effect of the generation of coherence in the
instantaneous eigenbasis during the cycle. Remarkably, we can also match the
performance of cycles found with the RL method of Ref. [69] that, as opposed
to our model-free approach, requires monitoring the full quantum state and
only optimizes the power.

Figure 1: Schematic representation of a quantum thermal machine controlled by
a computer agent. A quantum system (gray circle) can be coupled to a hot
(cold) bath at inverse temperature $\beta_{\text{H}}$ ($\beta_{\text{C}}$),
represented by the red (blue) square, enabling a heat flux $J_{\text{H}}(t)$
($J_{\text{C}}(t)$). The quantum system is controlled by the computer agent
through a set of experimental control parameters $\vec{u}(t)$, such as an
energy gap or an oscillator frequency, that control the power exchange $P(t)$,
and through a discrete control $d(t)=\\{\text{Hot},\text{Cold},\text{None}\\}$
that determines which bath is coupled to the quantum system.

### Setting: Black-box Quantum Thermal Machine

We describe a QTM by a quantum system, acting as a “working medium”, that can
exchange heat with a hot (H) or cold (C) thermal bath characterized by inverse
temperatures $\beta_{\text{H}}<\beta_{\text{C}}$ (Fig. 1). Our method can be
readily generalized to multiple baths, but we focus the description on two
baths here.

We can control the evolution of the quantum system and exchange work with it
through a set of time-dependent continuous control parameters $\vec{u}(t)$
that enter in the Hamiltonian $H[\vec{u}(t)]$ of the quantum system [96], and
through a discrete control $d(t)=\\{\text{Hot},\text{Cold},\text{None}\\}$
that determines which bath is coupled to the system. $J_{\text{H}}(t)$ and
$J_{\text{C}}(t)$ denote the heat flux flowing out respectively from the hot
and cold bath at time $t$.

Our method only relies on the following two assumptions:

  1. 1.

the RL agent can measure the heat fluxes $J_{\text{C}}(t)$ and
$J_{\text{H}}(t)$ (or their averages over a time period $\Delta t$);

  2. 2.

$J_{\text{C}}(t)$ and $J_{\text{H}}(t)$ are functions of the control history
${(\vec{u}(t-T),d(t-T)),...,(\vec{u}(t),d(t))}$, where $T$ is the timescale
over which the QTM remembers past controls.

In contrast to previous work [69], the RL optimization algorithm does not
require any knowledge of the microscopic model of the inner workings of the
quantum system, nor of its quantum state; it is only provided with the values
of the heat fluxes $J_{\text{C}}(t)$ and $J_{\text{H}}(t)$. These can be
either computed from a theoretical simulation of the QTM [69], or measured
directly from an experimental device whenever the energy change in the heat
bath can be monitored without influencing the energetics of the quantum system
(see e.g. experimental demonstrations [6, 7, 8, 9]). In this sense, our
quantum system is treated as a “black-box”, and our RL method is “model-free”.
Any theoretical model or experimental device satisfying these requirements can
be optimized by our method, including also classical stochastic thermal
machines. The timescale $T$ is finite because of energy dissipation and
naturally emerges by making the minimal assumption that the coupling of the
quantum system to the thermal baths drives the system towards a thermal state
within some timescale $T$. Such a timescale can be rigorously identified e.g.
within the weak system-bath coupling regime, and in the reaction coordinate
framework that can describe non-Markovian and strong-coupling effects [97]. In
a Markovian setting, $T$ is related to the inverse of the characteristic
thermalization rate.

The thermal machines we consider are the heat engine and the refrigerator. Up
to an internal energy contribution that vanishes after each repetition of the
cycle, the instantaneous power of a heat engine equals the extracted heat:

$P_{\text{heat}}(t)=J_{\text{C}}(t)+J_{\text{H}}(t),$ (1)

and the cooling power of a refrigerator is:

$P_{\text{cool}}(t)=J_{\text{C}}(t).$ (2)

The entropy production is given by

$\Sigma(t)=-\beta_{\text{C}}J_{\text{C}}(t)-\beta_{\text{H}}J_{\text{H}}(t),$
(3)

where we neglect the contribution of the quantum system’s entropy since it
vanishes after each cycle.

### Machine Learning Problem

Our goal is to identify cycles, i.e. periodic functions $\vec{u}(t)$ and
$d(t)$, that maximize a trade-off between power and efficiency on the long
run. Since power and efficiency cannot be simultaneously optimized, we use the
concept of Pareto-optimality [98, 99]. Pareto-optimal cycles are those where
power or efficiency cannot be further increased without sacrificing the other
one. The Pareto-front, defined as the collection of power-efficiency values
delivered by all Pareto-optimal cycles, represents all possible optimal trade-
offs. To find the Pareto-front, we define the reward function $r_{c}(t)$ as:

$r_{c}(t)=c\frac{P(t)}{P_{0}}-(1-c)\frac{\Sigma(t)}{\Sigma_{0}},$ (4)

where $P(t)$ is the power of a heat engine (Eq. 1) or cooling power of a
refrigerator (Eq. 2), $P_{0}$, $\Sigma_{0}$ are reference values to normalize
the power and entropy production, and $c\in[0,1]$ is a weight that determines
the trade-off between power and efficiency. As in Ref. [69], we are interested
in cycles that maximize the long-term performance of QTMs; we thus maximize
the return $\expectationvalue*{r_{c}}(t)$, where $\langle\cdot\rangle(t)$
indicates the exponential moving average of future values:

$\expectationvalue*{r_{c}}(t)=\kappa\int_{0}^{\infty}e^{-\kappa\tau}r_{c}(t+\tau)\,d\tau.$
(5)

Here $\kappa$ is the inverse of the averaging timescale, that will in practice
be chosen much longer than the cycle period, such that
$\expectationvalue*{r_{c}}(t)$ is approximately independent of $t$.

For $c=1$, we are maximizing the average power
$\expectationvalue*{r_{1}}=\expectationvalue*{P}/P_{0}$. For $c=0$, we are
minimizing the average entropy production
$\expectationvalue*{r_{0}}=-\expectationvalue{\Sigma}/\Sigma_{0}$, which
corresponds to maximizing the efficiency. For intermediate values of $c$, the
maximization of $\expectationvalue*{r_{c}}$ describes trade-offs between power
and efficiency (see “Optimizing the entropy production” in Materials and
Methods for details). Interestingly, if convex, it has been shown that the
full Pareto-front can be identified repeating the optimization of
$\expectationvalue{r_{c}}$ for many values of $c$ [98, 100].

## Results

Figure 2: (A): Schematic representation of the learning process. A computer
agent (blue box) chooses an action $a_{i}$ at time-step $i$ based on the
current state $s_{i}$ of the QTM (gray box) through the policy function
$\pi(a_{i}|s_{i})$. The action, that encodes the control ($\vec{u}(t_{i})$,
$d(t_{i})$), is passed to the QTM (lower arrow). The new state $s_{i+1}$,
composed of the time-series of the last $N$ actions, and the reward $r_{i+1}$
are returned to the agent (upper arrow), which uses this information to
improve $\pi(a|s)$ using the soft actor-critic algorithm, which learns also
the values function $Q(s,a)$. This process is reiterated until convergence of
the policy. (B-C): Schematic representation of the NN architectures used to
parameterize the policy (B) and the value function (C). The action time-series
in $s_{i}$ is processed using multiple 1D convolution blocks, each one halving
the length of the series. The final output is produced by fully connected
(f.c.) layers.

### Deep reinforcement learning for black-box quantum thermal machines

In RL, a computer agent must learn to master some task by repeated
interactions with some environment. Here we develop an RL approach where the
agent maximizes the return (5) and the environment is the QTM with its
controls (Fig. 2A). To solve the RL problem computationally, we discretize
time as $t_{i}=i\Delta t$. By time-discretizing the return (5), we obtain a
discounted return whose discount factor $\gamma=\exp(-\kappa\Delta t)$
determines the averaging timescale and expresses how much we are interested in
future or immediate rewards (see “Reinforcement Learning Implementation” in
Materials and Methods for details).

At each time step $t_{i}$, the agent employs a policy function $\pi(a|s)$ to
choose an action $a_{i}=\\{\vec{u}(t_{i}),d(t_{i})\\}$ based on the state
$s_{i}$ of the environment. Here, the policy function $\pi(a|s)$ represents
the probability of choosing action $a$, given that the environment is in state
$s$, $\vec{u}(t)$ are the continuous controls over the quantum system, and
$d(t_{i})\in\\{\text{Hot},\text{Cold},\text{None}\\}$ is a discrete control
that selects the bath the system is coupled to. All controls are considered to
be constant during time step of duration $\Delta t$. The aim of RL is to learn
an optimal policy function $\pi(a|s)$ that maximizes the return.

In order to represent a black-box quantum system whose inner mechanics are
unknown, we define the control history during a time interval of length $T$ as
the observable state:

$s_{i}=(a_{i-N},a_{i-N+1},\dots,a_{i-1}),$ (6)

where $N=T/\Delta t$. Therefore, the state of the quantum system is implicitly
defined by the sequence of the agent’s $N$ recent actions.

To find an optimal policy we employ the soft actor-critic algorithm, that
relies on learning also a value function $Q(s,a)$, generalized to a
combination of discrete and continuous actions [70, 71, 72, 73]. The policy
function $\pi(a|s)$ plays the role of an “actor” that chooses the actions to
perform, while a value function $Q(s,a)$ plays the role of a “critic” that
judges the choices made by the actor, thus providing feedback to improve the
actor’s behavior. We further optimize the method for a multi-objective setting
by introducing a separate critic for each objective, i.e. one value function
for the power, and one for the entropy production. This allow us to vary the
weight $c$ during training, thus enhancing convergence (see “Reinforcement
Learning Implementation” in Materials and Methods for details).

We learn the functions $\pi(a|s)$ and $Q(s,a)$ using a deep NN architecture
inspired by WaveNet, an architecture that was developed for processing audio
signals [101] (See Figs. 2B-C). We introduce a “convolution block” to
efficiently process the time-series of actions defining the state $s_{i}$. It
consists of a 1D convolution with kernel size and stride of $2$, such that it
halves the length of the input. It is further equipped with a residual
connection to improve trainability [102] (see “Reinforcement Learning
Implementation” in Materials and Methods for details). The policy
$\pi(a_{i}|s_{i})$ is described by a NN that takes the state $s_{i}$ as input,
and outputs parameters $\mu$ and $\sigma$ describing the probability
distribution from which action $a_{i}$ is sampled (Fig. 2B). The value
function $Q(s_{i},a_{i})$ is computed by feeding $(s_{i},a_{i})$ into a NN,
and outputting $Q(s_{i},a_{i})$ (Fig. 2C). Both $\pi(a_{i}|s_{i})$ and
$Q(s_{i},a_{i})$ process the state by feeding it through multiple convolution
blocks (upper orange boxes in Figs. 2B and 2C), each one halving the length of
the time-series, such that the number of blocks and of parameters in the NN is
logarithmic in $N$. Then a series of fully-connected layers produce the final
output.

The policy and value functions are determined by minimizing the loss functions
in Eqs. (39) and (49) using the ADAM optimization algorithm [103]. The
gradient of the loss functions is computed off-policy, over a batch of past
experience recorded in a replay buffer, using back-propagation (see
“Reinforcement Learning Implementation” in Materials and Methods for details).

### Pareto-optimal cycles for a superconducting qubit refrigerator

We first consider a refrigerator based on an experimentally realistic system:
a superconducting qubit coupled to two resonant circuits that behave as heat
baths [49] (Fig. 3A). Such a system was experimentally studied in the steady-
state in Ref. [6]. The system Hamiltonian is given by [49, 55, 63]:

$\hat{H}[u(t)]=-E_{0}\left[\Delta\hat{\sigma}_{x}+u(t)\hat{\sigma}_{z}\right],$
(7)

where $E_{0}$ is a fixed energy scale, $\Delta$ characterizes the minimum gap
of the system, and $u(t)$ is our control parameter. In this setup the coupling
to the baths, described by the commonly employed Markovian master equation
[104, 105, 106, 107], is fixed, and cannot be controlled. However, the qubit
is resonantly coupled to the baths at different energies. The $u$-dependent
coupling strength to the cold (hot) bath is described by the function
$\gamma^{(\text{C})}_{u}$ ($\gamma^{(\text{H})}_{u}$), respectively (Fig. 3F).
As in Ref. [63], the coupling strength is, respectively, maximal at $u=0$
($u=1/2$), with a resonance width determined by the “quality factor”
$Q_{\text{C}}$ ($Q_{\text{H}}$) (see “Physical model” in Materials and Methods
for details). This allows us to choose which bath is coupled to the qubit by
tuning $u(t)$.

Figure 3:  Training of the superconducting qubit refrigerator model to
optimize $\expectationvalue*{r_{c}}$ at $c=0.6$. (A): Schematic representation
of the energy levels of the qubit (horizontal black lines) that are controlled
by $u(t)$. The gray arrow represents the input power, while the colored arrows
represent the heat fluxes. (B): Return $\expectationvalue*{r_{c}}_{i}$
computed over past rewards (black curve), running average of the cooling power
$\expectationvalue*{P_{\text{cool}}}_{i}/P_{0}$ (green curve), and of the
negative entropy production $-\expectationvalue*{\Sigma}_{i}/\Sigma_{0}$
(orange curve), as a function of the training step. The dashed line represents
the value of the return found optimizing the period of a smoothed trapezoidal
cycle. (C): Value of the weight $c$ as a function of the step. It is varied
during training from $1$ to the final value $0.6$ to improve convergence. (D):
Actions chosen by the agent, represented by the value of $u$, as a function of
step, zoomed around the three black circles in panel (B). (E): Final
deterministic cycle found by the agent (thick black dots) and smoothed
trapezoidal cycle (thin dashed line) whose return is given by the dashed line
in panel (B), as a function of time. (F): coupling strength
$\gamma^{\text{(C)}}_{u}$ (blue curve) and $\gamma^{\text{(H)}}_{u}$ (red
curve) as a function of $u$ (on the y-axis). The parameters used for training
are $N=128$, $g_{\mathrm{H}}=g_{\mathrm{C}}=1$, $\beta_{\mathrm{H}}=10/3$,
$\beta_{\mathrm{C}}=2\beta_{\mathrm{H}}$, $Q_{\mathrm{H}}=Q_{\mathrm{C}}=4$,
$E_{0}=1$, $\Delta=0.12$, $\omega_{\mathrm{H}}=1.028$,
$\omega_{\mathrm{C}}=0.24$, $\mathcal{U}=[0,0.75]$, $\Delta t=0.98$,
$\gamma=0.997$, $P_{0}=6.62\cdot 10^{-4}$ and $\Sigma_{0}=0.037$.

In Fig. 3 we show an example of our training procedure to optimize the return
$\expectationvalue{r_{c}}$ at $c=0.6$ using $N=128$ steps determining the RL
state, and varying $c$ during training from $1$ to $0.6$ (Fig. 3C). In the
early stages of the training, the return $\expectationvalue*{r_{c}}_{i}$,
computed as in Eq. (28) but over past rewards, and the running averages of the
cooling power $\expectationvalue*{P_{\text{cool}}}_{i}$ and of the negative
entropy production $-\expectationvalue*{\Sigma}_{i}$ all start off negative
(Fig. 3B), and the corresponding actions are random (left panel of Fig. 3D).
Indeed, initially the RL agent has no experience controlling the QTM, so
random actions are performed, resulting in heating the cold bath, rather than
cooling it, and in a large entropy production. However, with increasing steps,
the chosen actions exhibit some structure (Fig. 3D), and the return
$\expectationvalue*{r_{c}}_{i}$ increases (Fig. 3B). While both the power and
the negative entropy production initially increase together, around step 100k
we see that $-\expectationvalue*{\Sigma}_{i}$ begins to decrease. This is a
manifestation of the fact that power and entropy production cannot be
simultaneously optimized. Indeed, the agent learns that in order to further
increase the return, it must “sacrifice” some entropy production to produce a
positive and larger cooling power. In fact, the only way to achieve positive
values of $\expectationvalue*{r_{c}}_{i}$ is to have a positive cooling power,
which inevitably requires producing entropy. Eventually all quantities in Fig.
3B reach a maximum value, and the corresponding final deterministic cycle
(i.e. the cycle generated by policy switching off stochasticity, see
“Reinforcement Learning Implementation” in Materials and Methods for details)
is shown in Fig. 3E as thick black dots.

For the same system, Ref. [63] proposed a smoothed trapezoidal cycle $u(t)$
oscillating between the resonant peaks at $u=0$ and $u=1/2$ and optimized the
cycle time (Fig. 3E, dashed line). While this choice outperformed a sine and a
trapezoidal cycle [49], the cycle found by our RL agent produces a larger
return (Fig. 3B). The optimal trapezoidal cycle found for $c=0.6$ is shown in
Fig. 3E as a dashed line (see “Comparing with other methods” in Materials and
Methods for details).

Figure 4:  Results for the optimization of the superconducting qubit
refrigerator model. (A): final value of the return
$\expectationvalue*{r_{c}}$, as a function of $c$, found using the RL method
(black and blue points), and optimizing the period of a trapezoidal cycle (red
dots). The error bars represent the standard deviation of the return computed
over $5$ independent training runs. (B): corresponding values of the final
average cooling power $\expectationvalue*{P_{\text{cool}}}$ and of the
coefficient of performance $\eta_{\text{cool}}$ found using the RL method
(black and blue dots), optimizing the trapezoidal cycle (red dots), and using
the RL method of Ref. [69] (purple cross). Results for each of the $5$
repetitions are shown as separate points to visualize the variability across
multiple trainings. (C-F): final deterministic cycles identified by the RL
method (thick black dots), as a function of time, corresponding to the blue
points in panels (A) and (B) (respectively for $c=1,0.8,0.6,0.4$ choosing the
training run with the largest return). The dashed line represents the
trapezoidal cycle that maximizes the return for the same value of $c$ [not
shown in panel (F) since no cycle yields a positive return]. The parameters
used for training are chosen as in Fig. 3.

Fig. 4 compares optimal cycles for different trade-offs between cooling power
and coefficient of performance $\eta_{\text{cool}}$, the latter defined as the
ratio between the average cooling power, and the average input power. This is
achieved by repeating the optimization for various values of $c$. To
demonstrate the robustness of our method, the optimization of
$\expectationvalue{r_{c}}$ was repeated $5$ times for each choice of $c$
(variability shown with error bars in Fig.4A, and as separate points in
Fig.4B). The RL method substantially outperforms the trapezoidal cycle by
producing larger final values of the return $\expectationvalue*{r_{c}}$ at all
values of $c$ (Fig. 4A), and by producing a better Pareto front (Fig. 4B). The
RL cycles simultaneously yield higher power by more than a factor of 10, and a
larger $\eta_{\text{cool}}$, for any choice of the power-efficiency trade-off.
The model-free RL cycles can also deliver the same power at a substantially
higher COP (roughly $10$ times larger) as compared to the cycle found with the
RL method of Ref. [69], which only optimizes the power. This is remarkable
since, as opposed to the current model-free method, the method in Ref. [69]
has access to the full quantum state of the system, and not only to the heat
currents (see “Comparing with other methods” in Materials and Methods for
details). This also shows that a large efficiency improvement can be achieved
by sacrificing very little power.

As expected, the period of the RL cycles increases as $c$ decreases and the
priority shifts from high power to high $\eta_{\text{cool}}$ (Figs. 4C-F,
black dots). However, the period is much shorter than the corresponding
optimized trapezoidal cycle (dashed line), and the optimal control sequence is
quite unintuitive, even going beyond the resonant point at $u=1/2$. As argued
in [49, 55, 63], the generation of coherence in the instantaneous eigenbasis
of the quantum system, occurring because $[\hat{H}(u_{1}),\hat{H}(u_{2})]\neq
0$ for $u_{1}\neq u_{2}$, causes power losses that increase with the speed of
the cycle. We find that we can interpret the power enhancement achieved by our
cycle as a mitigation of such detrimental effect: indeed, we find that
trapezoidal cycles operated at the same frequency as the RL cycle generate
twice as much coherence as the RL cycles (see “Generation of coherence” in
Materials and Methods for details). In either case, cycles with higher power
tend to generate more coherence.

Given the stochastic nature of RL, we also compared the cycles obtained across
the $5$ independent training runs, finding that cycles are typically quite
robust, displaying only minor changes (see Fig. 8 of Methods for four cycles
found in independent training runs corresponding to Figs. 4C-F).

### Pareto-optimal cycles for a quantum harmonic oscillator engine

We now consider a heat engine based on a collection of non-interacting
particles confined in a harmonic potential [43] (Fig. 5A). The Hamiltonian is
given by

$\hat{H}[u(t)]=\frac{1}{2m}\hat{p}^{2}+\frac{1}{2}m(u(t)w_{0})^{2}\hat{q}^{2},$
(8)

where $m$ is the mass of the system, $w_{0}$ is a reference frequency and
$\hat{p}$ and $\hat{q}$ are the momentum and position operators. The control
parameter $u(t)$ allows us to change the frequency of the oscillator. Here, at
every time step we let the agent choose which bath (if any) to couple to the
oscillator. The coupling to the baths, characterized by the thermalization
rates $\Gamma_{\alpha}$, is modeled using the Lindblad master equation as in
Ref. [43] (see “Physical model” in Materials and Methods for details). In
contrast to the superconducting qubit case, $c$ is held constant during
training.

Figure 5:  Results for the optimization of the harmonic oscillator heat engine
model. (A): Schematic representation of the energy levels of the particles
(black horizontal lines) trapped in a harmonic potential (parabolic curve)
whose amplitude is controlled by $u(t)$. The gray arrow represents the
extracted power, while the colored arrows represent the heat fluxes. (B):
final value of $\expectationvalue*{r_{c}}$, as a function of $c$, found using
the RL method (black and blue dots), and optimizing the Otto cycle (red dots).
The error bars represent the standard deviation of the return computed over
$5$ independent training runs. (C): corresponding values of the average power
$\expectationvalue*{P_{\text{heat}}}/P_{0}$ and of the efficiency
$\eta_{\text{heat}}$ found using the RL method (black and blue dots),
optimizing the Otto cycle (red dots), and using the RL method of Ref. [69]
(purple cross). Results for each of the $5$ repetitions are shown as separate
points to visualize the variability across multiple trainings. (D-E): final
deterministic cycle identified by the RL method (thick dots), as a function of
time, corresponding to the blue points in panels (B) and (C) (respectively
$c=1,0.5$ choosing the training run with the largest return). The color
corresponds to the discrete choice
$d=\\{\text{Hot},\text{Cold},\text{None}\\}$ (see legend). The dashed line
represents the Otto cycle that maximizes the return for the same value of $c$.
The parameters used for training are $N=128$,
$\Gamma^{(\text{H})}=\Gamma^{(\text{C})}=0.6$, $\beta_{\text{H}}=0.2$,
$\beta_{\text{C}}=2$, $w_{0}=2$, $\mathcal{U}=[0.5,1]$ (to enable a fair
comparison with Ref. [43]), $\Delta t=0.2$, $\gamma=0.999$, $P_{0}=0.175$ and
$\Sigma_{0}=0.525$.

Fig. 5 reports the results on the optimal trade-offs between extracted power
and efficiency $\eta_{\text{heat}}$, the latter defined as the ratio between
the extracted power and the input heat, in the same style of Fig. 4. In this
setup, we compare our RL-based results to the well-known Otto cycle. The
authors of Ref. [43] study this system by optimizing the switching times of an
Otto cycle, i.e. the duration of each of the 4 segments, shown as a dashed
lines in Figs. 5D-E, composing the cycle (see “Comparing with other methods”
in Materials and Methods for details).

The RL method produces cycles with a larger return and with a better power-
efficiency Pareto-front with respect to the Otto cycle (Fig. 5B,C). The cycles
found by the RL method significantly outperforms the Otto engine in terms of
delivered power. For $c=1$, a high-power cycle is found (Fig. 5D and
corresponding blue dots in Figs. 5B-C) but at the cost of a lower efficiency
than the Otto cycles. However, at $c=0.5$, the RL method finds a cycle that
matches the maximum efficiency of the Otto cycles, while delivering a $\sim
30\%$ higher power (Fig. 5E and corresponding blue dots in Figs. 5B-C)
Remarkably, our model-free RL method also finds cycles with nearly the same
power as the RL method of Ref. [69], but at almost twice the efficiency (see
“Comparing with other methods” in Materials and Methods for details). As in
Fig. 4, we see that a very small decrease in power can lead to a large
efficiency increase.

Interestingly, as shown in Figs. 5D-E, the cycles found by the RL agent share
many similarities with the Otto cycle: both alternate between the hot and cold
bath (orange and blue portions) with a similar period. However, there are some
differences: at $c=1$, the RL cycle ramps the value of $u$ while in contact
with the bath, eliminating the unitary stroke (Fig. 5D). Instead, at $c=0.5$,
the RL agent employs a unitary stroke that is quite different respect to a
linear ramping of $u$ (Fig. 5E, green dots). As in the superconducting qubit
case, the enhanced performance of the RL cycle may be interpreted as a
mitigation of quantum friction [93, 43].

Also in this setup, we verified that the discovered cycles are quite robust
across the $5$ independent training runs, displaying only minor changes (see
Fig. 9 of Methods for two cycles found in independent training runs
corresponding to Figs. 5D-E).

## Discussion

We introduced a model-free framework, based on Reinforcement Learning, to
discover Pareto-optimal thermodynamic cycles that describe the best possible
trade-off between power and efficiency of out-of-equilibrium quantum thermal
machines (heat engines and refrigerators). Our algorithm only requires
monitoring the heat fluxes of the QTM, making it a model-free approach. It can
therefore be used both for the theoretical optimization of known systems, and
potentially for the direct optimization of experimental devices for which no
model is known, and in the absence of any measurement performed on the quantum
system. Using state-of-the-art machine learning techniques, we demonstrate the
validity of our method applying it to two different prototypical setups. Our
black-box method discovered elaborate cycles that outperform previously
proposed cycles and are on par with a previous RL method that observes the
full quantum state [69]. Up to minor details, the cycles found by our method
are reproducible across independent training runs. Physically we find that
Otto cycles, commonly studied in literature, are not generally optimal, and
that optimal cycles balance a fast operation of the cycle, with the mitigation
of quantum friction.

Our method paves the way for a systematic use of RL in the field of quantum
thermodynamics. Future directions include investing larger systems to uncover
the impact of quantum many-body effects on the performance of QTMs, optimizing
systems in the presence of noise, and optimizing trade-offs that include power
fluctuations [108, 109, 99, 110].

## Methods

In this section we provide details on the optimization of the entropy
production, on the reinforcement learning implementation, on the physical
model used to describe the quantum thermal machines, on the training details,
on the convergence of the method, on the comparison with other methods, and on
the computation of the generation of coherence during the cycles. We also
provide access to the full code that was used to generate the results
presented in the manuscript, and the corresponding data.

Optimizing the entropy production

Here we discuss the relation between optimizing the power and the entropy
production, or the power and the efficiency. We start by noticing that we can
express the efficiency of a heat engine $\eta_{\text{heat}}$ and the
coefficient of performance of a refrigerator $\eta_{\text{cool}}$ in terms of
the averaged power and entropy production, i.e.

$\eta_{\nu}=\eta_{\nu}^{\text{(c)}}\,[1+\expectationvalue{\Sigma}/(\beta_{\nu}\expectationvalue*{P_{\nu}})]^{-1},$
(9)

where $\nu=\text{heat},\text{cool}$, $\eta_{\text{heat}}^{\text{(c)}}\equiv
1-\beta_{\text{H}}/\beta_{\text{C}}$ is the Carnot efficiency,
$\eta_{\text{cool}}^{\text{(c)}}\equiv\beta_{\text{H}}/(\beta_{\text{C}}-\beta_{\text{H}})$
is the Carnot coefficient of performance, and where we defined
$\beta_{\text{heat}}\equiv\beta_{\text{C}}$ and
$\beta_{\text{cool}}\equiv\beta_{\text{C}}-\beta_{\text{H}}$. We now show
that, thanks to this dependence of $\eta_{\nu}$ on
$\expectationvalue*{P_{\nu}}$ and $\expectationvalue*{\Sigma}$, if a cycle is
a Pareto-optimal trade-off between high power and high efficiency, then it is
also a Pareto-optimal trade-off between high power and low entropy-production
up to a change of $c$. This means that if we find all optimal trade-offs
between high power and low entropy-production (as we do with our method if the
Pareto-front is convex), we will have necessarily also found all Pareto-
optimal trade-offs between high power and high efficiency.

Mathematically, we want to prove that the cycles that maximize

$\expectationvalue*{G_{\nu}(c)}\equiv
c\expectationvalue*{P_{\nu}}+(1-c)\eta_{\nu}$ (10)

for some value of $c\in[0,1]$, also maximize the return in Eq. (5) for some
(possibly different) value of $c\in[0,1]$. To simplify the proof and the
notation, we consider the following two functions

$\displaystyle F(a,b,\theta)$
$\displaystyle=aP(\theta)-b\Sigma(P(\theta),\eta(\theta)),$ (11)
$\displaystyle G(a,b,\theta)$ $\displaystyle=aP(\theta)+b\eta(\theta),$

where $P(\theta)$ and $\eta(\theta)$ represent the power and efficiency of a
cycle parameterized by a set of parameters $\theta$, $a>0$ and $b>0$ are two
scalar quantities, and

$\Sigma(P,\eta)=\frac{\eta^{\text{(c)}}_{\nu}-\eta}{\eta}\beta_{\nu}P$ (12)

is obtained by inverting Eq. (9).

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

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

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

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

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

$0=\partial_{i}G(a_{1},b_{1},\theta_{1})=a_{1}\partial_{i}P(\theta_{1})+b_{1}\partial_{i}\eta(\theta_{1}).$
(15)

Now, let us compute the derivative in $\theta$ of $F(a_{2},b_{2},\theta_{1})$,
where $a_{2}>0$ and $b_{2}>0$ are two arbitrary positive coefficients. We have

$\partial_{i}F(a_{2},b_{2},\theta_{1})=(a_{2}-b_{2}\partial_{P}\Sigma)\partial_{i}P(\theta_{1})-(b_{2}\partial_{\eta}\Sigma)\partial_{i}\eta(\theta_{1}).$
(16)

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

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

thanks to Eq. (15) we have that

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

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. (17). The invertibility of Eq. (17) (i.e. a non-null determinant
of the matrix) is guaranteed by Eq. (13). We also have to make sure that if
$a_{1}>0$ and $b_{1}>0$, then also $a_{2}>0$ and $b_{2}>0$. To do this, we
invert Eq. (17), finding

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

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

To conclude the proof, we show that $\theta_{1}$ is a local maximum for
$F(a_{2},b_{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},\theta)$, we have that the Hessian matrix

$H^{(G)}_{ij}\equiv\partial_{ij}G(a_{1},b_{1},\theta_{1})=a_{1}\partial_{ij}P+b_{1}\partial_{ij}\eta$
(20)

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

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

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},$ (22)

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},\theta_{1})$ in the special point
$(a_{2},b_{2})$ previously identified, we substitute Eq. (19) into Eq. (21),
yielding

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

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. (13) and that $b_{1}>0$, we now need to show that
$Q_{ij}$ is positive semi-definite. Plugging Eq. (14) into Eq. (22) yields

$Q_{ij}=\beta_{[\nu]}\frac{\eta^{\text{(c)}}_{[\nu]}}{\eta^{2}}\partial_{i}\eta\,\partial_{j}\eta\,R_{ij},$
(24)

where

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

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 semidefinite 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_{[\nu]},\eta^{\text{(c)}}_{[\nu]},\eta>0$, we
have

$\sum_{ij}a_{i}Q_{ij}a_{j}=\beta_{[\nu]}\frac{\eta^{\text{(c)}}_{[\nu]}}{\eta^{2}}\sum_{ij}x_{i}R_{ij}x_{j}\geq
0,$ (26)

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. (25),
$2P/\eta$, 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. (15), 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}$ and $b_{2}$,
preserving the same exact optimization problem. This way, a value of
$c\in[0,1]$ can be identified.

Reinforcement Learning Implementation

As discussed in the main text, our goal is to maximize the return
$\expectationvalue{r_{c}}(t)$ defined in Eq. (5). To solve the problem within
the RL framework, we discretize time as $t_{i}=i\Delta t$. At every time-step
$t_{i}$, the aim of the agent is to learn an optimal policy that maximizes, in
expectation, the time-discretized return $\expectationvalue*{r_{c}}_{i}$. The
time-discrete reward and return functions are given by:

$\displaystyle r_{i+1}$ $\displaystyle=$ $\displaystyle\Delta
t^{-1}\int_{t_{i}}^{t_{i}+\Delta_{t}}r_{c}(t)\mathrm{d}t,$ (27)
$\displaystyle\expectationvalue*{r_{c}}_{i}$ $\displaystyle=$
$\displaystyle(1-\gamma)\sum_{j=0}^{\infty}\gamma^{j}r_{i+1+j}.$ (28)

Eq. (28) is the time-discrete version of Eq. (5), where the discount factor
$\gamma=\exp(-\kappa\Delta t)$ determines the averaging timescale and
expresses how much we are interested in future or immediate rewards.

To be precise, plugging Eq. (27) into Eq. (28) gives $\langle r_{c}\rangle(t)$
(up to an irrelevant constant prefactor) only in the limit of $\Delta t\to 0$.
However, also for finite $\Delta t$, both quantities are time-averages of the
reward, so they are equally valid definitions to describe a long-term trade-
off maximization.

As in Ref. [69], we use a generalization of the soft-actor critic (SAC)
method, first developed for continuous actions [70, 71], to handle a
combination of discrete and continuous actions [72, 73]. We further tune the
method to stabilize the convergence in a multi-objective scenario. We here
present an overview of our implementation of SAC putting special emphasis on
the differences with respect to the standard implementation. However, we refer
to [70, 71, 72, 73] for additional details. Our method, implemented with
PyTorch, is based on modifications and generalizations of the SAC
implementation provided by Spinning Up from OpenAI [111]. All code and data to
reproduce the experiments is available online (see Data Availability and Code
Availability sections).

The SAC algorithm is based on policy iteration, i.e. it consists of iterating
multiple times over two steps: a policy evaluation step, and a policy
improvement step. In the policy evaluation step, the value function of the
current policy is (partially) learned, whereas in the policy improvement step
a better policy is learned by making use of the value function. We now
describe these steps more in detail.

In typical RL problems, the optimal policy $\pi^{*}(s|a)$ is defined as the
policy that maximizes the expected return defined in Eq. (28), i.e.:

$\pi^{*}=\mathrm{arg}\max_{\pi}\,\mathop{\mathrm{E}_{\pi}}\limits_{s\sim\mu_{\pi}}\Big{[}\sum_{k=0}^{\infty}\gamma^{k}\,r_{k+1}\Big{|}s_{0}=s\Big{]},$
(29)

where $\mathrm{E}_{\pi}$ denotes the expectation value choosing actions
according to the policy $\pi$. The initial state $s_{0}=s$ is sampled from
$\mu_{\pi}$, i.e. the steady-state distribution of states that are visited by
$\pi$. In the SAC method, balance between exploration and exploitation [112]
is achieved by introducing an Entropy-Regularized maximization objective. In
this setting, the optimal policy $\pi^{*}$ is given by

$\pi^{*}=\mathrm{arg}\max_{\pi}\,\mathop{\mathrm{E}_{\pi}}\limits_{s\sim\mathcal{B}}\Big{[}\sum_{k=0}^{\infty}\gamma^{k}\,\Big{(}r_{k+1}+\alpha
H[\pi(\cdot|s_{k})]\Big{)}\Big{|}s_{0}=s\Big{]},$ (30)

where $\alpha\geq 0$ is known as the “temperature” parameter that balances the
trade-off between exploration and exploitation, and

$H[P]=\mathop{\mathrm{E}}\limits_{x\sim P}[-\log P(x)]$ (31)

is the entropy of the probability distribution $P$. Notice that we replaced
the unknown state distribution $\mu_{\pi}$ with $\mathcal{B}$, which is a
replay buffer populated during training by storing the observed one-step
transitions $(s_{k},a_{k},r_{k+1},s_{k+1})$.

Developing on Ref. [69], we generalize such approach to a combination of
discrete and continuous actions in the following way. Let us write an
arbitrary action $a$ as $a=(u,d)$, where $u$ is the continuous action and $d$
is the discrete action (for simplicity, we describe the case of a single
continuous action, though the generalization to multiple variables is
straightforward). From now on, all functions of $a$ are also to be considered
as functions of $u,d$. We decompose the joint probability distribution of the
policy as

$\pi(u,d|s)=\pi_{\mathrm{D}}(d|s)\cdot\pi_{\mathrm{C}}(u|d,s),$ (32)

where $\pi_{\mathrm{D}}(d|s)$ is the marginal probability of taking discrete
action $d$, and $\pi_{\mathrm{C}}(u|d,s)$ is the conditional probability
density of choosing action $u$, given action $d$ (D stands for “discrete”, and
C for “continuous”). Notice that this decomposition is an exact identity, thus
allowing us to describe correlations between the discrete and the continuous
action. With this decomposition, we can write the entropy of a policy as

$H[\pi(\cdot|s)]=H^{\pi}_{\text{D}}(s)+H^{\pi}_{\text{C}}(s),$ (33)

where

$\displaystyle H^{\pi}_{\text{D}}(s)$
$\displaystyle=H[\pi_{\mathrm{D}}(\cdot|s)],$ $\displaystyle
H^{\pi}_{\text{C}}(s)$
$\displaystyle=\sum_{d}\pi_{\mathrm{D}}(d|s)H[\pi_{\mathrm{C}}(\cdot|d,s)],$
(34)

correspond respectively to the entropy contribution of the discrete (D) and
continuous (C) part. These two entropies take on values in different ranges:
while the entropy of a discrete distribution with $|D|$ discrete actions is
non-negative and upper bounded by $\log(|D|)$, the (differential) entropy of a
continuous distribution can take on any value, including negative values
(especially for peaked distributions). Therefore, we introduce a separate
temperature for the discrete and continuous contributions replacing the
definition of the optimal policy in Eq. (30) with

$\pi^{*}=\mathrm{arg}\max_{\pi}\,\mathop{\mathrm{E}_{\pi}}\limits_{s\sim\mathcal{B}}\Big{[}\sum_{k=0}^{\infty}\gamma^{k}\,\Big{(}r_{k+1}+\alpha_{\text{D}}H^{\pi}_{\text{D}}(s_{k})\\\
+\alpha_{\text{C}}H^{\pi}_{\text{C}}(s_{k})\Big{)}\Big{|}s_{0}=s\Big{]},$ (35)

where $\alpha_{\text{C}}\geq 0$ and $\alpha_{\text{D}}\geq 0$ are two distinct
“temperature” parameters. This is one of the differences with respect to Refs.
[70, 71, 69]. Equation (35) defines our optimization objective. Accordingly,
we define the value function $Q^{\pi}(s,a)$ of a given policy $\pi$ as

$Q^{\pi}(s,a)=\text{E}_{\pi}\Big{[}r_{1}+\sum_{k=1}^{\infty}\gamma^{k}\,\Big{(}r_{k+1}+\alpha_{\text{D}}H^{\pi}_{\text{D}}(s_{k})\\\
+\alpha_{\text{C}}H^{\pi}_{\text{C}}(s_{k})\Big{)}\Big{|}s_{0}=s,a_{0}=a\Big{]}.$
(36)

Its recursive Bellman equation therefore reads

$Q^{\pi}(s,a)=\underset{s_{1}\atop
a_{1}\sim\pi(\cdot|s_{1})}{\text{E}}\Big{[}r_{1}+\gamma\Big{(}Q^{\pi}(s_{1},a_{1})+\alpha_{\text{D}}H^{\pi}_{\text{D}}(s_{1})\\\
+\alpha_{\text{C}}H^{\pi}_{\text{C}}(s_{1})\Big{)}\Big{|}s_{0}=s,a_{0}=a\Big{]}.$
(37)

As in Ref. [70, 71], we parameterize $\pi_{\mathrm{C}}(u|d,s)$ as a squashed
Gaussian policy, i.e. as the distribution of the variable

$\displaystyle\tilde{u}(\xi|d,s)$
$\displaystyle=u_{\text{a}}+\frac{u_{\text{b}}-u_{\text{a}}}{2}[1+\tanh\left(\mu(d,s)+\sigma(d,s)\cdot\xi)\right)],$
(38) $\displaystyle\xi$ $\displaystyle\sim\mathcal{N}(0,1),$

where $\mu(d,s)$ and $\sigma(d,s)$ represent respectively the mean and
standard deviation of the Gaussian distribution, $\mathcal{N}(0,1)$ is the
normal distribution with zero mean and unit variance, and where we assume that
$\mathcal{U}=[u_{\mathrm{a}},u_{\mathrm{b}}]$. This is the so-called
reparameterization trick.

We now describe the policy evaluation step. In the SAC algorithm, we learn two
value functions $Q_{\phi_{i}}(s,a)$ described by the learnable parameters
$\phi_{i}$, for $i=1,2$. $Q_{\phi}(s,a)$ is a function approximator, e.g. a
neural network. Since $Q_{\phi_{i}}(s,a)$ should satisfy the Bellman Eq. (37),
we define the loss function for $Q_{\phi_{i}}(s,a)$ as the mean square
difference between the left and right hand side of Eq. (37), i.e.

$L_{Q}(\phi_{i})=\mathop{\mathrm{E}}\limits_{(s,a,r,s^{\prime})\sim\mathcal{B}}\left[(Q_{\phi_{i}}(s,a)-y(r,s^{\prime}))^{2}\right],$
(39)

where

$y(r,s^{\prime})=r+\gamma\underset{a^{\prime}\sim\pi(\cdot|s^{\prime})}{\text{E}}\Big{[}\min_{j=1,2}Q_{\phi_{\text{targ},j}}(s^{\prime},a^{\prime})+\alpha_{\text{D}}H_{\text{D}}(s^{\prime})\\\
+\alpha_{\text{C}}H_{\text{C}}(s^{\prime})\Big{]}.$ (40)

Notice that in Eq. (40) we replaced $Q^{\pi}$ with
$\min_{j=1,2}Q_{\phi_{\mathrm{targ},j}}$, where $\phi_{\mathrm{targ},j}$, for
$j=1,2$, are target parameters which are not updated when minimizing the loss
function; instead, they are held fixed during backpropagation, and then they
are updated according to Polyak averaging, i.e.

$\phi_{\mathrm{targ},i}\leftarrow\rho_{\mathrm{polyak}}\phi_{\mathrm{targ},i}+(1-\rho_{\mathrm{polyak}})\phi_{i},$
(41)

where $\rho_{\mathrm{polyak}}$ is a hyperparameter. This change was shown to
improve learning [70, 71]. In order to evaluate the expectation value in Eq.
(40), we use the decomposition in Eq. (32) to write

$\mathop{\mathrm{E}}\limits_{a^{\prime}\sim\pi(\cdot|s^{\prime})}[\cdot]=\sum_{{d}^{\prime}}\pi_{\mathrm{D}}({d}^{\prime}|s^{\prime})\mathop{\mathrm{E}}\limits_{{u}^{\prime}\sim\pi_{\text{C}}(\cdot|d^{\prime},s^{\prime})}[\cdot],$
(42)

where we denote $a^{\prime}=(u^{\prime},d^{\prime})$. Plugging Eq. (42) into
Eq. (40) and writing the entropies explicitly as expectation values yields

$y(r,s^{\prime})=r+\gamma\sum_{{d}^{\prime}}\pi_{\text{D}}({d}^{\prime}|s^{\prime})\cdot\Big{(}\\\
\mathop{\mathrm{E}}\limits_{u^{\prime}\sim\pi_{\text{C}}(\cdot|{d}^{\prime},s^{\prime})}\left[\min_{j=1,2}Q_{\phi_{\text{targ},j}}(s^{\prime},d^{\prime},u^{\prime})-\alpha_{\text{C}}\log\pi_{\text{C}}(u^{\prime}|d^{\prime},s^{\prime})\right]\\\
-\alpha_{\text{D}}\log\pi_{\text{D}}(d^{\prime}|s^{\prime})\Big{)}.$ (43)

We then replace the expectation value over $u^{\prime}$ in Eq. (43) with a
single sampling $u^{\prime}\sim\pi_{\text{C}}(\cdot|{d}^{\prime},s^{\prime})$
(therefore one sampling for each discrete action) performed using Eq. (38).
This corresponds to performing a full average over the discrete action, and a
single sampling of the continuous action.

We now turn to the policy improvement step. Since we introduced two separate
temperatures, we cannot use the loss function introduced in Refs. [70, 71].
Therefore, we proceed in two steps. Let us define the following function

$Z_{\pi}(s)=-\mathop{\mathrm{E}}\limits_{{a\sim\pi(\cdot|s)}}\left[Q^{\pi^{\text{old}}}(s,a)\right]-\alpha_{\text{D}}H_{\text{D}}^{\pi}(s)-\alpha_{\text{C}}H^{\pi}_{\text{C}}(s),$
(44)

where $Q^{\pi^{\text{old}}}(s,a)$ is the value function of some given “old
policy” $\pi^{\text{old}}$, and $\pi$ is an arbitrary policy. First, we prove
that if a policy $\pi^{\text{new}}$ satisfies

$Z_{\pi^{\text{new}}}(s)\leq Z_{\pi^{\text{old}}}(s)$ (45)

for all values of $s$, then $\pi^{\text{new}}$ is a better policy than
$\pi^{\text{old}}$ as defined in Eq. (35). Next, we will use this property to
define a loss function that implements the policy improvement step. Equation
(45) implies that

$\mathop{\mathrm{E}}\limits_{{a\sim\pi^{\text{old}}(\cdot|s)}}\left[Q^{\pi^{\text{old}}}(s,a)\right]+\alpha_{\text{D}}H_{\text{D}}^{\pi^{\text{old}}}(s)+\alpha_{\text{C}}H^{\pi^{\text{old}}}_{\text{C}}(s)\leq\\\
\mathop{\mathrm{E}}\limits_{{a\sim\pi^{\text{new}}(\cdot|s)}}\left[Q^{\pi^{\text{old}}}(s,a)\right]+\alpha_{\text{D}}H_{\text{D}}^{\pi^{\text{new}}}(s)+\alpha_{\text{C}}H^{\pi^{\text{new}}}_{\text{C}}(s).$
(46)

We now use this inequality to show that $\pi^{\text{new}}$ is a better policy.
Starting from the Bellmann equation (37) for $Q^{\pi^{\text{old}}}$, we have
Eq. (47).

$Q^{\pi^{\text{old}}}(s,a)=\underset{s_{1}\atop
a_{1}\sim\pi^{\text{old}}(\cdot|s_{1})}{\text{E}}\Big{[}r_{1}+\gamma\Big{(}Q^{\pi^{\text{old}}}(s_{1},a_{1})+\alpha_{\text{D}}H^{\pi^{\text{old}}}_{\text{D}}(s_{1})+\alpha_{\text{C}}H^{\pi^{\text{old}}}_{\text{C}}(s_{1})\Big{)}\Big{|}s_{0}=s,a_{0}=a\Big{]}\leq\\\
\underset{s_{1}\atop
a_{1}\sim\pi^{\text{new}}(\cdot|s_{1})}{\text{E}}\Big{[}r_{1}+\gamma\Big{(}Q^{\pi^{\text{old}}}(s_{1},a_{1})+\alpha_{\text{D}}H^{\pi^{\text{new}}}_{\text{D}}(s_{1})+\alpha_{\text{C}}H^{\pi^{\text{new}}}_{\text{C}}(s_{1})\Big{)}\Big{|}s_{0}=s,a_{0}=a\Big{]}=\\\
\underset{s_{1}\atop
a_{1}\sim\pi^{\text{new}}(\cdot|s_{1})}{\text{E}}\Big{[}r_{1}+\gamma\Big{(}\alpha_{\text{D}}H^{\pi^{\text{new}}}_{\text{D}}(s_{1})+\alpha_{\text{C}}H^{\pi^{\text{new}}}_{\text{C}}(s_{1})\Big{)}\Big{|}s_{0}=s,a_{0}=a\Big{]}+\gamma\underset{s_{1}\atop
a_{1}\sim\pi^{\text{new}}(\cdot|s_{1})}{\text{E}}\Big{[}Q^{\pi^{\text{old}}}(s_{1},a_{1})\Big{|}s_{0}=s,a_{0}=a\Big{]}\\\
\leq\dots\leq Q^{\pi^{\text{new}}}(s,a).$ (47)

Using a strategy similar to that described in Refs. [112, 70], in Eq. (47) we
make a repeated use of inequality (46) and of the Bellmann equation for
$Q^{\pi^{\text{old}}}(s,a)$ to prove that the value function of
$\pi^{\text{new}}$ is better or equal to the value function of
$\pi^{\text{old}}$.

Let $\pi_{\theta}(a|s)$ be a parameterization of the policy function that
depends on a set of learnable parameters $\theta$. We define the following
loss function

$L_{\pi}(\theta)=\mathop{\mathrm{E}}\limits_{{s\sim\mathcal{B}\atop
a\sim\pi_{\theta}(\cdot|s)}}\left[-Q^{\pi^{\text{old}}}(s,a)-\alpha_{\text{D}}H_{\text{D}}^{\pi_{\theta}}(s)-\alpha_{\text{C}}H^{\pi_{\theta}}_{\text{C}}(s)\right].$
(48)

Thanks to Eqs. (44) and (45), this choice guarantees us to find a better
policy by minimizing $L_{\pi}(\theta)$ with respect to $\theta$. In order to
evaluate the expectation value in Eq. (48), as before we explicitly average
over the discrete action and perform a single sample of the continuous action,
and we replace $Q^{\pi^{\text{old}}}$ with $\min_{j}Q_{\phi_{j}}$. Recalling
the parameterization in Eq. (38), this yields

$\displaystyle\begin{multlined}L_{\pi}(\theta)=\underset{\begin{subarray}{c}s\sim\mathcal{B}\end{subarray}}{\text{E}}\Big{[}\sum_{d}\pi_{\text{D},\theta}(d|s)\Big{(}\alpha_{\text{D}}\log\pi_{\text{D},\theta}(d|s)+\\\
\alpha_{\text{C}}\log\pi_{{\text{C}},\theta}(\tilde{u}_{\theta}(\xi|d,s)|d,s)-\min_{j=1,2}Q_{\phi_{j}}(s,\tilde{u}_{\theta}(\xi|d,s),d)\Big{)}\Big{]},\end{multlined}L_{\pi}(\theta)=\underset{\begin{subarray}{c}s\sim\mathcal{B}\end{subarray}}{\text{E}}\Big{[}\sum_{d}\pi_{\text{D},\theta}(d|s)\Big{(}\alpha_{\text{D}}\log\pi_{\text{D},\theta}(d|s)+\\\
\alpha_{\text{C}}\log\pi_{{\text{C}},\theta}(\tilde{u}_{\theta}(\xi|d,s)|d,s)-\min_{j=1,2}Q_{\phi_{j}}(s,\tilde{u}_{\theta}(\xi|d,s),d)\Big{)}\Big{]},$
(49) $\displaystyle\quad\quad\xi\sim\mathcal{N}(0,1).$

We have defined and shown how to evaluate the loss functions $L_{Q}(\phi)$ and
$L_{\pi}(\theta)$ that allow us to determine the value function and the policy
[see Eqs. (39), (43) and (49)]. Now, we discuss how to automatically tune the
temperature hyperparameters $\alpha_{\text{D}}$ and $\alpha_{\text{C}}$. Ref.
[71] shows that constraining the average entropy of the policy to a certain
value leads to the same exact SAC algorithm with the addition of an update
rule to determine the temperatures. Let $\bar{H}_{\text{D}}$ and
$\bar{H}_{\text{C}}$ be respectively the fixed average values of the entropy
of the discrete and continuous part of the policy. We can then determine the
corresponding temperatures $\alpha_{\text{D}}$ and $\alpha_{\text{C}}$
minimizing the following two loss functions

$\displaystyle L_{\text{D}}(\alpha_{\text{D}})$
$\displaystyle=\alpha_{\text{D}}\underset{\begin{subarray}{c}s\sim\mathcal{B}\end{subarray}}{\text{E}}\left[H^{\pi}_{\text{D}}(s)-\bar{H}_{\text{D}}\right],$
(50) $\displaystyle L_{\text{C}}(\alpha_{\text{C}})$
$\displaystyle=\alpha_{\text{C}}\underset{\begin{subarray}{c}s\sim\mathcal{B}\end{subarray}}{\text{E}}\left[H^{\pi}_{\text{C}}(s)-\bar{H}_{\text{C}}\right].$

As usual, we evaluate the entropies by explicitly taking the average over the
discrete actions, and taking a single sample of the continuous action. To be
more specific, we evaluate $L_{\text{D}}$ by computing

$L_{\text{D}}(\alpha_{\text{D}})=\alpha_{\text{D}}\underset{\begin{subarray}{c}s\sim\mathcal{B}\end{subarray}}{\text{E}}\left[-\sum_{d}\pi_{\text{D}}(d|s)\log\pi_{\text{D}}(d|s)-\bar{H}_{\text{D}}\right],$
(51)

and $L_{\text{C}}$ by computing

$L_{\text{C}}(\alpha_{\text{C}})=\alpha_{\text{C}}\\\
\cdot\underset{\begin{subarray}{c}s\sim\mathcal{B}\end{subarray}}{\text{E}}\left[-\sum_{d}\pi_{\text{D}}(d|s)\underset{\begin{subarray}{c}u\sim\pi_{\text{C}}(\cdot|d,s)\end{subarray}}{\text{E}}[\log\pi_{\text{C}}(u|d,s)]-\bar{H}_{\text{C}}\right]$
(52)

and replacing the expectation value over $u$ with a single sample.

To summarize, the SAC algorithm consists of repeating over and over a policy
evaluation step, a policy improvement step, and a step where the temperatures
are updated. The policy evaluation step consists of a single optimization step
to minimize the loss functions $L_{Q}(\phi_{i})$ (for $i=1,2$), given in Eq.
(39), where $y(r,s^{\prime})$ is computed using Eq. (43). The policy
improvement step consists of a single optimization step to minimize the loss
function $L_{\pi}(\theta)$ given in Eq. (49). The temperatures are then
updated performing a single optimization step to minimize
$L_{\text{D}}(\alpha_{\text{D}})$ and $L_{\text{C}}(\alpha_{\text{C}})$ given
respectively in Eqs. (51) and (52). In all loss functions, the expectation
value over the states is approximated with a batch of experience sampled
randomly from the replay buffer $\mathcal{B}$.

We now detail how we parameterize $\pi(a|s)$ and $Q(s,a)$. The idea is to
develop an efficient way to process the state that can potentially be a long
time-series of actions. To this aim, we introduce a “convolution block” as a
building element for our NN architecture.

Figure 6: Schematic representation of the convolution block that takes as
input a 1D time-series of size $(C_{\text{in}},L_{\text{in}})$, where
$L_{\text{in}}$ is the length of the series and $C_{\text{in}}$ is the number
of channels, and produces an output of size
$(C_{\text{out}},L_{\text{in}}/2)$. In this image $L_{\text{in}}=4$. The
output is produced by stacking a 1D convolution of kernel size and stride of
$2$, and a non-linearity (left branch). A residual connection (right branch),
consisting only of linear operations, is added to improve trainability.

The convolution block, detailed in Fig. 6, takes an input of size
$(C_{\text{in}},L_{\text{in}})$, where $C_{\text{in}}$ is the number of
channels (i.e. the number of parameters determining an action at every time-
step) and $L_{\text{in}}$ is the length of the time-series, and produces an
output of size $(C_{\text{out}},L_{\text{out}}=L_{\text{in}}/2)$, thus halving
the length of the time-series. Notice that we include a skip connection (right
branch in Fig. 6) to improve trainability [102].

Using the decomposition in Eq. (32) and the parameterization in Eq. (38), the
quantities that need to be parameterized are the discrete probabilities
$\pi_{\text{D}}(d|s)$, the averages $\mu(d,s)$ and the variances
$\sigma(d,s)$, for $d=1,\dots,|D|$, $|D|=3$ being the number of discrete
actions.

Figure 7: Neural network architecture used to parameterize the policy
$\pi(\vec{u},d|s)$ (panel A) and to parameterize the value function
$Q(s,\vec{u},d)$ (panel B).

The architecture of the neural network that we use for the policy function is
shown in Fig. 7A. The state, composed of the time-series
$s_{i}=(a_{i-N},\dots,a_{i-1})$ which has shape
$(C_{\text{in}},L_{\text{in}}=N)$, is fed through a series of $\ln_{2}(N)$
convolutional blocks, which produce an output of length
$(C_{\text{out}},L=1)$. The number of input channels $C_{\text{in}}$ is
determined by stacking the components of $\vec{u}$ (which, for simplicity, is
a single real number $u$ in this appendix) and by using a one-hot encoding of
the discrete actions. We then feed this output, together with the last action
which has a privileged position, to a series of fully connected NNs with ReLU
activations. Finally, a linear network outputs $W(d|s)$, $\mu(d,s)$ and
$\log(\sigma(d,s))$, for all $d=1,\dots,|D|$. The probabilities
$\pi_{\text{D}}(d|s)$ are then produced applying the softmax operation to
$W(d|s)$.

We parameterize the value function $Q_{\phi}(s,u,d)$ as in Fig. 7B. As for the
policy function, the state $s$ is fed through $\ln_{2}(N)$ stacked convolution
blocks which reduce the length of the input to $(C_{\text{out}},L=1)$. This
output, together with the action $u$, is fed into a series of fully-connected
layers with ReLU activations. We then add a linear layer that produces $|D|$
outputs, corresponding to the value of $Q(s,u,d)$ for each $d=1,\dots,|D|$.

At last, we discuss a further change to the current method that we implemented
in the superconducting qubit refrigerator case to improve the converge. This
idea is the following. The return $\expectationvalue*{r_{c}}$ is a convex
combination of the power and of the negative entropy production. The first
term is positive when the system is delivering the desired power, while the
second term is strictly negative. Therefore, for $c$ close to $1$, the optimal
value of the return is some positive quantity. Instead, as $c$ decreases, the
optimal value of the return decreases, getting closer to zero (this can be
seen explicitly in Figs. 4A and 5B). However, a null return can also be
achieved by a trivial cycle that consists of doing nothing, i.e. of keeping
the control constant in time. Indeed, this yields both zero power, and zero
entropy production. Therefore, as $c$ decreases, it becomes harder and harder
for the RL agent to distinguish good cycles from these trivial solutions. We
thus modify our method to allow us to smoothly change the value of $c$ during
training from $1$ to the desired final value, which allows to tackle an
optimization problem by “starting from an easier problem” ($c=1$), and
gradually increasing its difficulty. This required the following modifications
to the previously described method.

We introduce two separate value functions, one for each objective (P for the
power, and $\Sigma$ for the entropy production)

$\displaystyle\begin{multlined}Q^{\pi}_{\text{P}}(s,a)=\text{E}_{\pi}\Big{[}r_{1}^{\text{(P)}}+\sum_{k=1}^{\infty}\gamma^{k}\,\Big{(}r^{(\text{P})}_{k+1}+\alpha_{\text{D}}H^{\pi}_{\text{D}}(s_{k})\\\
+\alpha_{\text{C}}H^{\pi}_{\text{C}}(s_{k})\Big{)}\Big{|}s_{0}=s,a_{0}=a\Big{]},\end{multlined}Q^{\pi}_{\text{P}}(s,a)=\text{E}_{\pi}\Big{[}r_{1}^{\text{(P)}}+\sum_{k=1}^{\infty}\gamma^{k}\,\Big{(}r^{(\text{P})}_{k+1}+\alpha_{\text{D}}H^{\pi}_{\text{D}}(s_{k})\\\
+\alpha_{\text{C}}H^{\pi}_{\text{C}}(s_{k})\Big{)}\Big{|}s_{0}=s,a_{0}=a\Big{]},$
(53)
$\displaystyle\begin{multlined}Q^{\pi}_{\Sigma}(s,a)=\text{E}_{\pi}\Big{[}r_{1}^{(\Sigma)}+\sum_{k=1}^{\infty}\gamma^{k}\,\Big{(}r^{(\Sigma)}_{k+1}+\alpha_{\text{D}}H^{\pi}_{\text{D}}(s_{k})\\\
+\alpha_{\text{C}}H^{\pi}_{\text{C}}(s_{k})\Big{)}\Big{|}s_{0}=s,a_{0}=a\Big{]},\end{multlined}Q^{\pi}_{\Sigma}(s,a)=\text{E}_{\pi}\Big{[}r_{1}^{(\Sigma)}+\sum_{k=1}^{\infty}\gamma^{k}\,\Big{(}r^{(\Sigma)}_{k+1}+\alpha_{\text{D}}H^{\pi}_{\text{D}}(s_{k})\\\
+\alpha_{\text{C}}H^{\pi}_{\text{C}}(s_{k})\Big{)}\Big{|}s_{0}=s,a_{0}=a\Big{]},$

where

$\displaystyle r_{i+1}^{\text{(P)}}$ $\displaystyle\equiv\frac{1}{\Delta
t}\int\limits_{t_{i}}^{t_{i}+\Delta t}\frac{P(\tau)}{P_{0}}\,d\tau,$
$\displaystyle r_{i+1}^{(\Sigma)}$ $\displaystyle\equiv\frac{1}{\Delta
t}\int\limits_{t_{i}}^{t_{i}+\Delta t}\frac{\Sigma(\tau)}{\Sigma_{0}}\,d\tau,$
(54)

represent respectively the normalized average power and average entropy
production during each time-step. Since the value functions in Eq. (53) are
identical to Eq. (36) up to a change of the reward, they separately satisfy
the same Bellmann equation as in Eq. (37), with $r_{1}$ replaced respectively
with $r_{1}^{\text{(P)}}$ and $r_{1}^{(\Sigma)}$. Therefore, we learn each
value functions minimizing the same loss function $L_{Q}$ given in Eq. (39),
with $r_{i}$ replaced with $r_{1}^{\text{(P)}}$ or $r_{1}^{(\Sigma)}$. Both
value functions are parameterized using the same architecture, but separate
and independent parameters. We now turn to the determination of the policy.
Comparing the definition of $r_{i}$ given in the main text with Eq. (54), we
see that $r_{i+1}=cr_{i+1}^{\text{(P)}}-(1-c)r_{i+1}^{(\Sigma)}$. Using this
property, and comparing Eq. (36) with Eq. (53), we see that

$Q^{\pi}(s,a)=cQ^{\pi}_{\text{P}}(s,a)-(1-c)Q^{\pi}_{\Sigma}(s,a).$ (55)

Therefore, we learn the policy minimizing the same loss function as in Eq.
(49), using Eq. (55) to compute the value function. To summarize, this method
allows us to vary $c$ dynamically during training. This requires learning two
value functions, one for each objective, and storing in the replay buffer the
two separate rewards $r_{i}^{\text{(P)}}$ and $r_{i}^{(\Sigma)}$.

At last, when we refer to “final deterministic cycle”, we are sampling from
the policy function “switching off the stochasticity”, i.e. choosing
continuous actions $u$ setting $\xi=0$ in Eq. (38), and choosing
deterministically the discrete action with the highest probability.

Physical model

As discussed in the main text, we describe the dynamics of the two analyzed
QTMs employing the Lindblad master equation that can be derived also for non-
adiabatic drivings [107], in the weak system-bath coupling regime performing
the usual Born-Markov and secular approximation [104, 105, 106] and neglecting
the Lamb-shift contribution. This approach describes the time-evolution of the
reduced density matrix of the quantum system, $\hat{\rho}(t)$, under the
assumption of weak system-bath interaction. Setting $\hbar=1$, the master
equation reads

$\frac{\partial}{\partial
t}{\hat{\rho}}(t)=-i\left[\hat{H}[\vec{u}(t)],\hat{\rho}(t)\right]+\sum\nolimits_{\alpha}\mathcal{D}^{(\alpha)}_{\vec{u}(t),d(t)}[\hat{\rho}(t)],$
(56)

where $\hat{H}[\vec{u}(t)]$ is the Hamiltonian of the quantum system that
depends explicitly on time via the control parameters $\vec{u}(t)$,
$[\cdot,\cdot]$ denotes the commutator, and
$\mathcal{D}^{(\alpha)}_{\vec{u}(t),d(t)}[\cdot]$, known as the dissipator,
describes the effect of the coupling between the quantum system and bath
$\alpha=\mathrm{H},\mathrm{C}$. We notice that since the RL agent produces
piece-wise constant protocols, we are not impacted by possible inaccuracies of
the master equation subject to fast parameter driving [113], provided that
$\Delta t$ is not smaller than the bath timescale. Without loss of generality,
the dissipators can be expressed as [105, 106]

$\mathcal{D}^{(\alpha)}_{\vec{u}(t),d(t)}=\lambda_{\alpha}[d(t)]\sum_{k}\gamma^{(\alpha)}_{k,\vec{u}(t)}\Big{(}\hat{A}^{(\alpha)}_{k,\vec{u}(t)}\hat{\rho}\hat{A}^{(\alpha)\dagger}_{k,\vec{u}(t)}\\\
-\frac{1}{2}\hat{A}^{(\alpha)\dagger}_{k,\vec{u}(t)}\hat{A}^{(\alpha)}_{k,\vec{u}(t)}\hat{\rho}-\frac{1}{2}\hat{\rho}\hat{A}^{(\alpha)\dagger}_{k,\vec{u}(t)}\hat{A}^{(\alpha)}_{k,\vec{u}(t)}\Big{)},$
(57)

where $\lambda_{\alpha}[d(t)]\in\\{0,1\\}$ are functions that determine which
bath is coupled the quantum system, $\hat{A}^{(\alpha)}_{k,\vec{u}(t)}$ are
the Lindblad operators, and $\gamma^{(\alpha)}_{k,\vec{u}(t)}$ are the
corresponding rates. In particular, $\lambda_{\text{H}}(\text{Hot})=1$,
$\lambda_{\text{C}}(\text{Hot})=0$, while $\lambda_{\text{H}}(\text{Cold})=0$,
$\lambda_{\text{C}}(\text{Cold})=1$, and
$\lambda_{\text{H}}(\text{None})=\lambda_{\text{C}}(\text{None})=0$. Notice
that both the Lindblad operators and the rates can depend on time through the
value of the control $\vec{u}(t)$. Their explicit form depends on the details
of the system, i.e. on the Hamiltonian describing the dynamics of the overall
system including the bath and the system-bath interaction. Below, we provide
the explicit form of $\hat{A}^{(\alpha)}_{k,\vec{u}(t)}$ and
$\gamma^{(\alpha)}_{k,\vec{u}(t)}$ used to model the two setups considered in
the manuscript. We adopt the standard approach to compute the instantaneous
power and heat currents [114]

$\displaystyle P(t)$
$\displaystyle\equiv-\mathop{\text{Tr}}\nolimits\left[\hat{\rho}(t)\,\frac{\partial}{\partial
t}\hat{H}[\vec{u}(t)]\right],$ (58) $\displaystyle J_{\alpha}(t)$
$\displaystyle\equiv\mathop{\text{Tr}}\nolimits\left[{\hat{H}}[\vec{u}(t)]\,\mathcal{D}^{(\alpha)}_{\vec{u}(t),d(t)}\right],$

that guarantees the validity of the first law of thermodynamics
$\partial{U}(t)/(\partial t)=-P(t)+\sum_{\alpha}J_{\alpha}(t)$, the internal
energy being defined as $U=\mathrm{Tr}[\hat{\rho}(t)\hat{H}[\vec{u}(t)]]$.

In the superconducting qubit refrigerator, we employ the model first put
forward in Ref. [49], and further studied in Refs. [55, 63]. In particular, we
consider the following Lindblad operators and corresponding rates (identifying
$k=\pm$):

$\displaystyle\hat{A}^{(\alpha)}_{+,u(t)}$
$\displaystyle=-i\ket*{e_{u(t)}}\bra*{g_{u(t)}},$
$\displaystyle\hat{A}^{(\alpha)}_{-,u(t)}$
$\displaystyle=+i\ket*{g_{u(t)}}\bra*{e_{u(t)}},$ (59)

where $\ket*{g_{u(t)}}$ and $\ket*{e_{u(t)}}$ are, respectively, the
instantaneous ground state and excited state of Eq. (7). The corresponding
rates are given by
$\gamma^{(\alpha)}_{\pm,u(t)}=S_{\alpha}[\pm\Delta\epsilon_{u(t)}]$, where
$\Delta\epsilon_{u(t)}$ is the instantaneous energy gap of the system, and

$S_{\alpha}(\Delta\epsilon)=\frac{g_{\alpha}}{2}\frac{1}{1+Q_{\alpha}^{2}(\Delta\epsilon/\omega_{\alpha}-\omega_{\alpha}/\Delta\epsilon)^{2}}\frac{\Delta\epsilon}{e^{\beta_{\alpha}\Delta\epsilon}-1}$
(60)

is the noise power spectrum of bath $\alpha$. Here $\omega_{\alpha}$,
$Q_{\alpha}$ and $g_{\alpha}$ are the base resonance frequency, quality factor
and coupling strength of the resonant circuit acting as bath
$\alpha=\text{H},\text{C}$ (see Refs. [49, 63] for details). As in Ref. [63],
we choose $\omega_{\text{C}}=2E_{0}\Delta$ and
$\omega_{\text{H}}=2E_{0}\sqrt{\Delta^{2}+1/4}$, such that the C (H) bath is
in resonance with the qubit when $u=0$ ($u=1/2$). The width of the resonance
is governed by $Q_{\alpha}$. The total coupling strength to bath $\alpha$,
plotted in Fig. 3F, is quantified by

$\gamma^{(\alpha)}_{u(t)}\equiv\gamma^{(\alpha)}_{+,u(t)}+\gamma^{(\alpha)}_{-,u(t)}.$
(61)

In the quantum harmonic oscillator based heat engine, following Ref. [43], we
describe the coupling to the baths through the Lindblad operators
$\hat{A}^{(\alpha)}_{+,u(t)}=\hat{a}_{u(t)}^{\dagger}$,
$\hat{A}^{(\alpha)}_{-,u(t)}=\hat{a}_{u(t)}$ and corresponding rates
$\gamma^{(\alpha)}_{+,u(t)}=\Gamma_{\alpha}\,n(\beta_{\alpha}u(t)\omega_{0})$
and
$\gamma^{(\alpha)}_{-,u(t)}=\Gamma_{\alpha}[1+n(\beta_{\alpha}u(t)\omega_{0})]$,
where we identify $k=\pm$.
$\hat{a}_{u(t)}=(1/\sqrt{2})\sqrt{m\omega_{0}u(t)}\,\hat{q}+i/\sqrt{m\omega_{0}u(t)}\,\hat{p}$
and $\hat{a}_{u(t)}^{\dagger}$ are respectively the (control dependent)
lowering and raising operators, $\Gamma_{\alpha}$ is a constant rate setting
the thermalization timescale of the system coupled to bath $\alpha$, and
$n(x)=[\exp(x)-1]^{-1}$ is the Bose-Einstein distribution.

Training details

We now provide additional practical details and the hyper parameters used to
produce the results of this manuscript.

In order to enforce sufficient exploration in the early stage of training, we
do the following. As in Ref. [111], for a fixed number of initial steps, we
choose random actions sampling them uniformly withing their range.
Furthermore, for another fixed number of initial steps, we do not update the
parameters to allow the replay buffer to have enough transitions.
$\mathcal{B}$ is a first-in-first-out buffer, of fixed dimension, from which
batches of transitions $(s_{k},a_{k},r_{k+1},s_{k+1},a_{k+1})$ are randomly
sampled to update the NN parameters. After this initial phase, we repeat a
policy evaluation, a policy improvement step and a temperature update step
$n_{\text{updates}}$ times every $n_{\text{updates}}$ steps. This way, the
overall number of updates coincides with the number of actions performed on
the QTM. The optimization steps for the value function and the policy are
performed using the ADAM optimizer with the standard values of $\beta_{1}$ and
$\beta_{2}$. The temperature parameters $\alpha_{\text{D}}$ and
$\alpha_{\text{C}}$ instead are determined using stochastic gradient descent
with learning rate $0.001$. To favor an exploratory behavior early in the
training, and at the same time to end up with a policy that is approximately
deterministic, we schedule the target entropies $\bar{H}_{\text{C}}$ and
$\bar{H}_{\text{D}}$. In particular, we vary them exponentially during each
step according to

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

where $a=\text{C},\text{D}$, $n_{\text{steps}}$ is the current step number,
and $\bar{H}_{a,\text{start}}$, $\bar{H}_{a,\text{end}}$ and
$\bar{H}_{a,\text{decay}}$ are hyperparameters. In the superconducting qubit
refrigerator case, we schedule the parameter $c$ according to a Fermi
distribution, i.e.

$c(n_{\text{step}})=c_{\text{end}}+(c_{\text{start}}-c_{\text{end}})\left[1+\exp\left(\frac{n_{\text{step}}-c_{\text{mean}}}{c_{\text{decay}}}\right)\right]^{-1}.$
(63)

In the harmonic oscillator engine case, to improve stability while training
for lower values of $c$, we do not vary $c$ during training, as we do in the
superconducting qubit refrigerator case. Instead, we discourage the agent from
never utilizing one of the two thermal baths by adding a negative reward if,
withing the last $N=128$ actions describing the state, less than $25$ describe
a coupling to either bath. In particular, if the number of actions
$N_{\alpha}$ where $d=\alpha$, with $\alpha=\text{Hot},\text{Cold}$ is less
than 25 in the state time-series, we sum to the reward the following penalty

$r_{penalty}=-1.4\frac{25-N_{\alpha}}{25}.$ (64)

This penalty has no impact on the final cycles where $N_{\alpha}$ is much
larger than $25$.

All hyperparameters used to produce the results of the superconducting qubit
refrigerator and of the harmonic oscillator heat engine are provided
respectively in Tables 1 and 2, where $c$ refers to the weight at which we are
optimizing the return.

Hyperparameter |  Qubit Refrigerator
---|---
Batch size | 512
Training steps | 500k
learning rate | 0.0003
$\mathcal{B}$ size | 280k
$\rho_{\text{polyak}}$ | 0.995
Channels per conv. block | (64, 64, 64, 128, 128, 128, 128)
Units per f.c. layer in $\pi$ | (256)
Units per f.c. layer in $Q^{\pi}$ | (256, 256)
Initial random steps | 5k
First update at step | 1k
$n_{\text{updates}}$ | 50
$\bar{H}_{\text{C},\text{start}}$ | 0
$\bar{H}_{\text{C},\text{end}}$ | -3.5
$\bar{H}_{\text{C},\text{decay}}$ | 440k
$c_{\text{start}}$ | 1
$c_{\text{end}}$ | $c$
$c_{\text{mean}}$ | 170k
$c_{\text{decay}}$ | 20k
Table 1: Hyperparameters used in numerical calculations relative to the superconducting qubit refrigerator that are not reported in the caption of Fig. 3. Hyperparameter |  Harmonic Engine
---|---
Batch size | 512
Training steps | 500k
learning rate | 0.0003
$\mathcal{B}$ size | 160k
$\rho_{\text{polyak}}$ | 0.995
Channels per conv. block |  (64, 64, 64, 128, 128, 128, 128)
Units per f.c. layer in $\pi$ | (256)
Units per f.c. layer in $Q^{\pi}$ | (256, 128)
Initial random steps | 5k
First update at step | 1k
$n_{\text{updates}}$ | 50
$\bar{H}_{\text{C},\text{start}}$ | -0.72
$\bar{H}_{\text{C},\text{end}}$ | -3.5
$\bar{H}_{\text{C},\text{decay}}$ | 144k
$\bar{H}_{\text{D},\text{start}}$ | $\ln 3$
$\bar{H}_{\text{D},\text{end}}$ | 0.01
$\bar{H}_{\text{D},\text{decay}}$ | 144k
Table 2: Hyperparameters used in numerical calculations relative to the
harmonic oscillator heat engine that are not reported in the caption of Fig.
5.

Convergence of the RL approach The training process presents some degree of
stochasticity, such as the initial random steps, the stochastic sampling of
actions from the policy function, and the random sampling of a batch of
experience from the replay buffer to compute an approximate gradient of the
loss functions. We thus need to evaluate the reliability of our approach.

As shown in the main text, specifically in Figs. 4 and 5, we ran the full
optimization $5$ times. Out of $65$ trainings in the superconducting qubit
refrigerator case, only $4$ failed, and out of the $55$ in the harmonic
oscillator engine, only $2$ failed, where by failed we mean that the final
return was negative. In such cases, we ran the training an additional time.

Figs. 4A and 5B display an error bar corresponding to the standard deviation,
at each value of $c$, computed over the $5$ repetitions. Instead, in Figs. 4B
and 5C we display one black dot for each individual training. As we can see,
the overall performance is quite stable and reliable.

Figure 8: Final deterministic cycle, identified in the superconducting qubit
refrigerator, at the fifth training. Same parameters and quantities are shown
as in Figs. 4C-F. Figure 9: Final deterministic cycle, identified in the
harmonic oscillator engine, at the fifth training. Same parameters and
quantities are shown as in Figs. 5D-E.

At last, we discuss the variability of the discovered cycles. The cycles shown
in Figs. 4C-F and 5D-E were chosen by selecting the largest return among the
$5$ repetitions. In Figs. 8 and 9 we display cycles discovered in the last of
the $5$ repetition, i.e. chosen without any post-selection. They correspond to
the same setups and parameters displayed in Figs. 4C-F and 5D-E. As we can
see, $5$ out of the $6$ displayed cycles are very similar to the ones
displayed in Figs. 4C-F and 5D-E, with a very slight variability. The only
exception is Fig. 8B, where the cycle has a visibly shorter period and
amplitude than the one shown in Fig. 4D. Despite this visible difference in
the cycle shape, the return of the cycle shown in Fig. 8B is $0.382$ compared
to $0.385$ of the cycle shown in Fig. 4B.

We therefore conclude that, up to minor changes, the cycles are generally
quite stable across multiple trainings.

Comparing with other methods

In Figs. 4 and 5 we compare the performance of our method respectively against
optimized trapezoidal cycles, and optimized Otto cycles. In both cases, we
also maximize the power using the RL method of Ref. [69]. We now detail how we
perform such comparison.

In the refrigerator based on a superconducting qubit, we consider the
trapezoidal cycle proposed in Ref. [49, 63], i.e. we fix

$u(t)=\frac{1}{4}\left(1+\frac{\tanh(a\,\cos\Omega t)}{\tanh(a)}\right)$ (65)

with $a=2$, and we optimize $\expectationvalue{r_{c}}$ with respect to
frequency $\Omega$. In the heat engine case based on a quantum harmonic
oscillator, we fix an Otto cycle as described in Ref. [43], i.e. a trapezoidal
cycle consisting of the 4 strokes shown in Figs. 5D-E as a dashed line, and we
optimize over the duration of each of the 4 strokes. In particular, we first
performed a grid search in the space of these four durations for $c=1$. After
identifying the largest power, we ran the Newton algorithm to further maximize
the return. We then ran the Newton algorithm for all other values of $c$.

The comparison with Ref. [69] was done using the source code provided in Ref.
[69], and using the same exact hyperparameters that were used in Ref. [69].

In particular, in the case of the refrigerator based on a superconducting
qubit, we re-ran the code using the hyperparameters reported in Table 1,
column “Figs. 3, 4”, of the Methods section of Ref. [69], and we trained for
the same number of steps (500k). We then evaluated its power and coefficient
of performance evaluating the deterministic policy (which typically has a
better performance). In the heat engine case based on a quantum harmonic
oscillator, we evaluated the performance of the cycle reported in Fig. 5a,c of
Ref. [69], whose training hyperparameters are reported in Table 1, column
“Fig. 5a”, of the Methods section of Ref. [69].

Generation of coherence

In order to quantify the coherence generated in the instantaneous eigenbasis
of the Hamiltonian in the refrigerator based on a superconducting qubit, we
evaluated the time average of relative entropy of coherence [115], defined as

$C(\hat{\rho}(t))=S(\hat{\rho}_{\mathrm{diag.}}(t))-S(\hat{\rho}(t)),$ (66)

where $S(\hat{\rho})=-\mathrm{Tr}[\hat{\rho}\ln\hat{\rho}]$ is the Von Neumann
entropy, and

$\hat{\rho}_{\mathrm{diag.}}(t)=\langle
g_{u(t)}|\hat{\rho}(t)|g_{u(t)}\rangle\cdot|g_{u(t)}\rangle\langle
g_{u(t)}|\\\ +\langle
e_{u(t)}|\hat{\rho}(t)|e_{u(t)}\rangle\cdot|e_{u(t)}\rangle\langle e_{u(t)}|$
(67)

is the density matrix, in the instantaneous eigenbasis $|g_{u(t)}\rangle$ and
$|e_{u(t)}\rangle$, with the off-diagonal terms canceled out.

We compute the time-average of the relative entropy of coherence generated by
the final deterministic cycle found by the RL agent, and compare it to the
coherence generated by a trapezoidal cycle operated at the same speed, i.e.
with the same period. As we can see in Table 3, the trapezoidal cycles
generate twice as much coherence as the RL cycles shown in Figs. 4C-F, i.e.
corresponding to $c=1,0.8,0.6,0.4$.

$c$ |  RL |  Trapez.
---|---|---
1 | 0.068 | 0.13
0.8 | 0.050 | 0.12
0.6 | 0.054 | 0.092
0.4 | 0.035 | 0.090
Table 3: Coherence generated by the final deterministic cycles identified by
the RL method (RL column) and generated by a trapezoidal cycle operated at the
same speed (Trapez. column) at the values of $c$ shown in the first column.
These values correspond to the cycles shown in Figs. 4C-F.

## Code and data availability

The code used to generate all results is available on GitHub
(https://github.com/PaoloAE/paper_rl_blackbox_thermal_machines). All raw data
that was generated with the accompanying code and that was used to produce the
results in the manuscript is available on Figshare
(https://doi.org/10.6084/m9.figshare.19180907).

## Acknowledgements

We are greatly thankful to Martí Perarnau-Llobet, Paolo Abiuso and Alberto
Rolandi for useful discussions and for suggesting to include the entropy
production in the return. We gratefully acknowledge 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, AA2-18).

## Competing interests

The authors declare no competing interests. P.A.E. and F.N. are authors of a
patent application containing aspects of this work (Application to the
European Patent Office, file number: 21 191 966.7).

## Author Contributions

P.A.E. and F.N. designed the research and method. P.A.E. wrote the computer
code and carried out the numerical calculations. P.A.E. and F.N. analysed the
data and wrote the manuscript.

## 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_ , Fundamental Theories of Physics (Springer International Publishing, 2019).
  * Ronzani  _et al._ [2018] A. Ronzani, B. Karimi, J. Senior, Y.-C. Chang, J. T. Peltonen, C.-D. Chen, and J. P. Pekola, Tunable photonic heat transport in a quantum heat valve, Nat. Phys. 14, 991 (2018).
  * Dutta  _et al._ [2019] B. Dutta, D. Majidi, A. G. Corral, P. A. Erdman, S. Florens, T. A. Costi, H. Courtois, and C. B. Winkelmann, Direct probe of the seebeck coefficient in a kondo-correlated single-quantum-dot transistor, Nano Lett. 19, 506 (2019).
  * Senior  _et al._ [2020] J. Senior, A. Gubaydullin, B. Karimi, J. T. Peltonen, J. Ankerhold, and J. P. Pekola, Heat rectification via a superconducting artificial atom, Commun. Phys. 3, 40 (2020).
  * Maillet  _et al._ [2020] O. Maillet, D. Subero, J. T. Peltonen, D. S. Golubev, and J. P. Pekola, Electric field control of radiative heat transfer in a superconducting circuit, Nat. Commun. 11, 4326 (2020).
  * Roßnagel  _et al._ [2016] J. Roßnagel, S. T. Dawkins, K. N. Tolazzi, O. Abah, E. Lutz, F. Schmidt-Kaler, and K. Singer, A single-atom heat engine, Science 352, 325 (2016).
  * Josefsson  _et al._ [2018] M. Josefsson, A. Svilans, A. M. Burke, E. A. Hoffmann, S. Fahlvik, C. Thelander, M. Leijnse, and H. Linke, A quantum-dot heat engine operating close to the thermodynamic efficiency limits, Nat. Nanotechnol. 13, 920 (2018).
  * Klatzow  _et al._ [2019] J. Klatzow, J. N. Becker, P. M. Ledingham, C. Weinzetl, K. T. Kaczmarek, D. J. Saunders, J. Nunn, I. A. Walmsley, R. Uzdin, and E. Poem, Experimental demonstration of quantum effects in the operation of microscopic heat engines, Phys. Rev. Lett. 122, 110601 (2019).
  * von Lindenfels  _et al._ [2019] D. von Lindenfels, O. Gräb, C. T. Schmiegelow, V. Kaushal, J. Schulz, M. T. Mitchison, J. Goold, F. Schmidt-Kaler, and U. G. Poschinger, Spin heat engine coupled to a harmonic-oscillator flywheel, Phys. Rev. Lett. 123, 080602 (2019).
  * Maslennikov  _et al._ [2019] G. Maslennikov, S. Ding, R. Hablätzel, J. Gan, A. Roulet, S. Nimmrichter, J. Dai, V. Scarani, and D. Matsukevich, Quantum absorption refrigerator with trapped ions, Nat. Commun. 10, 202 (2019).
  * Peterson  _et al._ [2019] J. P. S. Peterson, T. B. Batalhão, M. Herrera, A. M. Souza, R. S. Sarthour, I. S. Oliveira, and R. M. Serra, Experimental characterization of a spin quantum heat engine, Phys. Rev. Lett. 123, 240601 (2019).
  * Prete  _et al._ [2019] D. Prete, P. A. Erdman, V. Demontis, V. Zannier, D. Ercolani, L. Sorba, F. Beltram, F. Rossella, F. Taddei, and S. Roddaro, Thermoelectric conversion at 30 k in inas/inp nanowire quantum dots, Nano Lett. 19, 3033 (2019).
  * Horne  _et al._ [2020] N. V. Horne, D. Yum, T. Dutta, P. Hänggi, J. Gong, D. Poletti, and M. Mukherjee, Single-atom energy-conversion device with a quantum load, NPJ Quantum Inf. 6, 37 (2020).
  * Krantz  _et al._ [2019] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, A quantum engineer’s guide to superconducting qubits, Appl. Phys. Rev. 6, 021318 (2019).
  * Huang [1987] K. Huang,  _Statistical Mechanics, 2nd ed._ (Wiley, 1987).
  * 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).
  * 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).
  * 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).
  * Scandi and Perarnau-Llobet [2019] M. Scandi and M. Perarnau-Llobet, Thermodynamic length in open quantum systems, Quantum 3, 197 (2019).
  * 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).
  * 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).
  * Eglinton and Brandner [2022] J. Eglinton and K. Brandner, Geometric bounds on the power of adiabatic thermal machines, Phys. Rev. E 105, L052102 (2022).
  * Abiuso and Perarnau-Llobet [2020] P. Abiuso and M. Perarnau-Llobet, Optimal cycles for low-dissipation heat engines, Phys. Rev. Lett. 124, 110606 (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).
  * 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).
  * Arrachea  _et al._ [2007] L. Arrachea, M. Moskalets, and L. Martin-Moreno, Heat production and energy balance in nanoscale engines driven by time-dependent fields, Phys. Rev. B 75, 245420 (2007).
  * 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).
  * 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).
  * Campisi  _et al._ [2015] M. Campisi, J. Pekola, and R. Fazio, Nonequilibrium fluctuations in quantum heat engines: theory, example, and possible solid state experiments, New J. Phys. 17, 035012 (2015).
  * Dann and Kosloff [2020] R. Dann and R. Kosloff, Qunatum signatures in the quantum carnot cycle, New J. Phys. 22, 013055 (2020).
  * Molitor and Landi [2020] O. A. D. Molitor and G. T. Landi, Stroboscopic two-stroke quantum heat engines, Phys. Rev. A 102, 042217 (2020).
  * Shaghaghi  _et al._ [2022] V. Shaghaghi, G. M. Palma, and G. Benenti, Extracting work from random collisions: A model of a quantum heat engine, Phys. Rev. E 105, 034101 (2022).
  * Cavaliere  _et al._ [2022] F. Cavaliere, M. Carrega, G. De Filippis, V. Cataudella, G. Benenti, and M. Sassetti, Dynamical heat engines with non-markovian reservoirs, Phys. Rev. Res. 4, 033233 (2022).
  * Feldmann  _et al._ [1996] T. Feldmann, E. Geva, R. Kosloff, and P. Salamon, Heat engines in finite time governed by master equations, Am. J. Phys. 64, 485 (1996).
  * 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).
  * Rezek and Kosloff [2006] Y. Rezek and R. Kosloff, Irreversible performance of a quantum harmonic heat engine, New J. Phys. 8, 83 (2006).
  * Quan  _et al._ [2007] H. Quan, Y. Liu, C. Sun, and F. Nori, Quantum thermodynamic cycles and quantum heat engines, Phys. Rev. E 76, 031105 (2007).
  * Abah  _et al._ [2012] O. Abah, J. Roßnagel, G. Jacob, S. Deffner, F. Schmidt-Kaler, K. Singer, and E. Lutz, Single-ion heat engine at maximum power, Phys. Rev. Lett. 109, 203006 (2012).
  * 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).
  * Zhang  _et al._ [2014] K. Zhang, F. Bariani, and P. Meystre, Quantum optomechanical heat engine, Phys. Rev. Lett. 112, 150602 (2014).
  * Campisi and Fazio [2016] M. Campisi and R. Fazio, The power of a critical heat engine, Nat. Commun. 7, 11895 (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).
  * Kosloff and Rezek [2017] R. Kosloff and Y. Rezek, The quantum harmonic otto cycle, Entropy 19, 136 (2017).
  * Watanabe  _et al._ [2017] G. Watanabe, B. P. Venkatesh, P. Talkner, and A. del Campo, Quantum performance of thermal machines over many cycles, Phys. Rev. Lett. 118, 050601 (2017).
  * Deffner [2018] S. Deffner, Efficiency of harmonic quantum otto engines at maximal power, Entropy 20, 875 (2018).
  * Gelbwaser-Klimovsky  _et al._ [2018] D. Gelbwaser-Klimovsky, A. Bylinskii, D. Gangloff, R. Islam, A. Aspuru-Guzik, and V. Vuletic, Single-atom heat machines enabled by energy quantization, Phys. Rev. Lett. 120, 170601 (2018).
  * 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).
  * 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).
  * Das and Mukherjee [2020] A. Das and V. Mukherjee, Quantum-enhanced finite-time otto cycle, Phys. Rev. B 2, 033083 (2020).
  * Berry [2009] M. V. Berry, Transitionless quantum driving, J. Phys. A: Math. Theor. 42, 365303 (2009).
  * Deng  _et al._ [2013] Q.-h. Deng, J.and Wang, Z. Liu, P. Hänggi, and G. J., Boosting work characteristics and overall heat-engine performance via shortcuts to adiabaticity: Quantum and classical systems, Phys. Rev. E 88, 062122 (2013).
  * Torrontegui  _et al._ [2013] E. Torrontegui, S. Ibáñez, S. Martínez-Garaot, M. Modugno, A. del Campo, D. Guéry-Odelin, A. Ruschhaupt, X. Chen, and J. G. Muga, Shortcuts to adiabaticity, Adv. At., Mol., Opt. Phys. 62, 117 (2013).
  * del Campo  _et al._ [2014] A. del Campo, J. Goold, and M. Paternostro, More bang for your buck: Super-adiabatic quantum engines, Sci. Rep. 4, 6208 (2014).
  * Çakmak and Müstecaplıoğlu [2019] B. Çakmak and O. E. Müstecaplıoğlu, Spin quantum heat engines with shortcuts to adiabaticity, Phys. Rev. E 99, 032108 (2019).
  * Deng  _et al._ [2018] S. Deng, A. Chenu, P. Diao1, F. Li, S. Yu, I. Coulamy, A. del Campo, and H. Wu, Superadiabatic quantum friction suppression in finite-time thermodynamics, Sci. Adv. 18, eaar5909 (2018).
  * Funo  _et al._ [2019] K. Funo, N. Lambert, B. Karimi, J. P. Pekola, Y. Masuyama, and F. Nori, Speeding up a quantum refrigerator via counterdiabatic driving, Phys. Rev. B 100, 035407 (2019).
  * Villazon  _et al._ [2019] T. Villazon, A. Polkovnikov, and A. Chandran, Swift heat transfer by fast-forward driving in open quantum systems, Phys. Rev. A 100, 012126 (2019).
  * Khait  _et al._ [2022] I. Khait, J. Carrasquilla, and D. Segal, Optimal control of quantum thermal machines using machine learning, Phys. Rev. Res. 4, L012029 (2022).
  * 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).
  * Suri  _et al._ [2018] N. Suri, F. C. Binder, and S. Muralidharan, B. amd Vinjanampathy, Speeding up thermalisation via open quantum system variational optimisation, Eur. Phys. J. Spec. Top. 227, 203 (2018).
  * 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).
  * Erdman and Noé [2022] P. A. Erdman and F. Noé, Identifying optimal cycles in quantum thermal machines with reinforcement-learning, NPJ Quantum Inf. 8, 1 (2022).
  * Haarnoja  _et al._ [2018a] 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.
  * Haarnoja  _et al._ [2018b] 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 (2018b).
  * 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).
  * Haarnoja  _et al._ [2018c] T. Haarnoja, S. Ha, A. Zhou, J. Tan, G. Tucker, and S. Levine, Learning to walk via deep reinforcement learning, arXiv:1812.11103 (2018c).
  * 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).
  * An and Zhou [2019] Z. An and D. Zhou, Deep reinforcement learning for quantum gate control, EPL 126, 60002 (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).
  * 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._ [2022] R. Porotti, A. Essig, B. Huard, and F. Marquardt, Deep reinforcement learning for quantum state preparation with weak nonlinear measurements, Quantum 6, 747 (2022).
  * Marquardt [2021] F. Marquardt, Machine learning and quantum devices, SciPost Phys. Lect. Notes , 29 (2021).
  * Brown  _et al._ [2021] J. Brown, P. Sgroi, L. Giannelli, G. S. Paraoanu, E. Paladino, G. Falci, M. Paternostro, and A. Ferraro, Reinforcement learning-enhanced protocols for coherent population-transfer in three-level quantum systems, New J. Phys. 23, 093035 (2021).
  * Metz and Bukov [2023] F. Metz and M. Bukov, Self-correcting quantum many-body control using reinforcement learning with tensor networks, Nat. Mach. Intell. 5, 780 (2023).
  * 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).
  * 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).
  * 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).
  * 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).
  * Kosloff and Feldmann [2002] R. Kosloff and T. Feldmann, Discrete four-stroke quantum heat engine exploring the origin of friction, Phys. Rev. E 65, 055102 (2002).
  * Friedenberger and Lutz [2017] A. Friedenberger and E. Lutz, When is a quantum heat engine quantum?, EPL 120, 10002 (2017).
  * Brandner  _et al._ [2017] K. Brandner, M. Bauer, and S. U., Universal coherence-induced power losses of quantum heat engines in linear response, Phys. Rev. Lett. 119, 170602 (2017).
  * Lekscha  _et al._ [2018] J. Lekscha, H. Wilming, J. Eisert, and R. Gallego, Quantum thermodynamics with local control, Phys. Rev. E 97, 022142 (2018).
  * Strasberg  _et al._ [2016] P. Strasberg, G. Schaller, N. Lambert, and T. Brandes, Nonequilibrium thermodynamics in the strong coupling and non-markovian regime based on a reaction coordinate mapping, New J. Phys. 18, 073007 (2016).
  * Seoane and Solé [2016] L. F. Seoane and R. Solé, Multiobjective optimization 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.
  * 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).
  * Solon and Horowitz [2018] A. P. Solon and J. M. Horowitz, Phase transition in protocols minimizing work fluctuations, Phys. Rev. Lett. 120, 180605 (2018).
  * van den Oord  _et al._ [2016] A. van den Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalchbrenner, A. Senior, and K. Kavukcuoglu, Wavenet: A generative model for raw audio, arXiv:1609.03499 (2016).
  * He  _et al._ [2015] K. He, X. Zhang, S. Ren, and J. Sun, Deep residual learning for image recognition, arXiv:1512.03385 (2015).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv:1412.6980 (2014).
  * 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, 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).
  * Barato and Seifert [2015] A. C. Barato and U. Seifert, Thermodynamic uncertainty relation for biomolecular processes, Phys. Rev. Lett. 114, 158101 (2015).
  * 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. Res. 1, 033021 (2019).
  * Miller  _et al._ [2021] 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 (2021).
  * Achiam [2018] J. Achiam, Spinning Up in Deep Reinforcement Learning (2018).
  * Sutton and Barto [2018] R. S. Sutton and A. G. Barto,  _Reinforcement learning: An introduction_ (MIT press, 2018).
  * Dann  _et al._ [2018] R. Dann, A. Levy, and R. Kosloff, Time-dependent markovian quantum master equation, Phys. Rev. A 98, 052129 (2018).
  * Alicki [1979] R. Alicki, The quantum open system as a model of the heat engine, J. Phys. A: Math. Gen. 12, L103 (1979).
  * Baumgratz  _et al._ [2014] T. Baumgratz, M. Cramer, and M. Plenio, Quantifying coherence, Phys. Rev. Lett. 113, 140401 (2014).