# Steps Minimize Dissipation in Rapidly Driven Stochastic Systems

Steven Blaber<EMAIL_ADDRESS>Dept. of Physics, Simon Fraser University,
Burnaby, British Columbia V5A 1S6, Canada  Miranda D. Louwerse  Dept. of
Chemistry, Simon Fraser University, Burnaby, British Columbia V5A 1S6, Canada
David A. Sivak<EMAIL_ADDRESS>Dept. of Physics, Simon Fraser University,
Burnaby, British Columbia V5A 1S6, Canada

###### Abstract

Micro- and nano-scale systems driven by rapid changes in control parameters
(control protocols) dissipate significant energy. In the fast-protocol limit,
we find that protocols that minimize dissipation at fixed duration are
universally given by a two-step process, jumping to and from a point that
balances jump size with fast relaxation. Jump protocols could be exploited by
molecular machines or thermodynamic computing to improve energetic efficiency,
and implemented in nonequilibrium free-energy estimation to improve accuracy.

The birth of thermodynamics as a modern science can be traced to Sadi Carnot’s
study of the design principles for energetically efficient heat engines
described in _Reflections on the Motive Power of Fire_ Carnot et al. (1960).
In classical thermodynamics, minimum-dissipation protocols are important in
the design of macroscopic heat engines describing, for example, adiabatic (no
heat loss) and quasistatic (infinitely slow) compression of gas by a piston.
Nearly 200 years later, the field of stochastic thermodynamics Jarzynski
(2011); Seifert (2012) similarly studies the design principles governing the
ability to dynamically vary control parameters and perform work at minimum
energetic cost (minimum dissipation), but now in micro- and nano-scale
fluctuating systems. Minimum-dissipation protocols inform our understanding of
the design principles of biological molecular machines Brown and Sivak (2017,
2019) and are of practical use to single-molecule experiments Tafoya et al.
(2019), free-energy estimation Shenfeld et al. (2009); Geiger and Dellago
(2010); Blaber and Sivak (2020a); Chodera et al. (2011); Minh and Chodera
(2009), and thermodynamic computing Conte et al. (2019); Proesmans et al.
(2020).

In contrast to macroscopic systems that are well described by averages of
thermodynamic quantities, the stochastic fluctuations in microscopic systems
are large relative to the mean and cannot be ignored. The work done on a
stochastic system by a control protocol is a fluctuating quantity that depends
on the entire protocol history, making it particularly difficult to optimize.
General optimization requires minimizing over all possible paths through
control-parameter space, which cannot be solved in general Schmiedl and
Seifert (2007). Despite advances relating optimal-control to optimal-transport
theory, even numerical optimizations are still limited to relatively simple
systems Aurell et al. (2011).

Although general solutions are intractable, we have gained considerable
insight into minimum-dissipation protocols by considering simple systems. For
example, Schmiedl and Seifert Schmiedl and Seifert (2007) showed that for a
Brownian particle diffusing in a harmonic potential with time-dependent
minimum or stiffness, minimum-dissipation protocols exhibit jump
discontinuities. It was posited that jumps in control parameters are a general
feature of minimum-dissipation protocols, and they have since been observed in
a number of different systems Gomez-Marin et al. (2008); Then and Engel
(2008); Esposito et al. (2010).

More general insight can be gained by approximating the mean work in relevant
limits. For slow protocols, linear-response theory yields an approximation for
the mean work, from which the approximate minimum-dissipation protocol can be
calculated Sivak and Crooks (2012). Despite its success, the linear-response
formalism relies on near-equilibrium approximations that break down for fast
protocols, miss key features of the exact minimum-dissipation protocol (e.g.,
jumps in control parameters), and for short duration can perform worse than
naive (constant-velocity) protocols Blaber and Sivak (2020b).

While minimum-dissipation protocols for slowly driven systems are relatively
well understood, comparatively little is known about rapidly driven systems.
In this work we focus on fast protocols and find a universal design principle:
the minimum-dissipation protocol consists of jumps at the beginning and end of
the protocol, spending the entire duration at the control-parameter value that
optimally balances the _initial force-relaxation rate_ (IFRR) (7b) with the
jump size (13). Our results are physically intuitive, apply to a wide range of
physical systems, and generalize easily to multidimensional control. To
illustrate this, we calculate the minimum-dissipation protocols in a diverse
set of systems described by Fokker-Planck or master-equation dynamics with
single- (Fig. 1) or multi-dimensional control (Fig. 4). Combining our results
with known minimum-dissipation protocols in the slow limit Sivak and Crooks
(2012), we demonstrate that a simple interpolation scheme produces protocols
that reduce dissipation at all speeds (Fig. 3).

_Derivation_.—Consider a thermodynamic system described by dynamics of the
form

$\displaystyle\frac{\partial p_{\Lambda}(x,t)}{\partial
t}=L[x,\boldsymbol{\lambda}(t)]\,p_{\Lambda}(x,t)\ ,$ (1)

where $p_{\Lambda}(x,t)$ is the probability distribution over microstates $x$
at time $t$ given the control protocol $\Lambda$, and
$L[x,\boldsymbol{\lambda}(t)]$ is the operator describing the system’s time
evolution. $L$ is the Liouvillian for Hamiltonian, drift/diffusion operator
for Fokker-Planck, and the transition-rate matrix for master-equation
dynamics. The system is in contact with a heat bath at temperature $T$ such
that the equilibrium probability distribution over microstates $x$ at fixed
control parameters $\boldsymbol{\lambda}$ is
$\pi(x|\boldsymbol{\lambda})\equiv\exp\\{\beta[F(\boldsymbol{\lambda})-U(x,\boldsymbol{\lambda})]\\}$,
for internal energy $U(x,\boldsymbol{\lambda})$ and free energy
$F(\boldsymbol{\lambda})\equiv-k_{\rm B}T\,\ln\sum_{x}\exp[-\beta
U(x,\boldsymbol{\lambda})]$, where $\beta\equiv(k_{\rm B}T)$ for Boltzmann’s
constant $k_{\rm B}$. The average excess work $W_{\rm ex}\equiv W-\Delta F$ by
an external agent changing control parameters $\boldsymbol{\lambda}$ according
to protocol $\Lambda$ is

$\displaystyle\langle W_{\rm ex}\rangle_{\Lambda}=-\int_{0}^{\Delta
t}\mathrm{d}t\,\frac{\mathrm{d}\boldsymbol{\lambda}^{\rm
T}}{\mathrm{d}t}\langle\delta{\bf f}(t)\rangle_{\Lambda}\ ,$ (2)

where a bold symbol denotes a column vector and superscript ${\rm T}$ the
transpose. ${\bf f}\equiv-\partial U/\partial\boldsymbol{\lambda}$ are the
forces conjugate to the control parameters, and $\delta{\bf f}\equiv{\bf
f}-\langle{\bf f}\rangle_{\rm eq}$ the deviations from the equilibrium
averages. Angle brackets $\langle\cdots\rangle_{\Lambda}$ denote a
nonequilibrium ensemble average over the control-parameter protocol $\Lambda$.
Here we hold fixed the initial ($\boldsymbol{\lambda}_{\rm i}$) and final
($\boldsymbol{\lambda}_{\rm f}$) control parameters, consistent with
nonequilibrium free-energy estimation Ritort et al. (2002); Gore et al.
(2003); Shenfeld et al. (2009); Geiger and Dellago (2010); Kim et al. (2012);
Blaber and Sivak (2020b); Schindler et al. (2020); Aldeghi et al. (2018); Kuhn
et al. (2017); Ciordia et al. (2016); Wang et al. (2015); Gapsys et al.
(2015); Chodera et al. (2011); Minh and Chodera (2009) but distinct from
optimizations that constrain the initial and final probability distributions
Proesmans et al. (2020); Zhang (2020).

If the total duration $\Delta t$ is short compared to the system’s natural
relaxation time $\tau$ (a fast protocol), expanding the probability
distribution in $\Delta t/\tau$ around an initial equilibrium distribution
gives

$\displaystyle p_{\Lambda}(x,s)=\pi(x|\boldsymbol{\lambda}_{\rm
i})+p^{1}_{\Lambda}(x,s)\frac{\Delta
t}{\tau}+\mathcal{O}\left[\left(\frac{\Delta t}{\tau}\right)^{2}\right]\ ,$
(3)

for $s\equiv t/\Delta t$ and first-order coefficient $p^{1}_{\Lambda}(x,s)$.
Plugging (3) up to $\mathcal{O}(\Delta t/\tau)$ into (1) gives

$\displaystyle\frac{\partial p^{1}_{\Lambda}(x,s)}{\partial
s}\approx\mathcal{L}[x,\boldsymbol{\lambda}(s)]\,\pi(x|\boldsymbol{\lambda}_{\rm
i})\ ,$ (4)

with $\mathcal{L}\equiv\tau L$ the dimensionless time-evolution operator.
Solving for $p^{1}_{\Lambda}(x,s)$ and substituting into (3) yields

$\displaystyle p_{\Lambda}(x,s)\approx\pi(x|\boldsymbol{\lambda}_{\rm
i})+\frac{\Delta
t}{\tau}\int_{0}^{s}\mathrm{d}s^{\prime}\mathcal{L}[x,\boldsymbol{\lambda}(s^{\prime})]\,\pi(x|\boldsymbol{\lambda}_{\rm
i})\ .$ (5)

Multiplying by conjugate forces ${\bf f}$, integrating over microstates $x$,
and changing the time variable back to $t$ gives

$\displaystyle\langle{\bf f}(t)\rangle_{\Lambda}\approx\langle{\bf
f}\rangle_{\boldsymbol{\lambda}_{\rm
i}}+\int_{0}^{t}\mathrm{d}t^{\prime}\,{\bf R}_{\boldsymbol{\lambda}_{\rm
i}}[\boldsymbol{\lambda}(t^{\prime})]\ ,$ (6)

for the _initial force-relaxation rate_ (IFRR)

$\displaystyle{\bf R}_{\boldsymbol{\lambda}_{\rm i}}[\boldsymbol{\lambda}(t)]$
$\displaystyle\equiv\int\mathrm{d}x\ {\bf
f}(x)\,L[x,\boldsymbol{\lambda}(t)]\,\pi(x|\boldsymbol{\lambda}_{\rm i})$ (7a)
$\displaystyle=\frac{\mathrm{d}\langle{\bf
f}\rangle_{\boldsymbol{\lambda}_{\rm
i}}}{\mathrm{d}t}\bigg{|}_{\boldsymbol{\lambda}(t)}\ ,$ (7b)

the rate of change of the conjugate forces at the current control-parameter
values (averaged over the initial equilibrium distribution).

Within this approximation, the average excess work is

$\displaystyle\langle W_{\rm ex}\rangle_{\Lambda}\approx\langle W_{\rm
ex}\rangle_{\boldsymbol{\lambda}_{\rm i}}-\int_{0}^{\Delta
t}\mathrm{d}t\,\frac{\mathrm{d}\boldsymbol{\lambda}^{\rm
T}}{\mathrm{d}t}\int_{0}^{t}\mathrm{d}t^{\prime}\,{\bf
R}_{\boldsymbol{\lambda}_{\rm i}}[\boldsymbol{\lambda}(t^{\prime})]\ .$ (8)

The first RHS term is the excess work during an instantaneous jump between the
initial and final control-parameter values, which equals the _relative
entropy_ $k_{\rm B}TD(\pi_{\rm i}||\pi_{\rm f})\equiv k_{\rm
B}T\int\mathrm{d}x~{}\pi_{\rm i}\ln[\pi_{\rm i}/\pi_{\rm f}]$ between the
initial $\pi_{\rm i}\equiv\pi(x|\boldsymbol{\lambda}_{\rm i})$ and final
$\pi_{\rm f}\equiv\pi(x|\boldsymbol{\lambda}_{\rm f})$ equilibrium probability
distributions Large and Sivak (2019). Integrating (8) by parts gives our main
theoretical result: for sufficiently short duration, the excess work is

$\displaystyle\langle W_{\rm ex}\rangle_{\Lambda}\approx k_{\rm B}TD(\pi_{\rm
i}||\pi_{\rm f})-\int_{0}^{\Delta t}\mathrm{d}t\,{\bf
R}_{\boldsymbol{\lambda}_{\rm i}}^{\rm
T}[\boldsymbol{\lambda}(t)]\,[\boldsymbol{\lambda}_{{\rm
f}}-\boldsymbol{\lambda}(t)]\ .$ (9)

The second RHS term is the first-order correction in $\Delta t$, an
approximation for the saved work $W_{\rm save}\equiv k_{\rm B}TD(\pi_{\rm
i}||\pi_{\rm f})-W_{\rm ex}$ compared to an instantaneous protocol.

_Initial Force-Relaxation Rate_.—The IFRR can be intuitively understood by
considering one-dimensional exponential relaxation. For a discrete jump from
initial control-parameter value $\lambda_{\rm i}$ to an intermediate value
$\lambda$, an exponentially relaxing mean conjugate force obeys

$\displaystyle\langle f(t)\rangle_{\Lambda}=\langle f\rangle_{\lambda_{\rm
i}}+\left(\langle f\rangle_{\lambda}-\langle f\rangle_{\lambda_{\rm
i}}\right)e^{-t/\tau(\lambda)}\ ,$ (10)

where $\tau(\lambda)$ is the relaxation time of the conjugate force. The IFRR
is the initial slope of the mean conjugate force as it relaxes towards
equilibrium (7b), which for exponential relaxation is

$\displaystyle{R}_{{\lambda}_{\rm i}}(\lambda)=\frac{\langle
f\rangle_{\lambda_{\rm i}}-\langle f\rangle_{\lambda}}{\tau(\lambda)}\ .$ (11)

Under simple exponential relaxation, $\tau(\lambda)$ is the same relaxation
time defined in Ref. Sivak and Crooks, 2012 for slow protocols, thereby
connecting short- and long-duration control.

For a small control-parameter jump $\lambda-\lambda_{\rm i}$, static linear-
response theory, $\langle f\rangle_{\lambda_{\rm i}}-\langle
f\rangle_{\lambda}\approx\beta(\lambda-\lambda_{\rm i})\langle\delta
f^{2}\rangle_{\lambda_{\rm i}}$, implies that the IFRR further simplifies to

$\displaystyle{R}_{{\lambda}_{\rm i}}(\lambda)\approx\frac{\langle\delta
f^{2}\rangle_{\lambda_{\rm i}}(\lambda-\lambda_{\rm i})}{\tau(\lambda)}\ .$
(12)

The relaxation rate is zero at the initial control-parameter value and
increases with larger control-parameter jumps which drive the system further
from equilibrium.

_Minimum-dissipation protocols_.—Equation (9) allows for relatively
straightforward optimization to determine the _short-time efficient protocol_
(STEP), satisfying the Euler-Lagrange equation

$\displaystyle\frac{\partial}{\partial\boldsymbol{\lambda}}\left[{\bf
R}_{\boldsymbol{\lambda}_{\rm i}}^{\rm
T}\left(\boldsymbol{\lambda}\right)\left(\boldsymbol{\lambda}_{{\rm
f}}-\boldsymbol{\lambda}^{\rm
STEP}\right)\right]\bigg{|}_{\boldsymbol{\lambda}^{\rm STEP}}={\bf
R}_{\boldsymbol{\lambda}_{\rm i}}\left(\boldsymbol{\lambda}^{\rm STEP}\right)\
.$ (13)

As an algebraic equation, the solution is a point in control-parameter space,
thus the STEP consists of two jumps: a jump at the start from its initial
value to the optimal value $\boldsymbol{\lambda}^{\rm STEP}$, and at the end
from the optimal value to the final value. The STEP is a jump protocol
provided the time-evolution operator $L$ is independent of time derivatives of
the control parameters. For Fokker-Planck dynamics this is satisfied if the
system is driven by a (generally time-dependent) scalar potential.

To illustrate the two-step minimum-dissipation protocol we have calculated the
STEP for diverse model systems (Fig. 1). In the translating- and breathing-
trap systems described by Fokker-Planck dynamics (Supplemental Material I SM
), the STEP jumps halfway between the two endpoints, consistent with the
results of Ref. Schmiedl and Seifert, 2007. The single-spin-flip and two-state
binding/unbinding reaction systems are described by master-equation dynamics
(Supplemental Material II and III SM ), with STEPs that jump to points that
are respectively larger and smaller than halfway between the initial and final
control-parameter values. Specific jump sizes for the STEP depend on the
functional form of the IFRR, but the minimum-dissipation protocol always
consists of jumps to and from an intermediate control-parameter value.

Figure 1: Short-time efficient protocols (STEPs) for a single spin in a time-
dependent magnetic field (blue dots), Brownian particle in a translating
harmonic potential (black dashed), Brownian particle in a harmonic potential
with time-dependent stiffness (black dashed), and a two-state
binding/unbinding reaction system with variable binding and unbinding rates
controlled by the chemical-potential difference (red dash-dotted).

The STEP jumps to the point in control-parameter space that maximizes the
_short-time power savings_

$\displaystyle P_{\rm save}^{\rm st}(\boldsymbol{\lambda})\equiv{\bf
R}_{\boldsymbol{\lambda}_{\rm i}}^{\rm
T}(\boldsymbol{\lambda})(\boldsymbol{\lambda}_{{\rm f}}-\boldsymbol{\lambda})$
(14)

due to relaxation at intermediate $\boldsymbol{\lambda}$. The STEP balances
fast relaxation rate ${\bf R}_{\boldsymbol{\lambda}_{\rm i}}$ with large final
jump $\boldsymbol{\lambda}_{{\rm f}}-\boldsymbol{\lambda}$. The STEP spends
the duration $\Delta t$ relaxing at $\boldsymbol{\lambda}^{\rm STEP}$, so for
short duration $P_{\rm save}^{\rm st}(\boldsymbol{\lambda}^{\rm{STEP}})\Delta
t$ is the work saved relative to an instantaneous protocol.

To demonstrate the energetics of the STEP, consider the thermodynamic cycle
consisting of tightening and loosening a harmonic trap (Fig. 2). For a
quasistatic (infinitely slow) protocol, work equals the free-energy
difference, which exactly cancels for a cyclic process. An instantaneous
protocol has an additional contribution, which for tightening (loosening) the
trap equals the relative entropy between the open (closed) and closed (open)
states. The relative entropy is dissipated as heat during equilibration
between tightening and loosening the trap (outer vertical arrows). The STEP
spends the duration relaxing at an intermediate control-parameter value,
resulting in saved work approximated by the area of the rectangle with width
given by the final jump size $\lambda_{{\rm f}}-\lambda^{\rm STEP}$ and height
by ${R}_{{\lambda}_{\rm i}}(\lambda^{\rm STEP})\Delta t$. To maximize the
saved work (rectangle area) the STEP optimally balances the IFRR (determining
the height) with final jump size (width).

Figure 2: Thermodynamic cycle in the force vs. control parameter plane for the
breathing harmonic trap driven by instantaneous (red dotted), STEP (green
dash-dotted), and quasistatic (black dashed) protocols. Arrows denote protocol
direction for transitions between open ($\lambda_{\rm i}$) and closed
($\lambda_{\rm f}$) states, shown schematically. The area under each curve
gives the average work done by the respective protocol, the area under the
quasistatic curve is the free-energy difference $\Delta F_{{\rm
i}\rightarrow{\rm f}}=F_{\rm f}-F_{\rm i}$, the area between the instantaneous
(dotted) and quasistatic (dashed) curves is the relative entropy (e.g.,
$\langle W\rangle_{\lambda_{\rm i}}\ -\Delta F_{{\rm i}\rightarrow{\rm
f}}=k_{\rm B}TD(\pi_{\rm i}||\pi_{\rm f})$), and the area between the STEP
(dash-dotted) and instantaneous (dotted) curves is the saved work $\langle
W_{\rm save}\rangle$ (shaded rectangles). Control-parameter endpoints satisfy
$\lambda_{\rm i}/\lambda_{\rm f}=1/2$, with duration $\Delta t/\tau=2/5$ for
fastest relaxation time $\tau=1/(2\lambda_{\rm f})$.

For a single control parameter, if the duration is sufficiently short the
_gain_ $G_{\rm save}\equiv{\langle W_{\rm save}\rangle_{\Lambda}^{\rm
des}}/{\langle W_{\rm save}\rangle_{\Lambda}^{\rm naive}}$ in saved work by
the STEP is

$\displaystyle G_{\rm save}^{\rm
STEP}\approx\frac{\max_{\boldsymbol{\lambda}}[P_{\rm save}^{\rm
st}(\lambda)]}{\overline{P_{\rm save}^{\rm st}(\lambda)}}\ ,$ (15)

where an overbar denotes the spatial average $\overline{P_{\rm save}^{\rm
st}(\lambda)}\equiv(\Delta\lambda)^{-1}\int_{\lambda_{\rm i}}^{\lambda_{\rm
f}}\mathrm{d}\lambda~{}P_{\rm save}^{\rm st}(\lambda)$, “naive” the constant-
velocity protocol, and “des” a designed protocol. The gain from a STEP is
greatest when the power savings $P_{\rm save}^{\rm st}(\lambda)$ is sharply
peaked.

_Interpolated protocols_.—In order to design a protocol that performs well for
any duration, we combine the STEP (optimal for fast protocols) with the
minimum-dissipation protocol from linear-response theory (optimal for slow
protocols). The linear-response protocol is continuous and when driven by a
single control parameter proceeds at velocity
$\mathrm{d}\lambda/\mathrm{d}t\propto[\zeta(\lambda)]^{-1/2}$, where
$\zeta(\lambda)$ is the generalized friction coefficient Sivak and Crooks
(2012). We assume the shape of the protocol from linear-response theory
remains unchanged (i.e.,
$\mathrm{d}\lambda/\mathrm{d}t\propto[\zeta(\lambda)]^{-1/2}$) but with
initial jump $(\lambda^{\rm STEP}-\lambda_{\rm i})/(1+\Delta t/\tau)^{\alpha}$
and final jump $(\lambda_{\rm f}-\lambda^{\rm STEP})/(1+\Delta
t/\tau)^{\alpha}$, where the constant $\alpha$ controls the crossover from
slow to fast approximations. For our simple systems we empirically find
$\alpha=1$ performs relatively well.

Figure 3 shows the improvement from designed protocols relative to naive
(constant-velocity) for the model system of a breathing harmonic trap. The
difference between naive and designed work (Fig. 3a) shows the expected
asymptotic behavior in the short- and long-duration limits: scaling as $\Delta
t$ (slope of one) for small $\Delta t/\tau$ and $(\Delta t)^{-1}$ (slope of
negative one) for large $\Delta t/\tau$. Both the fast and slow designed
protocols perform worse than naive (the difference is negative) for large and
small $\Delta t/\tau$, respectively. The fast-protocol approximation (9)
breaks down for long duration because the conjugate-force relaxation rate
decreases as the system approaches equilibrium at $\boldsymbol{\lambda}$,
whereas (9) assumes a constant relaxation time. However, the interpolated
protocol performs well for any duration, and the largest work saved compared
to naive is for intermediate duration. The gain $G_{\rm save}$ quantifies the
percent increase in saved work from designed protocols relative to naive,
where a gain greater than one indicates the designed does less work than the
naive. For this system, the largest gain in saved work occurs in the fast
limit (small $\Delta t/\tau$) for the STEP, interpolated, and exact optimal
protocols.

Figure 3: Benefit in the breathing harmonic trap from designed protocols
relative to the naive (constant-velocity) protocol, as a function of the
duration $\Delta t$ scaled by the fastest integral relaxation time $\tau$. The
different designed (“des”) protocols include the exact optimal Schmiedl and
Seifert (2007) (“opt”, solid black), linear-response optimized (“slow opt”,
dashed blue), STEP (“fast opt”, red dots), and interpolated optimal protocol
(“inter opt”, dash-dotted green). (a) Difference between the work done by the
naive (constant-velocity) and designed protocols. (b) Gain $G_{\rm
save}\equiv{\langle W_{\rm save}\rangle_{\Lambda}^{\rm des}}/{\langle W_{\rm
save}\rangle_{\Lambda}^{\rm naive}}$ in saved work. Solid red line in (b)
denotes the short-duration limit (15). Control-parameter endpoints satisfy
$\lambda_{\rm i}/\lambda_{\rm f}=16$, and the interpolated protocol uses
$\alpha=1$ and fastest integral relaxation time $\tau=1/(2\lambda_{\rm i})$
Blaber and Sivak (2020b).

_Multi-dimensional control_.—Optimization of multi-dimensional control
protocols has seen a recent surge in interest, primarily driven by possible
improvements to nonequilibrium free-energy estimates in fast-switching
simulations Chipot and Lelièvre (2011); Dellago and Hummer (2014). Previous
calculations of minimum-dissipation protocols which observed jumps were
limited to one-dimensional optimization. A significant advantage of the
present approximation is that it permits simple multidimensional control-
protocol optimization. Equation (13) implies that for multidimensional control
the STEP consists of jumps to and from the control-parameter point
$\boldsymbol{\lambda}^{\rm STEP}$.

To illustrate, we consider a nine-spin Ising model with frustrated boundary
conditions (Fig. 4) Rotskoff et al. (2017); Venturoli et al. (2009). We use a
2D control parameter $\boldsymbol{h}=(h_{\rm b},h_{\rm g})$ of magnetic fields
applied to the mid-edge spins (Fig. 4a) which initially hold the system in the
spin-down state and reverse during the protocol, driving the system to invert
its magnetization. Supplemental Material IV SM  gives model details.

Figure 4: a) Nine-spin ferromagnetic Ising model (internal black spins) with
fixed boundary conditions (external gray spins). The multi-dimensional control
parameter is two external magnetic fields, $h_{\rm b}$ (blue) applied to
horizontal-edge spins and $h_{\rm g}$ (green) applied to vertical-edge spins.
b) Short-time power savings (14) as function of control parameters ($h_{\rm
b}$,$h_{\rm g}$). Red line: naive protocol; red star:
$\boldsymbol{h}^{\rm{STEP}}$ (13). c) Work difference between designed and
naive protocols (dotted red), and its short-duration approximation (9) (solid
red). d) Gain $G_{\rm save}\equiv{\langle W_{\rm save}\rangle_{\Lambda}^{\rm
des}}/{\langle W_{\rm save}\rangle_{\Lambda}^{\rm naive}}$ in saved work for
multi-dimensional STEP relative to naive (dotted red), and its short-duration
limit (15) (solid red). Control-parameter endpoints are $\boldsymbol{h}_{\rm
i}=(-2,-2)$ and $\boldsymbol{h}_{\rm f}=(2,2)$, with duration $\Delta t$ and
fastest relaxation time $\tau=N/k_{0}$, for $N=9$ spins and single-spin flip
attempt rate $k_{0}$.

The power saving (14) vanishes at initial and final control-parameter values,
respectively corresponding to zero relaxation rate and zero final jump size
(Fig. 4b). By jumping past control-parameter regions with small power saving,
the STEP outperforms the naive protocol for short duration, as quantified by
the difference between naive and designed work (Fig. 4c) and the gain in saved
work (Fig. 4d). Indeed, for short duration the STEP more than doubles the work
saved by the naive protocol (i.e., has gain greater than two).

_Discussion_.—We have developed an approximation for work in the fast-protocol
limit (9) that permits straightforward optimization (13) simply from the
initial force-relaxation rate (IFRR), Eq. (7b). We find that jumps are a
universal feature of minimum-dissipation protocols in this fast limit, which
we illustrate with several model systems under single- (Fig. 1) or multi-
dimensional control (Fig. 4). Jumps minimize dissipation for fast protocols
because the relaxation rate is approximately constant, with no diminishing
returns from spending the entire duration at a single control-parameter value.
Therefore, the STEP jumps between the fixed control-parameter endpoints to
spend the entire duration at the control-parameter value that maximizes the
product of the IFRR and the subsequent jump size (14). This breaks down for
slow protocols since with sufficient time at a single control-parameter value,
the relaxation rate decreases over time; indeed, in the slow limit the
minimum-dissipation protocol is continuous Sivak and Crooks (2012). We combine
these two seemingly disparate limits with a simple interpolation scheme,
producing protocols that perform well for any duration (Fig. 3).

One important application of minimum-dissipation protocols is to free-energy
estimation, which aids the design of novel ligands for targeted protein
binding Schindler et al. (2020); Aldeghi et al. (2018); Kuhn et al. (2017);
Ciordia et al. (2016); Wang et al. (2015); Gapsys et al. (2015); Chodera et
al. (2011). Quite generally, the accuracy of free-energy estimates decreases
with increasing dissipation Ritort et al. (2002); Gore et al. (2003); Shenfeld
et al. (2009); Geiger and Dellago (2010); Kim et al. (2012); Blaber and Sivak
(2020b). Based on the results of Ref. Schmiedl and Seifert, 2007, jump
protocols have been used to reduce dissipation and improve free-energy
estimates Geiger and Dellago (2010), but previously it was unknown whether
jumps would always reduce dissipation in these more complex systems, and there
was no simple procedure to find the optimal jump size. The present formalism
demonstrates that jumps are a general feature and gives a simple optimization
procedure applicable to multidimensional control. This makes protocol
optimization tractable for a considerably expanded range of systems.

Although we focused on systems with known equations of motion, the IFRR (7b)
and short-time power savings (14) are easily estimated without detailed
dynamical information: the system only needs to be equilibrated at a single
control-parameter value; the protocol can be very short; the average converges
with few samples; and the STEP is found using standard optimization techniques
applied to (14). The STEP can thus be computed relatively inexpensively,
easing determination of minimum-dissipation protocols in rapidly driven
complex chemical and biological systems. This opens the door to improve the
energetic efficiency in thermodynamic computing Conte et al. (2019); Proesmans
et al. (2020) and the accuracy of nonequilibrium free-energy estimates in
simulations and single-molecule experiments Tafoya et al. (2019); Blaber and
Sivak (2020b); Shenfeld et al. (2009); Ritort et al. (2002).

This work is supported by an SFU Graduate Deans Entrance Scholarship (SB), an
NSERC CGS Doctoral fellowship (MDL), an NSERC Discovery Grant and Discovery
Accelerator Supplement (DAS), and a Tier-II Canada Research Chair (DAS), and
was enabled in part by support provided by WestGrid (www.westgrid.ca) and
Compute Canada Calcul Canada (www.computecanada.ca). The authors thank John
Bechhoefer, Jannik Ehrich, and Joseph Lucero (SFU Physics) for enlightening
feedback on the manuscript.

## References

  * Carnot et al. (1960) S. Carnot, E. Clapeyron, and R. Clausius, in _Reflections on the Motive Power of Fire: And Others Papers on the Second Law of Thermodynamics_ , edited by E. Mendoza (New York: Dover Publications, 1960).
  * Jarzynski (2011) C. Jarzynski, Annu. Rev. Condens. Matter Phys. 2, 329 (2011).
  * Seifert (2012) U. Seifert, Rep. Prog. Phys. 75, 126001 (2012).
  * Brown and Sivak (2017) A. I. Brown and D. A. Sivak, Physics in Canada 73 (2017).
  * Brown and Sivak (2019) A. I. Brown and D. A. Sivak, Chem. Rev. 120, 434 (2019).
  * Tafoya et al. (2019) S. Tafoya, S. J. Large, S. Liu, C. Bustamante, and D. A. Sivak, Proc. Natl. Acad. Sci. U.S.A 116, 5920 (2019).
  * Shenfeld et al. (2009) D. K. Shenfeld, H. Xu, M. P. Eastwood, R. O. Dror, and D. E. Shaw, Phys. Rev. E 80, 046705 (2009).
  * Geiger and Dellago (2010) P. Geiger and C. Dellago, Phys. Rev. E 81, 021127 (2010).
  * Blaber and Sivak (2020a) S. Blaber and D. A. Sivak, Phys. Rev. E 101, 022118 (2020a).
  * Chodera et al. (2011) J. D. Chodera, D. L. Mobley, M. R. Shirts, R. W. Dixon, K. Branson, and V. S. Pande, Curr. Opin. Struc. Biol. 21, 150 (2011).
  * Minh and Chodera (2009) D. D. Minh and J. D. Chodera, J. Chem. Phys. 131, 134110 (2009).
  * Conte et al. (2019) T. Conte, E. DeBenedictis, N. Ganesh, T. Hylton, J. P. Strachan, R. S. Williams, A. Alemi, L. Altenberg, G. Crooks, J. Crutchfield, et al., arXiv preprint arXiv:1911.01968 (2019).
  * Proesmans et al. (2020) K. Proesmans, J. Ehrich, and J. Bechhoefer, Phys. Rev. Lett. 125, 100602 (2020).
  * Schmiedl and Seifert (2007) T. Schmiedl and U. Seifert, Phys. Rev. Lett. 98, 108301 (2007).
  * Aurell et al. (2011) E. Aurell, C. Mejía-Monasterio, and P. Muratore-Ginanneschi, Phys. Rev. Lett. 106, 250601 (2011).
  * Gomez-Marin et al. (2008) A. Gomez-Marin, T. Schmiedl, and U. Seifert, J. Chem. Phys. 129, 024114 (2008).
  * Then and Engel (2008) H. Then and A. Engel, Phys. Rev. E 77, 041105 (2008).
  * Esposito et al. (2010) M. Esposito, R. Kawai, K. Lindenberg, and C. Van den Broeck, EPL 89, 20003 (2010).
  * Sivak and Crooks (2012) D. A. Sivak and G. E. Crooks, Phys. Rev. Lett. 108, 190602 (2012).
  * Blaber and Sivak (2020b) S. Blaber and D. A. Sivak, J. Chem. Phys. 153, 244119 (2020b).
  * Ritort et al. (2002) F. Ritort, C. Bustamante, and I. Tinoco, Proc. Natl. Acad. Sci. U.S.A 99, 13544 (2002).
  * Gore et al. (2003) J. Gore, F. Ritort, and C. Bustamante, Proc. Natl. Acad. Sci. U.S.A 100, 12564 (2003).
  * Kim et al. (2012) S. Kim, Y. W. Kim, P. Talkner, and J. Yi, Phys. Rev. E 86, 041130 (2012).
  * Schindler et al. (2020) C. E. Schindler, H. Baumann, A. Blum, D. Böse, H.-P. Buchstaller, L. Burgdorf, D. Cappel, E. Chekler, P. Czodrowski, D. Dorsch, et al., J. Chem. Inf. Model (2020).
  * Aldeghi et al. (2018) M. Aldeghi, V. Gapsys, and B. L. de Groot, ACS Cent. Sci. 4, 1708 (2018).
  * Kuhn et al. (2017) B. Kuhn, M. Tichý, L. Wang, S. Robinson, R. E. Martin, A. Kuglstatter, J. Benz, M. Giroud, T. Schirmeister, R. Abel, et al., J. Med. Chem 60, 2485 (2017).
  * Ciordia et al. (2016) M. Ciordia, L. Pérez-Benito, F. Delgado, A. A. Trabanco, and G. Tresadern, J. Chem. Inf. Model 56, 1856 (2016).
  * Wang et al. (2015) L. Wang, Y. Wu, Y. Deng, B. Kim, L. Pierce, G. Krilov, D. Lupyan, S. Robinson, M. K. Dahlgren, J. Greenwood, et al., J. Am. Chem. Soc. 137, 2695 (2015).
  * Gapsys et al. (2015) V. Gapsys, S. Michielssens, J. H. Peters, B. L. de Groot, and H. Leonov, in _Molecular Modeling of Proteins_ , edited by A. Kukol (Springer New York, New York, NY, 2015), pp. 173–209.
  * Zhang (2020) Y. Zhang, EPL 128, 30002 (2020).
  * Large and Sivak (2019) S. J. Large and D. A. Sivak, J. Stat. Mech.: Theory Exp. 2019, 083212 (2019).
  * (32) See Supplemental Material at [URL will be inserted by publisher] for Model Details.
  * Chipot and Lelièvre (2011) C. Chipot and T. Lelièvre, SIAM J. Appl. Math 71, 1673 (2011).
  * Dellago and Hummer (2014) C. Dellago and G. Hummer, Entropy 16, 41 (2014).
  * Rotskoff et al. (2017) G. M. Rotskoff, G. E. Crooks, and E. Vanden-Eijnden, Phys. Rev. E 95, 012148 (2017).
  * Venturoli et al. (2009) M. Venturoli, E. Vanden-Eijnden, and G. Ciccotti, J Math Chem 45, 188 (2009).
  * Glauber (1963) R. J. Glauber, J. Math. Phys. 4, 294 (1963).