# Optimal finite-time heat engines under constrained control

Zhuolin Ye<EMAIL_ADDRESS>Institut für Theoretische Physik,
Universität Leipzig, Postfach 100 920, D-04009 Leipzig, Germany  Federico
Cerisola  Department of Materials, University of Oxford, Parks Road, Oxford
OX1 3PH, United Kingdom  Department of Physics and Astronomy, University of
Exeter, Stocker Road, Exeter EX4 4QL, UK  Paolo Abiuso  ICFO—Institut de
Ciències Fotòniques, The Barcelona Institute of Science and Technology, 08860
Castelldefels (Barcelona), Spain  Department of Applied Physics, University of
Geneva, 1211 Geneva, Switzerland  Janet Anders  Department of Physics and
Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, UK  Institut
für Physik und Astronomie, University of Potsdam, 14476 Potsdam, Germany
Martí Perarnau-Llobet  Department of Applied Physics, University of Geneva,
1211 Geneva, Switzerland  Viktor Holubec<EMAIL_ADDRESS>Department of Macromolecular Physics, Faculty of Mathematics and Physics,
Charles University, V Holešovičkách 2, CZ-180 00 Praha, Czech Republic

###### Abstract

We optimize finite-time stochastic heat engines with a periodically scaled
Hamiltonian under experimentally motivated constraints on the bath temperature
$T$ and the scaling parameter $\lambda$. We present a general geometric proof
that maximum-efficiency protocols for $T$ and $\lambda$ are piecewise
constant, alternating between the maximum and minimum allowed values. When
$\lambda$ is restricted to a small range and the system is close to
equilibrium at the ends of the isotherms, a similar argument shows that this
protocol also maximizes output power. These results are valid for arbitrary
dynamics. We illustrate them for an overdamped Brownian heat engine, which can
experimentally be realized using optical tweezers with stiffness $\lambda$.

##  I Introduction

The unprecedented improvement in experimental control over microscopic
Brownian [1] and quantum systems [2, 3, 4] has induced a revolution in the
study of heat engines [5, 6]. It aims to generalize equilibrium and finite-
time thermodynamics [7, 8, 9, 10, 11, 12, 13, 14, 15] to the nanoscale, where
thermal and quantum fluctuations render thermodynamic variables such as work
and heat stochastic [16]. Intense effort is devoted to uncover optimal
performance of stochastic heat engines [17, 18, 19, 20, 21, 22, 23, 24, 25,
26, 27, 28, 29, 30, 31, 32, 16, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42].
However, optimal control protocols are only known under approximations of fast
[34, 35, 36] or slow [37, 38, 39, 40, 41, 42] driving, or for specific
microscopic models: engines based on overdamped Brownian particles in harmonic
[24], log-harmonic [43], or slowly-varying potentials [44], and for an
underdamped harmonic Brownian heat engine [45]. Furthermore, most these exact
results have been obtained under constraints on the state of the working
medium of the engine [46], instead of experimentally motivated constraints on
the control parameters [47, 48]. An exception is Ref. [49], showing that
reaching maximum efficiency of slowly driven cyclic heat engines requires
control over the scaling of the full Hamiltonian to avoid heat leakages.

In this letter, we optimize finite-time thermodynamic cycles under constraints
on control parameters such as trap stiffness of optical tweezers $\lambda$ and
bath temperature $T$. We show that, different from constraining the response
such as the width $\sigma$ of the phase distribution, constraining the control
allows for surprisingly simple and general derivation of maximum-efficiency
and maximum-power protocols. Besides other stark differences, for constrained
control of Brownian heat engines, these protocols significantly outperform the
protocol optimised for power and efficiency under constraints on $\sigma$
[24].

##  II Setup

Following Ref. [49], we assume that the Hamiltonian of the system that serves
as a working medium of the stochastic heat engine is of the form,

$H(x,t)=\lambda(t)f(x),$ (1)

where the control parameter $\lambda(t)$ periodically expands and shrinks the
energy spectrum in time, and $f(x)$ is an arbitrary function of the system
degrees of freedom $x$ such that the equilibrium partition function $Z(t)=\int
dx\exp[-H(x,t)/(k_{\rm B}T)]$ is finite for all $k_{\rm B}T\geq 0$. This class
of Hamiltonians generalizes the well-known ‘breathing’ parabola model [24] for
an overdamped particle trapped in a parametrically driven harmonic potential.
It also includes semi-classical two (or multi) level systems with controlled
gaps between the individual energy levels [16], and quantum spins, where the
control parameter is an externally controlled magnetic field [17].

We connect the system to a heat bath and periodically alter its temperature
$T(t)$ with the same finite period $t_{\rm p}$ as $\lambda(t)$. The parameters
under experimental control are thus $\lambda(t)$ and $T(t)$, and we assume
that experimental conditions restrict their minimum and maximum values [50].
Our aim is to find optimal $t_{\rm p}$-periodic protocols for the constrained
control parameters

$\lambda(t)\in[\lambda_{-},\lambda_{+}],\quad T(t)\in[T_{-},T_{+}].$ (2)

Figure 1: The maximum-efficiency protocol (3) under the constraints in Eq. (2)
(dashed line) compared to a suboptimal cycle (dotted line). Heat flows into
the system during the red branches of the cycles ($d\sigma>0$) and is released
to a heat sink during the blue branches ($d\sigma<0$).

##  III Maximum efficiency

Our first main result is a general geometric proof that the maximum-efficiency
finite-time cycle under the constraints (2) is a Carnot-Otto cycle composed of
two isotherms/isochores interconnected by two adiabats. The maximum-efficiency
protocol $\\{T(t),\lambda(t)\\}$ is thus piecewise constant:

$\\{T(t),\lambda(t)\\}_{\eta}=\begin{cases}\\{T_{+},\lambda_{+}\\},&0<t<t_{+}\\\
\\{T_{-},\lambda_{-}\\},&t_{+}<t<t_{\rm p}\end{cases}.$ (3)

And the maximum efficiency is given by

$\eta=1-\lambda_{-}/\lambda_{+}.$ (4)

The proof relies just on the definition of heat and it is thus independent of
the details of the system dynamics, including the times $t_{+}$ and $t_{\rm
p}$. It holds both for situations when the heat bath is memoryless (Markovian)
and non-Markovian. The non-equilibrium dynamics of the system communicating
with a Markovian bath can be described by Fokker-Planck or master equations
for the probability density for $x$ [51]. Except for a few exactly solvable
settings [16, 51], these equations are usually hard to solve analytically for
non-quasi-static time-dependent protocols. However, in the non-Markovian case,
a corresponding closed deterministic description might not be available at all
[52]. Then one has to resort to stochastic descriptions, such as a generalized
Langevin equation, making even a numerical optimization challenging. The
derivation also holds in situations with a non-equilibrium bath, such as in
recently intensely studied cyclic active Brownian heat engines [53, 54, 55,
56].

Let us now derive Eqs. (3)-(4). Under reasonable assumptions, any periodic
variation of the control parameters eventually induces a periodic average
response of the system, $\sigma(t)=\langle f[x(t)]\rangle$. This ensemble
average is a functional of $T(t)$ and $\lambda(t)$ specified by dynamical
equations of the system. Due to the factorized structure of the Hamiltonian
(1), the average internal energy of the system $\langle H(x,t)\rangle$ is
given by $\lambda(t)\sigma(t)$. Decomposing its infinitesimal change into a
component corresponding to the external variation of the control $\lambda$
(work) and the rest (heat) [6, 5], it follows that output work and input heat
increments are given by đ$W_{\rm out}(t)=-\sigma(t)d\lambda(t)$ and
đ$Q(t)=\lambda(t)d\sigma(t)$, respectively. Per cycle, the engine transforms
the fraction $\eta=W_{\rm
out}/Q_{\mathrm{in}}=1-Q_{\mathrm{out}}/Q_{\mathrm{in}}$ of the heat $Q_{\rm
in}=\int_{0}^{t_{\rm p}}\lambda(t)\theta[d\sigma(t)]d\sigma(t)$ from the heat
source into output work $W_{\rm out}=-\int_{0}^{t_{\rm
p}}\sigma(t)d\lambda(t)$, and dumps the remaining heat $Q_{\rm out}=Q_{\rm
in}-W_{\rm out}=\int_{0}^{t_{\rm p}}\lambda(t)\theta[-d\sigma(t)]d\sigma(t)$
into the heat sink ($\theta(\bullet)$ denotes the Heaviside step function
reflecting that heat flows on average into the system when $\sigma$ increases,
i.e., $d\sigma>0$).

Consider now the $\lambda-\sigma$ diagram of the cycle depicted in Fig. 1. We
seek the shape of the cycle which yields maximum efficiency $\eta$ under the
constraints (2) on the boundary values of the control parameters $\lambda$ and
$T$ 111 A similar optimisation problem is often solved in courses on classical
thermodynamics to show that maximum efficiency of an equilibrium cycle under
the constrains $T(t)\in[T_{-},T_{+}]$ on the bath temperature is the Carnot
efficiency. However, in our case, the system can be arbitrarily far form
equilibrium.. The cycle must run clockwise to secure that
$Q_{\mathrm{in}}>Q_{\mathrm{out}}$. Next, we note that maximizing $\eta$
amounts to minimizing the ratio $Q_{\mathrm{out}}/Q_{\mathrm{in}}$. For given
boundary values $\sigma_{\pm}$ of $\sigma$, this is obviously achieved by
setting $\lambda=\lambda_{+}$ when $d\sigma>0$ and $\lambda=\lambda_{-}$ when
$d\sigma<0$. In such a case, $Q_{\mathrm{in}}=\lambda_{+}\Delta\sigma$,
$Q_{\mathrm{out}}=\lambda_{-}\Delta\sigma$, and the efficiency is given by Eq.
(4). The increase in the system response $\Delta\sigma=\sigma_{+}-\sigma_{-}$,
which can be a complicated functional of the protocol $\\{T(t),\lambda(t)\\}$,
cancelled out. Eq. (4) is thus valid for arbitrary $\sigma_{\pm}$, and it
represents the maximum efficiency of a heat engine based on Hamiltonian (1)
under the constraints (2). The corresponding maximum-efficiency protocol for
$\lambda$ forms a rectangle ranging from $\lambda_{-}$ to $\lambda_{+}$ in the
$\lambda-\sigma$ diagram regardless the cycle duration and dynamical equations
of the system. The only constraint on these control parameters is that the
cycle runs in the $\lambda-\sigma$ diagram clockwise.

When not driven, a system out of equilibrium relaxes towards the equilibrium
state corresponding to the instantaneous values of the fixed control
parameters. When the control parameters change periodically, the resulting
non-equilibrium state of the system can no longer relax to equilibrium but it
lags behind the quasi-static cycle specified by the instantaneous values of
the control parameters. In our setting, $\sigma(t)$ lags behind $\sigma^{\rm
eq}(t)=\int dxf(x)\exp\\{-\lambda(t)f(x)/[k_{\rm B}T(t)]\\}/Z(t)$. In Sec. B1
of the Appendix, we show that $\sigma^{\rm eq}(t)$ is a monotonously
increasing function of $T/\lambda$. Denoting as $t_{+}$ the duration of the
$\lambda=\lambda_{+}$ branch, clockwise cycles with $\Delta\sigma>0$ are thus
obtained for temperature protocols $T(t)$ which obey (i) $\dot{T}(t)\geq 0$
when $\lambda=\lambda_{+}$, (ii) $\dot{T}(t)\leq 0$ when
$\lambda=\lambda_{-}$, and (iii) $T(t_{+}-)/\lambda_{+}>T(t_{\rm
p}-)/\lambda_{-}$. The last condition implies that the maximum efficiency (4)
obeys the standard 2nd law inequality $\eta\leq 1-T(t_{\rm p}-)/T(t_{+}-)\leq
1-T_{-}/T_{+}$. It saturates for the ‘compression ratio’
$\lambda_{-}/\lambda_{+}=T_{-}/T_{+}$. Even for a finite cycle time $t_{\rm
p}$, output power in this case vanishes. This is because $\sigma^{\rm eq}(t)$
becomes constant, yielding an infinitesimal quasi-static cycle, producing a
vanishing output work. In the maximum-efficiency protocol (3), we use the
specific protocol for $T(t)$ that maximizes the upper bound on $\eta$. In Sec.
B1 of the Appendix, we argue that this temperature protocol also maximizes the
output work of the engine regardless $\lambda(t)$ because it yields the
largest temperate differences between the bath and the system when they
exchange heat. However, we reiterate that the maximum efficiency (4) can be
achieved for an arbitrary protocol for $T(t)$ that obeys the above conditions
(i)-(iii). This freedom in $T(t)$ can be exploited in setups where precise
control of the bath (effective) temperature is difficult, such as in active
Brownian heat engines [54].

The adiabatic branches, where the protocol $\\{T(t),\lambda(t)\\}$ changes
between the boundary values, can be realised using several qualitatively
different approaches [31]. First, one can disconnect the system from the heat
bath, which might be impractical for microscopic engines. Second, one can keep
the system in thermal contact with the bath and vary the control parameters
$T$ and $\lambda$ in such a way that the response $\sigma$ does not change
[58]. The advantage of this approach is that it allows one to circumvent some
of the shortcomings of overdamped thermodynamics [59], where the heat fluxes
through the momentum degrees of freedom are neglected. Finally, one may
realise the adiabatic branches much faster than the relaxation time of the
parameter $\sigma$ [60]. In the specific maximum efficiency protocol (3), we
employ the last possibility, because it minimizes the cycle time $t_{\rm p}$
and thus yields the largest output power $P\equiv W_{\rm out}/t_{\rm p}$.
Besides, it allows for a direct comparison with the maximum-efficiency
protocols derived for Brownian heat engines under constraints on $\sigma$
[24]. However, other realisations of the adiabatic branches yield the same
maximum efficiency (4). Finally, we reiterate that also the choice of the
durations $t_{+}$ and $t_{\rm p}-t_{+}$ of the isotherms in the protocol (3)
does not affect the maximum $\eta$.

##  IV Maximum output power

If the durations of the isotherms are long enough compared to the relaxation
time of the system, i.e., $\Delta\sigma$ is close to its equilibrium value,
and the compression ratio $\lambda_{-}/\lambda_{+}$ is large, the maximum-
efficiency protocol (3) also yields maximum output work $W_{\rm{out}}$ and
power $P$ under the constrained control (2). This is our second main result.
To prove it, consider the generally unreachable geometric loose upper bound on
the output work $\max W_{\rm out}=\Delta\lambda\max\Delta\sigma^{\rm
eq}=(\lambda_{+}-\lambda_{-})[\sigma^{\rm eq}(T_{+}/\lambda_{-})-\sigma^{\rm
eq}(T_{-}/\lambda_{+})]$, which follows from the broadly valid assumption
$\max\Delta\sigma<\max\Delta\sigma^{\rm eq}$ and the insight that $W_{\rm
out}$ is given by the area enclosed by the cycle in the $\lambda-\sigma$
diagram. Expanding $\max W_{\rm out}$ in $\Delta\lambda$ yields $\max W_{\rm
out}=\Delta\lambda[\sigma^{\rm eq}(T_{+}/\lambda_{+})-\sigma^{\rm
eq}(T_{-}/\lambda_{-})]+\mathcal{O}(\Delta\lambda^{2})$. Up to the leading
order in $\Delta\lambda$ and under the condition that the system has relaxed
at the ends of the two isotherms to equilibrium, this upper bound is saturated
by the protocol (3), which completes the proof. We note that: (i) The
condition $\Delta\sigma=\Delta\sigma^{\rm eq}$ does not mean that the cycle is
slow as the system has to be close to equilibrium at the ends of the two
isotherms only and can be arbitrarily far from equilibrium otherwise. (ii)
This condition allows one to analytically calculate the whole probability
distribution for the output work regardless of additional details of the
system dynamics [61, 16]. Interestingly, for semi-classical systems, piece-
wise constant protocols with two or more branches also maximize output power
in the opposite limiting regime, when the cycle time is much shorter than the
system relaxation time [62, 35, 34].

| $\lambda_{\rm pwc}(t)$ | $\lambda_{\mathrm{pwl}}(t)$ | $\lambda_{\mathrm{slow}}(t)$ | $\lambda_{\mathrm{S}}(t)$
---|---|---|---|---
$t<t_{+}$ | $a$ | $a+bt$ | $\frac{a}{(1+bt)^{2}}$ | $\frac{T_{+}}{2\sigma_{\rm S}(a,b,t)}-\frac{\sqrt{b}-\sqrt{a}}{\mu t_{+}\sqrt{\sigma_{\rm S}(a,b,t)}}$
$t>t_{+}$ | $c$ | $c+dt$ | $\frac{c}{(1+dt)^{2}}$ | $\frac{T_{-}}{2\sigma_{\rm S}(a,b,t)}+\frac{\sqrt{b}-\sqrt{a}}{\mu t_{-}\sqrt{\sigma_{\rm S}(a,b,t)}}$
Table 1: Considered classes of protocols with free parameters $\\{a,b,c,d\\}$
determined by the target of the optimization. The resulting protocols are in
general discontinuous at times $t_{+}$ and $t_{\rm p}$. The protocols
$\\{\lambda_{\rm pwc},\lambda_{\rm S}\\}$ have two and $\\{\lambda_{\rm
pwl},\lambda_{\rm slow}\\}$ have four free parameters. The piecewise constant
protocol $\lambda_{\rm pwc}(t)$ is a variant of the maximum-efficiency
protocol (3), where $\lambda(t)$ do not have to reach the boundary values
$\lambda_{-}$ and $\lambda_{+}$. The piecewise linear protocol
$\lambda_{\mathrm{pwl}}(t)$ is perhaps the simplest one that can be
implemented in the lab. The protocol $\lambda_{\mathrm{slow}}(t)$ minimizes
the irreversible looses during isothermal branches under close-to-equilibrium
conditions. Such protocols can be derived for Brownian heat engines with
Hamiltonians of the form $\lambda(t)x^{n}/n$ (for details, see Sec. D in the
Appendix). The protocol $\lambda_{\rm S}(t)$ maximizes both power and
efficiency under the constraint that $\sigma(0)=\sigma(t_{\rm p})\equiv a$ and
$\sigma(t_{+})\equiv b$ [24]. The corresponding response $\sigma_{\rm S}(t)$
is given by Eq. (C4) in the Appendix. Note that the protocols $\\{\lambda_{\rm
pwl},\lambda_{\rm slow}\\}$ reduce to $\lambda_{\mathrm{pwc}}(t)$ for $b=d=0$.

Beyond these regimes, $W_{\rm{out}}$ and $P$ strongly depend on all details of
the dynamics through $\sigma(t)$ and cycle time $t_{\rm p}$. While
$W_{\rm{out}}$ and $P$ are still optimised by the temperature protocol and the
choice of fast adiabats in (3), optimal protocols for $\lambda(t)$ under the
constraints (2) are no longer piecewise constant and they have to be
identified for each system separately. Similarly as the derivation of maximum-
efficiency and maximum-power protocols under constraints on the system state
[24, 43, 45, 46], this often involves functional optimisation or extensive
numerical work which are both nontrivial tasks.

Figure 2: Numerical optimization of the output power of the Brownian heat
engine illustrates that the maximum-efficiency protocol (3) also yields
maximum power when the compression ratio $\lambda_{-}/\lambda_{+}$ is large
and the durations $t_{+}=t_{-}=1$ of the two isotherms are comparable to the
relaxation time $1/(2\mu\lambda_{\pm})$ for $\sigma$. (a) Powers (in units of
the ultimate maximum power $P^{*}$ for $\lambda_{\rm pwl}$) and (b) the
corresponding efficiencies obtained using the maximum-efficiency protocol
$\lambda_{\eta}$ (3) and the four classes of protocols in Tab. 1 optimised for
power. For $\lambda_{-}/\lambda_{+}\geq 0.59$ all protocols except for
$\lambda_{\rm S}$ coincide. (c) and (d) show the protocols and the resulting
response for $\lambda_{-}/\lambda_{+}=0.4$. (e) the relative differences
$\delta X=(X_{\lambda_{\rm pwl}}-X_{\lambda_{\rm pwc}})/X_{\lambda_{\rm pwl}}$
of power ($X=P$) and efficiency ($X=\eta$) for $\lambda_{\rm pwl}$ and
$\lambda_{\rm pwc}$. (f) shows the optimal values of paramters $b$ and $d$ for
$\lambda_{\rm pwl}$. Parameters used: $k_{\rm B}T_{+}$=1, $k_{\rm
B}T_{-}$=0.25, $t_{+}=t_{-}=1$, $\lambda_{+}=0.5$, $\mu=1$.

To illustrate the main features of maximum-power protocols for constrained
control, we now consider a specific Brownian heat engine based on an
overdamped particle diffusing in a controlled harmonic potential. This model
describes experimental realizations of microscopic heat engines using optical
tweezers [60, 63, 64]. Besides, the corresponding maximum-efficiency and
maximum-power protocols under the constrained response is known [24], allowing
for a direct comparison with our results. The Hamiltonian (1) now reads
$H(x,t)=\lambda(t)x^{2}/2$, with $x$ the position of the particle. The
response of the system $\sigma(t)=\langle x^{2}/2\rangle$ is proportional to
the position variance and it obeys the first order differential equation

$d\sigma(t)/dt=-2\mu\lambda(t)\sigma(t)+\mu k_{\rm B}T(t).$ (5)

It can be solved analytically for an arbitrary mobility $\mu$ of the particle
and any protocol for $\\{T(t),\lambda(t)\\}$ [16]. In what follows, we measure
time, length and energy in units of $(2\mu\lambda_{+})^{-1}$, $\sqrt{k_{\rm
B}T_{+}/(2\lambda_{+})}$, and $k_{\rm B}T_{+}$, respectively.

We ask which protocol for $\lambda$ yields the largest output work and power
under the constraints (2). Even though the model (5) is exactly solvable, the
optimal $\lambda(t)$ has to be found numerically, e.g., by the method
described in Ref. [65]. To keep the optimization transparent, we instead
consider a specific set of families of protocols for the isothermal strokes
with free parameters and then numerically optimize over them. When such
classes are chosen suitably, the resulting suboptimal performance will be
close to the global optimum [66, 67]. We consider the protocols in Tab. 1.
Besides using the protocol for temperature and adiabatic branches form (3), we
also fix during the optimisation durations of the two isotherms and thus the
cycle time. They can be further optimised once the optimal variation of
$\lambda$ is known.

Main results of the optimisation procedure are summarized in Fig. 2. (i) With
increasing minimum compression ratio $\lambda_{-}/\lambda_{+}$ allowed by the
constraints (2), maximum power for all considered protocols in (a) is first
constant and then, at an optimal compression ratio $r^{*}$, decreases. The
decreasing part corresponds to protocols which span between the allowed
boundary values, i.e., $\max\lambda(t)=\lambda(0+)=\lambda_{+}$ and
$\min\lambda(t)=\lambda(t_{+}+)=\lambda_{-}$. On the other hand, at the
plateau, the boundary values of the protocol are chosen within the bounds (2)
to keep the optimal compression ratio $r^{*}$. (ii) Values of maximum power
obtained for the protocols which have enough free parameters are
indistinguishable within our numerical precision. As the corresponding
optimized protocols seem to have minimum possible curvature
$\ddot{\lambda}(t)$, we conclude that the maximum-power protocol is piecewise
linear. (iii) Only the protocol $\lambda_{\rm S}$, optimised for constrained
response $\sigma$, yields notably smaller power than other protocols. (iv) In
agreement with our above discussion, for large enough values of
$\lambda_{-}/\lambda_{+}\geq 0.59$, the optimized parameters for the protocols
$\lambda_{\rm pwl}$, $\lambda_{\rm pwc}$, and $\lambda_{\rm slow}$ are
$b=d=0$, $a=\lambda_{+}$, and $c=\lambda_{-}$, reducing them to the maximum-
efficiency protocol (3). (v) The maximum power of the protocols
$\lambda_{\eta}$ and $\lambda_{\rm pwl}$ differ just by 1%. In Fig. D of the
Appendix, we further show that for branches durations that further optimize
output power, the relative difference is always below 12%. This suggests that
the maximum-efficiency protocol (3) can be used to obtain maximum power
approximately in situations when $\lambda_{-}$ is a free parameter over which
one can further optimise $W_{\rm out}$ and $P$.

Such optimisation seems a natural task for various experimental platforms
where the main limitation is on the maximum strength of the potential. The
maximum power regime of the maximum-efficiency protocol can be, to a large
extend, investigated analytically. First, assuming again that durations of the
isotherms are long enough that the system is close to equilibrium at times
$t_{+}$ and $t_{\rm p}$, we have
$W_{\rm{out}}=(\lambda_{+}-\lambda_{-})[\sigma^{\rm
eq}(T_{+}/\lambda_{+})-\sigma^{\rm eq}(T_{-}/\lambda_{-})]$. For
$f(x)=\lvert{x}\rvert^{n}$ in Eq. (1), we then find that the optimal
compression ratio is $\lambda_{-}/\lambda_{+}=\sqrt{T_{-}/T_{+}}$, which leads
to the output work
$W_{\mathrm{out}}=k_{\mathrm{B}}T_{+}(2\eta_{\mathrm{CA}}-\eta_{\mathrm{C}})/n$
and Curzon-Ahlborn efficiency $\eta=\eta_{\mathrm{CA}}=1-\sqrt{T_{-}/T_{+}}$
(see Sec. B1 of the Appendix for details). For other than power-law
Hamiltonians, the efficiency at maximum power can differ from
$\eta_{\mathrm{CA}}$ but it can still be determined numerically regardless
details of dynamical equations for the system (for details, see Sec. B1 in the
Appendix). Relaxing the assumption of slow (but not quasi-static) isotherms,
the optimization of $W_{\mathrm{out}}$ with respect to $\lambda_{-}$ requires
specification of the dynamics. In Fig. 3, we show that the efficiency at
maximum power of the Brownian heat engine described by Eq. (5) and driven by
the maximum-efficiency protocol (3) is bounded between the Curzon-Ahlborn
efficiency, achieved for slow isotherms, and the efficiency
$2-\sqrt{4-2\eta_{\mathrm{C}}}<\eta_{\mathrm{CA}}$, reached in the limit of
infinitely fast cycles (see Sec. B in the Appendix for details).

Figure 3: Efficiency at maximum output power ($\eta$, red axis) and the
corresponding optimal compression ratio ($\lambda^{*}_{-}/\lambda_{+}$, blue
axis) for the maximum-efficiency protocol (3) as functions of the temperature
ratio $T_{-}/T_{+}$ for six values of cycle duration $t_{\rm p}$ (coloured
curves) and $t_{-}=t_{+}$. The cycle time, $t_{\rm p}$, is measured in units
of the system relaxation time during the hot isotherm. The corresponding ratio
during the cold isotherm reads $t_{\rm p}\,\lambda^{*}_{-}/\lambda_{+}$. We
used the same parameters as in Fig. 2.

##  V Conclusion

The main strength of the presented derivations of the maximum-efficiency and
maximum-power protocols under constrained control is their simplicity and
unprecedented generality. A possible extension of these results to more
complicated Hamiltonians is sketched in Sec. A of the Appendix. However, more
general extensions remain to be explored in future work. Nevertheless, the
validity of our results for Brownian heat engines is already of experimental
relevance as they are often realized using optical tweezers with strict bounds
on the trap stiffness $\lambda$ – too small $\lambda$ leads to losing the
Brownian particle while too large $\lambda$ can induce its overheating and
even evaporation. At the same time, the achievable trap stiffnesses are well
above $10^{-6}$ N/m [31]. For a spherical Brownian particle with the radius of
$10^{-6}$ m in water, the Stokes law predicts the mobility of $\mu\approx
0.5\times 10^{8}$ m/Ns, leading to the relaxation time $1/(2\mu\lambda)$ of
the response $\sigma$ on the order of $10^{-2}$ s. The assumption that the
durations of the isotherms are comparable to the response relaxation time,
used in our derivation of the maximum-power protocol, is thus, in this setup,
natural. Besides, we believe that extensions of our results can find
applications in more involved optimization tasks, e.g., performed using
machine learning algorithms [68, 69] or geometric methods [70, 71], as well as
in quantum setups [40, 72, 73].

Acknowledgements. ZY is grateful for the sponsorship of China Scholarship
Council (CSC) under Grant No. 201906310136. FC gratefully acknowledges funding
from the Fundational Questions Institute Fund (FQXi-IAF19-01). VH gratefully
acknowledges support by the Humboldt foundation and by the Czech Science
Foundation (project No. 20-02955J). PA is supported by “la Caixa” Foundation
(ID 100010434, Grant No. LCF/BQ/DI19/11730023), and by the Government of Spain
(FIS2020-TRANQI and Severo Ochoa CEX2019-000910-S), Fundacio Cellex, Fundacio
Mir-Puig, Generalitat de Catalunya (CERCA, AGAUR SGR 1381. MPL acknowledges
financial support from the Swiss National Science Foundations (Ambizione grant
PZ00P2-186067). JA acknowledges funding from the Engineering and Physical
Sciences Research Council (EPSRC) (EP/R045577/1) and thanks the Royal Society
for support.

## References

  * Pesce  _et al._ [2020] G. Pesce, P. H. Jones, O. M. Maragò, and G. Volpe, Optical tweezers: theory and practice, Eur. Phys. J. Plus 135, 949 (2020).
  * Wineland [2013] D. J. Wineland, Nobel lecture: Superposition, entanglement, and raising schrödinger’s cat, Rev. Mod. Phys. 85, 1103 (2013).
  * Haroche [2013] S. Haroche, Nobel lecture: Controlling photons in a box and exploring the quantum to classical boundary, Rev. Mod. Phys. 85, 1083 (2013).
  * Myers  _et al._ [2022] N. M. Myers, O. Abah, and S. Deffner, Quantum thermodynamic devices: from theoretical proposals to experimental reality, AVS Quantum Sci. 4, 027101 (2022).
  * Sekimoto [2010] K. Sekimoto,  _Stochastic energetics_ , Vol. 799 (Springer, Berlin, 2010).
  * Seifert [2012] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Rep. Prog. Phys. 75, 126001 (2012).
  * Curzon and Ahlborn [1975] F. L. Curzon and B. Ahlborn, Efficiency of a carnot engine at maximum power output, Am. J. Phys. 43, 22 (1975).
  * Rubin [1979] M. H. Rubin, Optimal configuration of a class of irreversible heat engines. i, Phys. Rev. A 19, 1272 (1979).
  * Salamon  _et al._ [1980] P. Salamon, A. Nitzan, B. Andresen, and R. S. Berry, Minimum entropy production and the optimization of heat engines, Phys. Rev. A 21, 2115 (1980).
  * Andresen  _et al._ [1984a] B. Andresen, R. S. Berry, M. J. Ondrechen, and P. Salamon, Thermodynamics for processes in finite time, Acc. Chem. Res. 17, 266 (1984a).
  * Andresen  _et al._ [1984b] B. Andresen, P. Salamon, and R. S. Berry, Thermodynamics in finite time, Phys. Today 37, 62 (1984b).
  * Mozurkewich and Berry [1981] M. Mozurkewich and R. S. Berry, Finite-time thermodynamics: Engine performance improved by optimized piston motion, Proc. Natl. Acad. Sci. U.S.A. 78, 1986 (1981).
  * Van den Broeck [2005] C. Van den Broeck, Thermodynamic efficiency at maximum power, Phys. Rev. Lett. 95, 190602 (2005).
  * Jiménez de Cisneros and Hernández [2007] B. Jiménez de Cisneros and A. C. Hernández, Collective working regimes for coupled heat engines, Phys. Rev. Lett. 98, 130602 (2007).
  * Hoffmann [2008] K. H. Hoffmann, An introduction to endoreversible thermodynamics, AAPP Phys. Math. Nat. Sci. 86, 1 (2008).
  * Holubec and Ryabov [2021] V. Holubec and A. Ryabov, Fluctuations in heat engines, Phys. A Math. Theor. 55, 013001 (2021).
  * Geva and Kosloff [1992] E. Geva and R. Kosloff, A quantum-mechanical heat engine operating in finite time. a model consisting of spin-1/2 systems as the working fluid, J. Chem. Phys. 96, 3054 (1992).
  * Parmeggiani  _et al._ [1999] A. Parmeggiani, F. Jülicher, A. Ajdari, and J. Prost, Energy transduction of isothermal ratchets: Generic aspects and specific examples close to and far from equilibrium, Phys. Rev. E 60, 2127 (1999).
  * Hondou and Sekimoto [2000] T. Hondou and K. Sekimoto, Unattainability of carnot efficiency in the brownian heat engine, Phys. Rev. E 62, 6021 (2000).
  * 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).
  * Astumian and Hänggi [2002] R. D. Astumian and P. Hänggi, Brownian motors, Phys. Today 55, 33 (2002).
  * Parrondo and de Cisneros [2002] J. M. R. Parrondo and B. J. de Cisneros, Energetics of brownian motors: a review, Appl. Phys. A 75, 179 (2002).
  * Reimann [2002] P. Reimann, Brownian motors: noisy transport far from equilibrium, Phys. Rep. 361, 57 (2002).
  * Schmiedl and Seifert [2007a] T. Schmiedl and U. Seifert, Efficiency at maximum power: An analytically solvable model for stochastic heat engines, EPL 81, 20003 (2007a).
  * Tu [2008] Z. C. Tu, Efficiency at maximum power of feynman’s ratchet as a heat engine, J. Phys. A Math. Theor. 41, 312003 (2008).
  * Esposito  _et al._ [2009a] M. Esposito, K. Lindenberg, and C. Van den Broeck, Universality of efficiency at maximum power, Phys. Rev. Lett. 102, 130602 (2009a).
  * Esposito  _et al._ [2009b] M. Esposito, K. Lindenberg, and C. V. den Broeck, Thermoelectric efficiency at maximum power in a quantum dot, EPL 85, 60010 (2009b).
  * Esposito  _et al._ [2010a] M. Esposito, R. Kawai, K. Lindenberg, and C. Van den Broeck, Efficiency at maximum power of low-dissipation carnot engines, Phys. Rev. Lett. 105, 150603 (2010a).
  * 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).
  * Raz  _et al._ [2016] O. Raz, Y. Subaşı, and R. Pugatch, Geometric heat engines featuring power that grows with efficiency, Phys. Rev. Lett. 116, 160601 (2016).
  * Holubec and Ryabov [2018] V. Holubec and A. Ryabov, Cycling tames power fluctuations near optimum efficiency, Phys. Rev. Lett. 121, 120601 (2018).
  * Brandner and Saito [2020] K. Brandner and K. Saito, Thermodynamic geometry of microscopic heat engines, Phys. Rev. Lett. 124, 040602 (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).
  * 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).
  * Erdman  _et al._ [2019] P. A. Erdman, V. Cavina, R. Fazio, F. Taddei, and V. Giovannetti, Maximum power and corresponding efficiency for two-level heat engines and refrigerators: optimality of fast cycles, New J. Phys. 21, 103049 (2019).
  * Das and Mukherjee [2020] A. Das and V. Mukherjee, Quantum-enhanced finite-time otto cycle, Phys. Rev. Research 2, 033083 (2020).
  * Esposito  _et al._ [2010b] M. Esposito, R. Kawai, K. Lindenberg, and C. Van den Broeck, Efficiency at maximum power of low-dissipation carnot engines, Phys. Rev. Lett. 105, 150603 (2010b).
  * Abiuso and Perarnau-Llobet [2020] P. Abiuso and M. Perarnau-Llobet, Optimal cycles for low-dissipation heat engines, Phys. Rev. Lett. 124, 110606 (2020).
  * Humphrey  _et al._ [2002] T. E. Humphrey, R. Newbury, R. P. Taylor, and H. Linke, Reversible quantum brownian heat engines for electrons, Phys. Rev. Lett. 89, 116801 (2002).
  * Abiuso  _et al._ [2020] P. Abiuso, H. J. Miller, M. Perarnau-Llobet, and M. Scandi, Geometric optimisation of quantum thermodynamic processes, Entropy 22, 1076 (2020).
  * 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).
  * 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).
  * Holubec [2014a] V. Holubec, An exactly solvable model of a stochastic heat engine: optimization of power, power fluctuations and efficiency, J. Stat. Mech. Theory Exp. 2014, P05022 (2014a).
  * [44] See Supplemental Information at http://link.aps.org/.
  * Dechant  _et al._ [2017] A. Dechant, N. Kiesel, and E. Lutz, Underdamped stochastic heat engine at maximum efficiency, EPL 119, 50003 (2017).
  * Zhang [2020] Y. Zhang, Optimization of stochastic thermodynamic machines, J. Stat. Phys. 178, 1336 (2020).
  * Abiuso  _et al._ [2022] P. Abiuso, V. Holubec, J. Anders, Z. Ye, F. Cerisola, and M. Perarnau-Llobet, Thermodynamics and optimal protocols of multidimensional quadratic brownian systems, J. Phys. Commun. 6, 063001 (2022).
  * Zhong and DeWeese [2022] A. Zhong and M. R. DeWeese, Limited-control optimal protocols arbitrarily far from equilibrium, arXiv preprint arXiv:2205.08662 (2022).
  * Bauer  _et al._ [2016] M. Bauer, K. Brandner, and U. Seifert, Optimal performance of periodically driven, stochastic heat engines under limited control, Phys. Rev. E 93, 042112 (2016).
  * Plata  _et al._ [2019] C. A. Plata, D. Guéry-Odelin, E. Trizac, and A. Prados, Optimal work in a harmonic trap with bounded stiffness, Phys. Rev. E 99, 012140 (2019).
  * Risken [1996] H. Risken, Fokker-planck equation, in  _The Fokker-Planck Equation: Methods of Solution and Applications_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 1996) pp. 63–95.
  * Holubec  _et al._ [2022] V. Holubec, A. Ryabov, S. A. M. Loos, and K. Kroy, Equilibrium stochastic delay processes, New J. Phys. 24, 023021 (2022).
  * Krishnamurthy  _et al._ [2016] S. Krishnamurthy, S. Ghosh, D. Chatterji, R. Ganapathy, and A. K. Sood, A micrometre-sized heat engine operating between bacterial reservoirs, Nat. Phys. 12, 1134 (2016).
  * Holubec  _et al._ [2020] V. Holubec, S. Steffenoni, G. Falasco, and K. Kroy, Active brownian heat engines, Phys. Rev. Research 2, 043262 (2020).
  * Kumari  _et al._ [2020] A. Kumari, P. S. Pal, A. Saha, and S. Lahiri, Stochastic heat engine using an active particle, Phys. Rev. E 101, 032109 (2020).
  * Gronchi and Puglisi [2021] G. Gronchi and A. Puglisi, Optimization of an active heat engine, Phys. Rev. E 103, 052134 (2021).
  * Note [1] A similar optimisation problem is often solved in courses on classical thermodynamics to show that maximum efficiency of an equilibrium cycle under the constrains $T(t)\in[T_{-},T_{+}]$ on the bath temperature is the Carnot efficiency. However, in our case, the system can be arbitrarily far form equilibrium.
  * Martínez  _et al._ [2015] I. A. Martínez, E. Roldán, L. Dinis, D. Petrov, and R. A. Rica, Adiabatic processes realized with a trapped brownian particle, Phys. Rev. Lett. 114, 120601 (2015).
  * Arold  _et al._ [2018] D. Arold, A. Dechant, and E. Lutz, Heat leakage in overdamped harmonic systems, Phys. Rev. E 97, 022131 (2018).
  * Blickle and Bechinger [2012] V. Blickle and C. Bechinger, Realization of a micrometre-sized stochastic heat engine, Nat. Phys. 8, 143 (2012).
  * Holubec and Ryabov [2017] V. Holubec and A. Ryabov, Diverging, but negligible power at carnot efficiency: Theory and experiment, Phys. Rev. E 96, 062107 (2017).
  * Cavina  _et al._ [2018] V. Cavina, A. Mari, A. Carlini, and V. Giovannetti, Optimal thermodynamic control in open quantum systems, Phys. Rev. A 98, 012139 (2018).
  * Martínez  _et al._ [2015] I. A. Martínez, É. Roldán, L. Dinis, D. Petrov, J. M. R. Parrondo, and R. A. Rica, Brownian carnot engine, Nat. Phys. 12, 67 (2015).
  * Martínez  _et al._ [2017] I. A. Martínez, É. Roldán, L. Dinis, and R. A. Rica, Colloidal heat engines: a review, Soft Matter 13, 22 (2017).
  * Then and Engel [2008] H. Then and A. Engel, Computing the optimal protocol for finite-time processes in stochastic thermodynamics, Phys. Rev. E 77, 041105 (2008).
  * Holubec and Ryabov [2015] V. Holubec and A. Ryabov, Efficiency at and near maximum power of low-dissipation heat engines, Phys. Rev. E 92, 052125 (2015).
  * Gingrich  _et al._ [2016] T. R. Gingrich, G. M. Rotskoff, G. E. Crooks, and P. L. Geissler, Near-optimal protocols in complex nonequilibrium transformations, Proc. Natl. Acad. Sci. U.S.A. 113, 10263 (2016).
  * 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).
  * Khait  _et al._ [2022] I. Khait, J. Carrasquilla, and D. Segal, Optimal control of quantum thermal machines using machine learning, Phys. Rev. Research 4, L012029 (2022).
  * Dechant and Sakurai [2019] A. Dechant and Y. Sakurai, Thermodynamic interpretation of wasserstein distance, arXiv preprint arXiv:1912.08405 (2019).
  * Van Vu and Hasegawa [2021] T. Van Vu and Y. Hasegawa, Geometrical bounds of the irreversibility in markovian systems, Phys. Rev. Lett. 126, 010601 (2021).
  * Carrega  _et al._ [2022] M. Carrega, L. M. Cangemi, G. De Filippis, V. Cataudella, G. Benenti, and M. Sassetti, Engineering dynamical couplings for quantum thermodynamic tasks, PRX Quantum 3, 010323 (2022).
  * Pancotti  _et al._ [2020] N. Pancotti, M. Scandi, M. T. Mitchison, and M. Perarnau-Llobet, Speed-ups to isothermality: Enhanced quantum thermal machines through control of the system-bath coupling, Phys. Rev. X 10, 031015 (2020).
  * Holubec [2014b] V. Holubec,  _Non-equilibrium Energy Transformation Processes: Theoretical Description at the Level of Molecular Structures_, Springer Theses (Springer International Publishing, 2014).
  * Schmiedl and Seifert [2007b] T. Schmiedl and U. Seifert, Optimal finite-time processes in stochastic thermodynamics, Phys. Rev. Lett. 98, 108301 (2007b).
  * Holubec and Ryabov [2016] V. Holubec and A. Ryabov, Maximum efficiency of low-dissipation heat engines at arbitrary power, Journal of Statistical Mechanics: Theory and Experiment 2016, 073204 (2016).

##  Appendix A Maximum-efficiency protocol for multi-term Hamiltonians

Consider a heat engine with a working fluid described by the Hamiltonian

$H(x,t)=\sum_{i}\lambda_{i}(t)f_{i}(x)$ (6)

with control parameters $\lambda_{i}(t)$, $i=1,\dots,N$. As in the main text,
we now aim to derive the finite-time protocol for the constrained control
parameters, $\lambda_{i}(t)\in(\lambda_{i}^{-},\lambda_{i}^{+})$, which would
yield maximum efficiency of the engine. It will turn out that if the
compression ratios $\lambda_{i}^{-}/\lambda_{i}^{+}$ for all the control
parameters equal, the geometric argument from the main text still applies.

The heat increment is for the Hamiltonian (6) given by
đ$Q=\sum_{i}\lambda_{i}(t)d\sigma_{i}(t)$ with the response functions
$\sigma_{i}(t)=\langle f_{i}(x)\rangle$. For arbitrary fixed maximum changes
$\Delta\sigma_{i}$ in the response functions during the cycle, geometric upper
and lower bounds on $Q_{\rm in}$ and $Q_{\rm out}$ and thus on efficiency are
achieved by clockwise rectangular cycles in the individual
$\lambda_{i}-\sigma_{i}$ diagrams. These hypothetical cycles yield the
following geometric upper bound on efficiency:

$\eta=1-\frac{Q_{\rm out}}{Q_{\rm in}}\leq
1-\frac{\sum_{i}\Delta\sigma_{i}\lambda_{i}^{-}}{\sum_{i}\Delta\sigma_{i}\lambda_{i}^{+}}.$
(7)

We use the term ‘geometric’ to stress that this bound follows from the
analysis of the cycle in the $\lambda-\sigma$ diagram, without considering the
relation between the protocol $[\lambda_{1}(t),\dots,\lambda_{N}(t)]$ and the
response $[\sigma_{1}(t),\dots,\sigma_{N}(t)]$ imposed by dynamical equations
of the working fluid. This means that the given set of $\Delta\sigma_{i}$
might not be achievable by the piecewise constant protocol and thus the bound
in (7) is loose. Furthermore, we seek an optimal protocol constrained just by
the conditions on $\lambda_{i}$ and the upper bound in (7) in general strongly
depends on the fixed values of $\Delta\sigma_{i}$. For the single-term
Hamiltonian $H(x,t)=\lambda(t)f(x)$ used in the main text, this has not been
an issue because then $\Delta\sigma$ in the nominator and denominator in (7)
cancel out and the upper bound becomes independent of the details of the
dynamics. The optimal protocol for efficiency is then the piecewise constant
protocol for $\lambda(t)$ because it saturates the geometric upper bound. To
sum up, the bound in (7) allows to derive the maximum-efficiency protocol only
if it happens to be independent of $\Delta\sigma_{i}$. In the opposite case,
the optimal protocol cannot be determined without considering the dynamical
equations and performing the corresponding functional optimisation.

Let us now investigate when the upper bound in (7) becomes independent of the
system response, $\Delta\sigma_{i}$. Defining the set of ‘probabilities’
$p_{i}=\Delta\sigma_{i}\lambda_{i}^{+}/\sum_{i}\Delta\sigma_{i}\lambda_{i}^{+}$,
the ratio in the upper bound in (7) can be rewritten as the average

$\frac{\sum_{i}\Delta\sigma_{i}\lambda_{i}^{-}}{\sum_{i}\Delta\sigma_{i}\lambda_{i}^{+}}=\sum_{i}p_{i}\frac{\lambda_{i}^{-}}{\lambda_{i}^{+}}.$
(8)

This expression becomes independent of $\sigma_{i}$ only if all the
compression ratios $\lambda_{i}^{-}/\lambda_{i}^{+}$ are equal. In such a
case, the maximum-efficiency protocol is thus a piecewise constant protocol
for each of $\lambda_{i}$ and yields the efficiency

$\eta=1-\lambda_{i}^{-}/\lambda_{i}^{+}.$ (9)

Besides this result, the probabilistic interpretation (8) of the upper bound
in (7) also yields the dynamics independent (but in general loose) _upper
bound_ on efficiency,

$\eta\leq 1-\min_{i}\frac{\lambda_{i}^{-}}{\lambda_{i}^{+}}.$ (10)

To close this section, we note that a piecewise constant protocol for
$\lambda_{i}$ will always yield the efficiency
$1-(\sum_{i}\Delta\sigma_{i}\lambda_{i}^{-})/(\sum_{i}\Delta\sigma_{i}\lambda_{i}^{+}),$
with values of $\Delta\sigma_{i}$ induced by the dynamical equations of the
system. Within the class of piecewise constant protocols, the upper bound (10)
is then tight if the constraints on all the control parameters $\lambda_{i}$
allow to achieve the minimum compression ratio
$\min_{i}\frac{\lambda_{i}^{-}}{\lambda_{i}^{+}}$. Furthermore, for such
protocols, Eq. (8) also implies the _lower bound_ on the efficiency

$\eta\geq 1-\max_{i}\frac{\lambda_{i}^{-}}{\lambda_{i}^{+}},$ (11)

which is always tight.

##  Appendix B Properties of maximum-efficiency protocol

In this section, we provide further details concerning the maximum-efficiency
protocol for the Hamiltonian, $H(x,t)=\lambda(t)f(x)$, discussed in the main
text. First, we argue that the maximum-efficiency protocol that yields maximum
output work for the given piecewise constant $\lambda(t)$ requires piecewise
constant variation of temperature. Then, we investigate output power of the
maximum-efficiency protocol as a function of the lower bound on the control
parameter $\lambda(t)$.

###  B.1 Temperature protocol

In the main text, we have shown that the maximum-efficiency protocol for the
control parameter $\lambda(t)$ is piecewise constant and the corresponding
efficiency $\eta=1-\lambda_{-}/\lambda_{+}$. The only condition on the
temperature protocol was that the cycle is performed clockwise in the
$\lambda-\sigma$ diagram. Nevertheless, in order to allow the engine to
operate at Carnot efficiency and to maximize its output work, we have chosen
the protocol

$\\{T(t),\lambda(t)\\}_{\eta}=\begin{cases}\\{T_{+},\lambda_{+}\\},&0<t<t_{+}\\\
\\{T_{-},\lambda_{-}\\},&t_{+}<t<t_{\rm p}\end{cases}.$ (12)

For this choice of $T(t)$, the working medium of the engine operates with the
largest possible temperature gradient during the whole cycle. This maximizes
the heat flux through the engine, which can be utilised to yield the maximum
amount of work $W_{\rm out}=\eta Q_{\rm in}$. Besides, the engine efficiency
$\eta$ is also known to increase with the bath temperature difference [see
also Fig. 7 (c)-(f)].

Let us now provide an alternative and more technical argument that the choice
of $T(t)$ in Eq. (12) maximizes the output work. We restrict this argument to
the maximum-efficiency protocol for $\lambda$ in Eq. (12). However,
generalizations to other protocols are straightforward. The main idea is that
connecting the system to the hottest possible bath when $\dot{\sigma}>0$ and
to the coldest possible bath when $\dot{\sigma}<0$ maximizes the extent of the
cycle in the $\sigma$ direction in the $\sigma-\lambda$ diagram and thus also
$W_{\rm out}$.

For the protocol (12), the output work is given by

$W_{\rm out}=\Delta\lambda\Delta\sigma,$ (13)

with $\Delta\lambda=\lambda_{+}-\lambda_{-}$ and the maximum change in the
response parameter during the cycle $\Delta\sigma=\sigma_{+}-\sigma_{-}$. To
maximize $W_{\rm out}$, we thus need to maximize $\Delta\sigma$. To this end,
it is reasonable to assume that

$\Delta\sigma\leq\Delta\sigma^{\rm eq},$ (14)

where $\Delta\sigma^{\rm eq}=\max\sigma^{\rm eq}-\min\sigma^{\rm eq}$ is the
maximum change in the response parameter $\sigma$ during the cycle with
isochoric branches (constant $\lambda$) longer than the system relaxation
time. This assumption is in particular valid for arbitrary overdamped
dynamics, where $\sigma$ always converges to its equilibrium value ($k_{\rm
B}$ denotes the Boltzmann constant)

$\sigma^{\rm eq}(t)=\sum_{x}f(x)\frac{\exp\left\\{-\lambda(t)f(x)/[k_{\rm
B}T(t)]\right\\}}{\sum_{x}\exp\left\\{-\lambda(t)f(x)/[k_{\rm
B}T(t)]\right\\}},$ (15)

corresponding to the instantaneous values of the control parameters
$\\{T(t),\lambda(t)\\}$. Noticing that $\sigma^{\rm eq}(t)=U(t)/\lambda(t)$,
where $U(t)=\left<H(x,t)\right>$ is the thermodynamic internal energy of the
system, the positivity of heat capacity

$C_{\rm v}=\frac{\partial{U}}{\partial{T}}=\frac{\partial{\sigma^{\rm
eq}}}{\partial(T/\lambda)}>0$ (16)

implies that $\sigma^{\rm eq}$ is a monotonously increasing function of the
ratio $T/\lambda$.

From Fig. 1 in the main text, it follows that $\max\sigma^{\rm eq}$ and
$\min\sigma^{\rm eq}$ are the values of $\sigma^{\rm eq}$ at the ends of the
isochores with $\lambda=\lambda_{+}$ and $\lambda=\lambda_{-}$, respectively.
The upper bound on $\Delta\sigma$ is thus given by

$\max\sigma^{\rm eq}-\min\sigma^{\rm eq}=\sigma^{\rm
eq}(T_{+}/\lambda_{+})-\sigma^{\rm eq}(T_{-}/\lambda_{-}).$ (17)

It is attained for slow isochores when $T=T_{+}$ for $\lambda=\lambda_{+}$ and
$T=T_{-}$ for $\lambda=\lambda_{-}$. As long as $\dot{\sigma}^{\rm eq}>0$ for
$\lambda=\lambda_{+}$ and $\dot{\sigma}^{\rm eq}<0$ for $\lambda=\lambda_{-}$
(so that the used definitions of input and output heat hold), details of the
temperature protocol during the isochores in this limit do not alter the value
of $\Delta\sigma^{\rm eq}$ and thus $W_{\rm
out}=\Delta\lambda\Delta\sigma^{\rm eq}$. However, these details become
important for finite-time cycles.

A typical dynamical equation for an overdamped degree of freedom has the form

$\dot{\sigma}(t)=t_{\rm R}^{-1}[\sigma^{\rm eq}-\sigma(t)].$ (18)

For constant values of control parameters $T(t)$ and $\lambda(t)$, which enter
the relaxation time $t_{\rm R}$ and the equilibrium state $\sigma^{\rm eq}(t)$
defined in Eq. (15), this equation describes an exponential relaxation of
$\sigma$ to $\sigma^{\rm eq}$ (for a specific example, see Sec. B.2.2). For a
cyclic variation of the control parameters, $\sigma$ follows $\sigma^{\rm eq}$
like a donkey chasing a carrot [74]. As a result, $\sigma\leq\sigma^{\rm eq}$
and $\dot{\sigma}\geq 0$ for $\lambda=\lambda_{+}$, when $\sigma^{\rm eq}$
increases to $\max\sigma^{\rm eq}$, and $\sigma\geq\sigma^{\rm eq}$ and
$\dot{\sigma}\leq 0$ for $\lambda=\lambda_{-}$, when $\sigma^{\rm eq}$
decreases to $\min\sigma^{\rm eq}$. The change in the response
$\Delta\sigma=\int_{0}^{t_{+}}\dot{\sigma}dt=-\int_{t_{+}}^{t_{\rm
p}}\dot{\sigma}dt$ and thus it can be maximized by maximizing (minimizing) the
instantaneous rate of change of the response, $\dot{\sigma}$, during the first
(second) isochore. From Eq. (18), it follows that this is achieved by setting
$\sigma^{\rm eq}=\max\sigma^{\rm eq}$ during the fist isochore and
$\sigma^{\rm eq}=\min\sigma^{\rm eq}$ during the second one. Altogether, this
suggests that the piecewise constant temperature protocol in Eq. (12) yields
maximum $\Delta\sigma$ and thus output work $W_{\rm out}$ (13) for arbitrary
cycle duration.

Figure 4: Efficiency at maximum output work obtained using the Hamiltonian
$H=\lambda(t)(|x|^{n}/n-\ln|x|)$ as a function of $T_{-}/T_{+}$. Parameters
used: $k_{\rm B}T_{+}=1$ and $\lambda_{+}=0.5$.

###  B.2 Maximum output work of maximum-efficiency protocol

Let us now turn to the task of maximizing the output work $W_{\rm
out}=(\lambda_{+}-\lambda_{-})\Delta\sigma$ with respect to $\lambda_{-}$.
Analytical results can be obtained in the limits of slow and fast isotherms.

####  B.2.1 Slow isotherms

When the duration of the isotherms is longer than the relaxation time of the
response $\sigma$, one can approximate $\sigma_{+}$ and $\sigma_{-}$ in
$\Delta\sigma$ by their equilibrium values. Using Eq. (13), the output work
then reads

$W_{\rm out}=\Delta\lambda\Delta\sigma^{\rm eq}.$ (19)

Equation (16) implies that the partial derivative of $\sigma^{\rm eq}$ with
respect to the control parameter $\lambda$ ($T$ is constant) is given by

$\frac{\partial}{\partial\lambda}\sigma^{\rm
eq}(T/\lambda)=-\frac{T}{\lambda^{2}}C_{\rm v}.$ (20)

The condition on the extreme of $W_{\rm out}$ (19) with respect to
$\lambda_{-}$ thus reads:

$\begin{split}\frac{\partial{W_{\rm
out}}}{\partial\lambda_{-}}&=(\lambda_{+}-\lambda_{-})\frac{T_{-}}{\lambda_{-}^{2}}C_{\rm
v}(T_{-}/\lambda_{-})\\\
&-\frac{U(T_{+}/\lambda_{+})}{\lambda_{+}}+\frac{U(T_{-}/\lambda_{-})}{\lambda_{-}}=0,\end{split}$
(21)

where we additionally used the relation $\sigma^{\rm eq}=U/\lambda$ between
$\sigma^{\rm eq}$ and the internal energy $U$.

For power law Hamiltonians of the form $H=\lambda|x|^{n}/n$ where $C_{\rm
v}=k_{\rm B}/n$ and $U=k_{\rm B}T/n$, this equation can be solved explicitly.
The resulting optimal compression ratio is given by
$\lambda_{-}/\lambda_{+}=\sqrt{T_{-}/T_{+}}$. The corresponding efficiency at
the maximum output work is given by the Curzohn-Ahlborn efficiency,

$\eta=1-\frac{\lambda_{-}}{\lambda_{+}}=1-\sqrt{\frac{T_{-}}{T_{+}}}\equiv\eta_{\rm
CA},$ (22)

and the maximum output work is (Carnot efficiency $\eta_{\rm
C}=1-T_{-}/{T_{+}}$)

$W_{\rm out}=\frac{k_{\rm B}T_{+}}{n}(2\eta_{\rm CA}-\eta_{\rm C}).$ (23)

Let us now consider the asymmetric Hamiltonian
$H=\lambda(t)(|x|^{n}/n-\ln|x|)$. In this case, the internal energy and heat
capacity are given by

$\displaystyle U$ $\displaystyle=$ $\displaystyle\frac{k_{\rm
B}T+\lambda\left[1+\ln{\frac{\lambda}{nk_{\rm
B}T}}-\psi^{(0)}\left(\frac{\lambda+k_{\rm B}T}{nk_{\rm
B}T}\right)\right]}{n},$ (24) $\displaystyle C_{\rm v}$ $\displaystyle=$
$\displaystyle\frac{nk_{\rm B}T(k_{\rm
B}T-\lambda)+\lambda^{2}\psi^{(1)}\left(\frac{\lambda+k_{\rm B}T}{nk_{\rm
B}T}\right)}{n^{2}k^{2}_{\rm B}T^{2}},$ (25)

where $\psi^{(m)}(z)$ denotes the polygamma function of order $m$. In this
case, Eq. (21) is transcendental and we solved it numerically. In Fig. 4, we
show the resulting efficiency at the maximum output work as a function of
$T_{-}/T_{+}$. Even though the resulting efficiency is still close to
$\eta_{\rm CA}$, it can be both slightly larger and smaller than that.

####  B.2.2 Fast isotherms

Let us now assume that the duration of the isothermal branches are much
shorter than the system relaxation time. In such a case, the work optimisation
cannot be done without specifying the dynamical equation for the response
$\sigma$. To this end, we assume that it obeys the overdamped equation

$\dot{\sigma}(t)=t_{\rm R}^{-1}[\sigma^{\rm eq}-\sigma(t)]$ (26)

with the equilibrium value $\sigma^{\rm eq}$ and relaxation time $t_{R}$
determined by the values of the control parameters $\\{T(t),\lambda(t)\\}$ at
time $t$. The most prominent examples of systems described by this formula is
a two-level system [35] and an overdamped particle trapped in a harmonic
potential [24].

Solving Eq. (26) for the maximum-efficiency protocol (12), we find that

$\sigma(t)=\begin{cases}\sigma_{0}{\rm e}^{-\frac{t}{t_{R}^{+}}}+\sigma_{\rm
eq}^{+}(1-{\rm e}^{-\frac{t}{t_{R}^{+}}}),&0<t<t_{+}\\\ \sigma_{1}{\rm
e}^{-\frac{t-t_{+}}{t_{R}^{-}}}+\sigma_{\rm eq}^{-}(1-{\rm
e}^{-\frac{t-t_{+}}{t_{R}^{-}}}),&t_{+}<t<t_{\rm p}\end{cases},$ (27)

where $\sigma_{0}\equiv\sigma(0)$ and $\sigma_{1}\equiv\sigma(t_{+})$ are
determined by the condition that $\sigma(t)$ must be a continuous function of
time. The variables corresponding to the first (second) isotherm are denoted
by max (min). It turns out that

$\Delta\sigma=\sigma_{+}-\sigma_{-}=\sigma_{1}-\sigma_{0}=\\\
\Delta\sigma^{\rm eq}\frac{\sinh\left(t_{+}/t_{\rm
R}^{+}\right)\sinh\left(t_{-}/t_{\rm R}^{-}\right)}{\sinh\left(t_{+}/t_{\rm
R}^{+}+t_{-}/t_{\rm R}^{-}\right)}\leq\Delta\sigma^{\rm eq}.$ (28)

Substituting this result into the expression for the output work (13) and
expanding the result up to the leading order in the ratios of duration of the
individual isotherms to the corresponding relaxation times, $t_{+}/t_{\rm
R}^{+}$ and $t_{-}/t_{\rm R}^{-}$, we find that

$W_{\rm out}=\Delta\lambda\Delta\sigma^{\rm eq}\frac{\frac{t_{+}t_{-}}{t_{\rm
R}^{+}t_{\rm R}^{-}}}{\frac{t_{+}}{t_{\rm R}^{+}}+\frac{t_{-}}{t_{\rm
R}^{-}}}.$ (29)

To maximize the output work, we need to choose a specific model to determine
the dependence of the equilibrium values of response and relaxation times on
the control parameters. To this end, we consider the paradigmatic model of
stochastic thermodynamics,
$\dot{\sigma}(t)=-2\,\mu\,\lambda(t)\,\sigma(t)+\,\mu\,k_{\rm B}T$, describing
an overdamped Brownian particle with mobility $\mu$ in a harmonic trap [24,
75]. In this case, $\sigma^{\rm eq}=T/(2\mu\lambda)$ and
$t_{R}=1/(2\mu\lambda)$, and the maximum output work (29) is produced for

$\frac{\lambda_{-}}{\lambda_{+}}=\sqrt{(\alpha+1)\left(\alpha+1-\eta_{\rm
C}\right)}-\alpha,$ (30)

where $\alpha\equiv t_{+}/t_{-}$. The corresponding efficiency reads

$\eta=1-\frac{\lambda_{-}}{\lambda_{+}}=\alpha+1-\sqrt{(\alpha+1)(\alpha+1-\eta_{\rm
C})},$ (31)

which reduces to $\eta_{\rm CA}$ for $\alpha\to 0$ and $\eta_{\rm C}/2$ for
$\alpha\to\infty$. Assuming that $\alpha=1$ ($t_{+}=t_{-}$), Eq. (31) is given
by the formula

$\eta=2-\sqrt{4-2\eta_{\rm C}}=\frac{\eta_{\rm C}}{2}+\frac{\eta^{2}_{\rm
C}}{16}+\mathcal{O}(\eta^{3}_{\rm C})$ (32)

used in the main text. The corresponding expansion for the Curzohn-Ahlborn
efficiency, $\eta_{\rm CA}\approx\frac{\eta_{\rm C}}{2}+\frac{2\eta_{\rm
C}^{2}}{16}$, has identical linear and twice larger quadratic term. Last but
not least, with respect to $\alpha$, the output power $W_{\rm out}/t_{\rm p}$
using Eq. (29) develops a peak at
$\alpha=\alpha^{*}=\sqrt{\frac{\lambda_{-}}{\lambda_{+}}}<1$. This also
contradicts the situation with constrained $\sigma$, where maximum power is
attained when the durations of the isotherms are equal ($\alpha=\alpha^{*}=1$)
[24].

##  Appendix C Numerical results

In this section, we give further details concerning the numerical optimisation
of output work (and power) described in the main text. All the numerical
results were obtained using the model [24, 75]

$\dot{\sigma}(t)=-2\,\mu\,\lambda(t)\,\sigma(t)+\,\mu\,k_{\rm B}T$ (33)

for an overdamped Brownian particle with mobility $\mu$ in a harmonic trap
with stiffness $\lambda(t)$, which well describes experiments with optical
tweezers. This equation is exactly solvable for arbitrary protocol
$\\{T(t),\lambda(t)\\}$. However, the solutions can involve exponentials of
very large or small numbers, which can lead to numerical instabilities
inducing large looses of precision, and thus they have to be treated with
care. To secure that our solutions to Eq. (33) are always enough precise, we
have solved it in our analysis also numerically.

For $\lambda(t)$, we considered the maximum-efficiency protocol (12) and the
four classes of protocols in Tab. I in the main text. For $T(t)$, we always
used the protocol in Eq. (12). For these protocols, we optimised the output
power and efficiency

$P=\frac{Q_{\rm in}-Q_{\rm out}}{t_{\rm p}},\quad\eta=1-\frac{Q_{\rm
out}}{Q_{\rm in}},$ (34)

with the input heat $Q_{\rm in}=\int_{0}^{t_{+}}dt\dot{\sigma}(t)\lambda(t)$
and the output heat $Q_{\rm out}=-\int_{t_{+}}^{t_{\rm
p}}dt\dot{\sigma}(t)\lambda(t)$,  as functions of the parameters
$\\{a,b,c,d\\}$ (see Tab. I in the main text) under different constraints.

In Sec. C.1, we verify that the protocol $\lambda_{\rm S}$ from Ref. [24]
yields both maximum power and efficiency under constrains on the response
$\sigma$. In Sec. C.2, we provide numerical results illustrating that the
protocol (12) yields maximum efficiency under constrains on the response
$\lambda$. Furthermore, we present additional evidence that the maximum-
efficiency protocol (12) optimized with respect to $\lambda_{-}$ can be used
to approximately obtain also the maximum power.

###  C.1 Constrained response

Figure 5: Optimal performance for fixed boundary values of the response:
$\sigma_{-}\equiv\sigma(0)=\sigma(t_{\rm p})$ and
$\sigma_{+}\equiv\sigma(t_{+})=0.5$. (a) Maximum power (in units of the
ultimate maximum power $P^{*}$ for $\lambda_{\rm pwl}$) and (b) maximum
efficiency as functions of $\sigma_{-}/\sigma_{+}$. Lines corresponding to
$\lambda_{\rm S}$ (cyan dashed) and $\lambda_{\rm slow}$ (blue dotted)
perfectly overlap. The maximum-efficiency protocol (12) and the piecewise
constant protocol $\lambda_{\rm pwc}$ are in this case equal. We used the same
parameters as in Fig. 2 in the main text.

When the values of the response (position variance) $\sigma$ at the ends of
the two isotherms are fixed, i.e., $\sigma_{-}\equiv\sigma(0)=\sigma(t_{\rm
p})$ and $\sigma_{+}\equiv\sigma(t_{+})$, the protocol which yields both
maximum efficiency and power reads [24]

$\lambda_{\rm S}=\begin{cases}\frac{T_{+}}{2\sigma_{\rm
S}}-\frac{\sqrt{\sigma_{+}}-\sqrt{\sigma_{-}}}{\mu t_{+}\sqrt{\sigma_{\rm
S}}},~{}~{}0<t<t_{+}\\\ \frac{T_{-}}{2\sigma_{\rm
S}}+\frac{\sqrt{\sigma_{+}}-\sqrt{\sigma_{-}}}{\mu t_{-}\sqrt{\sigma_{\rm
S}}},~{}~{}t_{+}<t<t_{\rm p}\end{cases}$ (35)

with

$\sigma_{\rm
S}=\begin{cases}\frac{\sigma_{-}}{2}\left[1+(\sqrt{\frac{\sigma_{+}}{\sigma_{-}}}-1)\frac{t}{t_{+}}\right]^{2},0<t<t_{+}\\\
\frac{\sigma_{+}}{2}\left[1+(\sqrt{\frac{\sigma_{-}}{\sigma_{+}}}-1)\frac{t-t_{+}}{t_{-}}\right]^{2},t_{+}<t<t_{\rm
p}\end{cases}.$ (36)

However, this protocol is no longer optimal when one imposes just maximum and
minimum values on the response, i.e., $\sigma(t)\in[\sigma_{-},\sigma_{+}]$.
Then, our analysis shows that the maximum-efficiency and maximum-power
protocol is still of the above form, but with
$\sigma_{-}<\sigma(0)=\sigma(t_{\rm p})<\sigma(t_{+})<\sigma_{+}$.

In Fig. 5, we show the maximum power (a) and maximum efficiency (b) for the
trial protocols under the constraint $\sigma_{-}\equiv\sigma(0)=\sigma(t_{\rm
p})$ and $\sigma_{+}\equiv\sigma(t_{+})$. As expected, power and efficiency
corresponding to the protocol $\lambda_{\rm S}$ are largest from all the
protocols. In particular, the figure demonstrates that the linear protocol,
which was found to maximize output power for constrained $\lambda$, yields
smaller output power than $\lambda_{\rm S}$. And the piecewise constant
protocol yields smaller efficiency than $\lambda_{\rm S}$. Nevertheless, it is
interesting to note that the performance of the protocol
$\lambda_{\mathrm{slow}}(t)$, which optimizes both output power and efficiency
for slow driving (see Sec. D for details), is for the chosen parameters
indistinguishable from that of $\lambda_{\rm S}$. This means that the chosen
cycle is slow enough. Finally, for small enough cycles (small
$\sigma_{-}/\sigma_{+}$) performance of all the protocols equal.

Figure 6: (a) Maximum efficiency and (b) the corresponding power (in units of
the ultimate maximum power $P^{*}$ for $\lambda_{\rm pwl}$) as functions of
$\lambda_{-}/\lambda_{+}$.  All protocols except for $\lambda_{\rm S}$
perfectly overlap. We used the same parameters as in Fig. 2 in the main text.
Figure 7:  The relative differences $\delta X=(X_{\lambda_{\rm
pwl}}-X_{\lambda_{\eta}})/X_{\lambda_{\rm pwl}}$ of (a) maximum power ($X=P$)
optimized with respect to $\lambda_{-}$ and (b) the corresponding efficiency
($X=\eta$) for the linear protocol $\lambda_{\rm pwl}$ and the maximum-
efficiency protocol (12) for different values of $T_{-}/T_{+}$ and
$t_{-}/t_{+}$. (c)-(f) show the corresponding values of maximum power and
efficiency. The piecewise constant protocol $\lambda_{\rm pwc}$ and the
maximum-efficiency protocol $\lambda_{\eta}$ (12) are in this case equal. We
used the same parameters as in Fig. 2 in the main text.

###  C.2 Constrained control

Let us now turn to the optimisation problem with constraints on the control:
$\lambda(t)\in[\lambda_{-},\lambda_{+}]$. As we know from the main text, in
this case the maximum-efficiency protocol is piecewise constant and maximum-
power protocol is piecewise linear. The numerical proof of the latter is given
in Fig. 2 in the main text. In Fig. 6, we show the numerical proof of the
former: the optimal performance of all the trial protocols in the figure is
the same as all of them collapse to the maximum-efficiency protocol
(parameters $b$ and $d$ for all the protocols equal to 0).

In Fig. 2 in the main text, we show that maximum powers of the protocols
$\lambda_{\eta}$ (12) and $\lambda_{\rm pwl}$ differ just by 1%. This suggests
that the maximum-efficiency protocol (12) can be used to estimate the overall
maximum power optimized with respect to $\lambda_{-}$. In Fig. 7, we provide
additional evidence supporting this conclusion. We show that for a broad range
of values of $T_{-}/T_{+}$ and $t_{-}/t_{+}$, the relative difference in
maximum power for the piecewice linear and (piecewise constant) maximum-
efficiency protocol is small. Furthermore, from panels (c)-(f) we conclude
that the optimal ratio $t_{-}/t_{+}$ is between 1 and 2, which is in agreement
with the results of Sec. B.2.2 [see Eq. (32) below]. Thus, for branch
durations that optimize output power, the relative difference $\delta P$ in
(a) is always below 12%, and it decreases with the temperature ratio.

##  Appendix D Optimal driving for systems close to equilibrium

In this section, we consider optimisation of a slowly driven heat engine based
on an overdamped Brownian particle trapped in the power-law potential
$H=\lambda(t)x^{n}/n$ with $n=2,4,\dots$. We use the temperature protocol from
Eq. (12) and impose fix values of the response $\sigma$ (or, in the slow
driving limit equivalently also the control $\lambda$) at the ends of the two
isotherms. Dynamics of the particle position is described by the Langevin
equation

$\dot{x}=-\mu\lambda(t)x^{n-1}+\sqrt{2D(t)}\xi(t),$ (37)

where $D(t)=\mu k_{\rm B}T(t)$ denotes the diffusion coefficient. From Eq.
(37) and its formal solution

$x(t)=-\mu\int dt\,\lambda(t)x^{n-1}(t)+\sqrt{2D(t)}\int dt\,\xi(t),$ (38)

we find that $\left<x(t)\xi(t)\right>=\sqrt{D/2}$ and thus

$\frac{d}{dt}\left<x^{2}(t)\right>=-2\mu\lambda(t)\left<x^{n}(t)\right>+2D.$
(39)

Let us now assume that the control parameters $\\{T(t),\lambda(t)\\}$ vary
slowly with respect to the relaxation time of the system, such that, during
the limit cycle, the system is always close to equilibrium, and solve this
equation up to the first order in $\dot{\lambda}(t)$. To this end, we consider
the ansatz $\left<x^{2}(t)\right>=\left<x^{2}(t)\right>_{0}$ and
$\left<x^{n}(t)\right>=\left<x^{n}(t)\right>_{0}+\left<x^{n}(t)\right>_{\dot{\lambda}}$,
where

$\left<x^{m}(t)\right>_{0}=\int_{-\infty}^{\infty}dxx^{m}\,\frac{\exp(-\mu\frac{\lambda
x^{n}}{nD})}{Z},$ (40)

with the partition function $Z=2[nD/\mu\lambda(t)]^{1/n}\Gamma(1+1/n)$, is the
value of the moment $\left<x^{m}(t)\right>$ corresponding to the infinitely
slow driving, and $\left<x^{n}(t)\right>_{\dot{\lambda}}$ is the correction
proportional to $\dot{\lambda}$. We find that

$\displaystyle\left<x^{n}(t)\right>_{0}$ $\displaystyle=$
$\displaystyle\frac{D(t)}{\mu\lambda(t)},$ (41)
$\displaystyle\left<x^{2}(t)\right>_{0}$ $\displaystyle=$
$\displaystyle\left[\frac{nD(t)}{\mu\lambda(t)}\right]^{2/n}\frac{\Gamma(3/n)}{\Gamma(1/n)},$
(42)

and

$\left<x^{n}(t)\right>_{\dot{\lambda}}=-\frac{1}{2\mu\lambda(t)}\frac{d}{dt}\left<x^{2}(t)\right>_{0}.$
(43)

We reiterate that this solution is valid only for protocols
$\\{T(t),\lambda(t)\\}$ which are changing slowly with respect to the
relaxation time of the system so that the system is during the whole cycle
close to equilibrium. However, as we know from the previous discussion, both
the piecewise constant maximum-efficiency protocol for constrained control and
the optimal protocol (35) for the constrained response contain
discontinuities, where $\\{T(t),\lambda(t)\\}$ changes abruptly. To be able to
use the slow driving approximation for the derivation of optimal cyclic
protocols, we thus need to additionally assume that during these jumps, the
system is not driven far from equilibrium. To this end, we assume that the
ratio $\lambda(t)/T(t)$ in the Boltzmann factor is during the jumps at the
ends of the isotherms constant. This additional assumption fixes the state of
the system $\sigma$ at the ends of the isotherms and thus the present
optimization scheme is only suitable for the optimization under the
constrained response. Let us now proceed with the optimization.

Work done on the system during the time interval $t_{\rm i}\leq t\leq t_{\rm
f}$ for the given Hamiltonian reads

$W=\frac{1}{n}\int_{t_{\rm i}}^{t_{\rm
f}}dt\,\dot{\lambda}(t)\left<x^{n}(t)\right>\equiv W(t_{\rm i},t_{\rm f}).$
(44)

Having fixed the state of the system at the ends of the isotherms, it is
enough to maximize the work during these branches. For an isothermal process,
the work Eq. (44) can be written as $W=\Delta F+W_{\rm irr}$, where the first
term, denoting the nonequilibrium free energy difference [24], comes from
$\left<x^{n}(t)\right>_{0}$, and the second term reads

$W_{\rm irr}=\frac{1}{n}\int_{t_{\rm i}}^{t_{\rm
f}}dt\,\dot{\lambda}(t)\left<x^{n}(t)\right>_{\dot{\lambda}}=\frac{1}{n^{2}\mu}\left(\frac{nD}{\mu}\right)^{2/n}\\\
\times\frac{\Gamma(3/n)}{\Gamma(1/n)}\int_{t_{\rm i}}^{t_{\rm
f}}dt\,\dot{\lambda}(t)^{2}\lambda(t)^{-2(1+n)/n}.$ (45)

As $\Delta F$ is fixed by the imposed boundary conditions on the state of the
system $\sigma$, to maximize the output work $-W$ means to minimize the
irreversible work $W_{\rm irr}$ as a functional of $\lambda(t)$. This leads to
the Euler-Lagrange equation

$\ddot{\lambda}(t)\lambda(t)-\frac{1+n}{n}\dot{\lambda}(t)^{2}=0,$ (46)

which has the general solution

$\lambda_{\rm slow}(t)=\frac{a}{(1+bt)^{n}}.$ (47)

We thus come to an interesting conclusion that the optimal slow protocol for
the constrained response scales with the same exponent as the potential. The
values of $a$ and $b$ can be expressed in terms of the boundary conditions for
$\lambda_{\rm slow}(t)$, i.e., $\lambda_{\rm slow}(t_{\rm
i})\equiv\lambda_{\rm i}$ and $\lambda_{\rm slow}(t_{\rm f})\equiv\lambda_{\rm
f}$. The optimal slow protocol (47) then reads

$\lambda_{\rm slow}(t)=\frac{\lambda(t_{\rm
i})}{\left[1+\left(\sqrt[n]{\frac{\lambda(t_{\rm i})}{\lambda(t_{\rm
f})}}-1\right)\frac{t-t_{\rm i}}{t_{\rm f}-t_{\rm i}}\right]^{n}}.$ (48)

And the corresponding irreversible work and input work are given by

$\displaystyle W_{\rm
irr}=\frac{\frac{\Gamma(3/n)}{\Gamma(1/n)}\left[\frac{nD}{\mu\lambda_{\rm
i}}\right]^{2/n}\left(\sqrt[n]{\frac{\lambda_{\rm i}}{\lambda_{\rm
f}}}-1\right)^{2}}{\mu(t_{\rm f}-t_{\rm i})},$ (49) $\displaystyle
W=\frac{\frac{\Gamma(3/n)}{\Gamma(1/n)}\left[\frac{nD}{\mu\lambda_{\rm
i}}\right]^{2/n}\left(\sqrt[n]{\frac{\lambda_{\rm i}}{\lambda_{\rm
f}}}-1\right)^{2}}{\mu(t_{\rm f}-t_{\rm
i})}-\frac{D}{n\mu}\ln\frac{\lambda_{\rm i}}{\lambda_{\rm f}}.$ (50)

These results are valid for the individual isothermal branches of the cycle.
Importantly, the obtained optimized values of the irreversible work are
correct up to the order $1/(t_{\rm f}-t_{\rm i})$, which is their exact
dependence on the process duration [24]. These results are thus exact even
though they were obtained from the approximate optimal protocol. According to
Refs. [24, 76], these irreversible works determine the optimal performance of
the engine under the constrains on $\sigma$, i.e., they give the maximum
output work $W_{\rm out}=-W(0,t_{+})-W(t_{+},t_{p})$ and efficiency
$\eta=W_{\rm out}/[T_{h}\Delta S-W_{\rm irr}(0,t_{+})]$ ($\Delta S$ is the
increase entropy of the system during the hot isotherm). Also this performance
is thus from the approximate analysis based on the slow driving obtained
exactly.