Entropic equality for worst-case work at any protocol speed

Average values of quantities are not always typical values. In non-equilibrium nano and quantum systems this is often the case, with, for example, the work output of a protocol having a significant probability of deviating from the average. Hence, in these important systems, statements about averages have limited use when it comes to predicting what will happen in any given trial; the fluctuations need to be discussed explicitly. Two key relations concerning fluctuations in work, Crooks' theorem [1] and Jarzynski's Equality [2], have been studied extensively theoretically and experimentally. These theorems hold for any speed of changing the Hamiltonian, the thermalization can be partial or negligible during the protocol. Amongst other things the theorems can be used to determine free energies of equilibrium states from non-equilibrium experiments.

A recently developed alternative approach is single-shot statistical mechanics [3–12], inspired by single-shot information theory [13, 14]. The focus is on statements that are guaranteed to be true in every trial, rather than on average behaviors. For example, one can ask whether a process's work output is guaranteed to exceed some threshold value (such as an activation energy), or whether a process's work cost is guaranteed not to exceed some threshold value (beyond which the system may break from dissipating heat). These statements concern the worst-case work of a process. A key realization is that the optimal worst-case work is determined not by the von Neumann/Shannon entropy of the initial state, but rather the max entropy, which is the logarithm of the number of non-zero eigenvalues of the density matrix. Thus, which entropy one should use in statements about optimal work depends on which property of the work probability distribution one is interested in.

Single-shot statistical mechanics began with almost no a priori relation to fluctuation theorems. Promising links were made in [6, 15]. In [15] it was shown that Crooks' theorem can be used to make statements about worst-case work. This created the beginnings of a bridge between Crooks' theorem and results in single-shot papers. We here complete the bridge, showing that key expressions concerning optimal worst-case work from [3, 5, 6] follow from Crooks' theorem plus some extra thought. To our knowledge this is a new and unexpected result. We moreover generalize them by giving an equality for the worst-case work. The equality governs time-varying Hamiltonian protocols, including fast ones, assuming Markovian heat baths and a weak restriction on the strength of any initial quench. The initial state is taken to be diagonal in the energy eigenbasis, though not necessarily thermal, and energy coherences can arise during the protocol. Figure 1 illustrates the set-up. The equality has the form 'worst-case work = penalty—optimum,' with the penalty term guaranteed to be non-negative such that the optimum can be derived by setting the penalty to zero. We believe this concrete link to fluctuation theorems will give significant impetus to single-shot statistical mechanics, allowing it to harness results from the highly developed fluctuation theorem approach. As a demonstration we apply the result to an electron box experiment described in the fluctuation theorem formalism [16–18]. We begin by defining the general set-up.

Zoom In Zoom Out Reset image size

The physical scenario is depicted in figure 1. For concreteness, we shall use the following explicit model (and later discuss other possible models). A protocol will be a sequence of elementary changes: (i) changes of the Hamiltonian and (ii) thermalizations. We shall initially assume there is a finite number of such steps (but later show that the continuum limit is well-defined and corresponds to a master equation, at least in the discrete-classical case). The Hamiltonian is parameterized by ${\lambda }_{m}$, with m an integer that labels the step.

1. Hamiltonian changes map ${\lambda }_{m}$ to ${\lambda }_{m+1}$. We follow [19] in supposing there is an energy measurement in the instantaneous energy eigenbasis at the beginning and end of each Hamiltonian-changing step. In a given realization the system then evolves from $| {i}_{m},{\lambda }_{m}\rangle$ to $| {i}_{m}^{\prime },{\lambda }_{m+1}\rangle$, where i_m labels the energy eigenstate. This costs work given by the energy difference: ${w}_{m}=E(| {i}_{m},{\lambda }_{m+1}\rangle )-E(| {i}_{m}^{\prime },{\lambda }_{m}\rangle )$. An important special case is ${i}_{m}={i}_{m}^{\prime }$, which arises in the quasi-static (quantum adiabatic) limit, as well as if the energy eigenbasis is constant and only the energy eigenvalues change; this can be termed the discrete-classical case.

2. Thermalizations map ${i}_{m}^{\prime }$ to ${i}_{m+1}$, cost no work, and preserve the Hamiltonian: $| {i}_{m}^{\prime },{\lambda }_{m+1}\rangle \to | {i}_{m+1},{\lambda }_{m+1}\rangle$. For notational simplicity let us label this as $| i\rangle \to | j\rangle$ with energy ${E}_{i}\to {E}_{j}$. We do not assume that the system thermalizes fully but that the hopping probabilities respect thermal detailed balance: $\tfrac{p(| i\rangle \to \,| j\rangle )}{p(| j\rangle \to \,| i\rangle )}=\,{{\rm{e}}}^{-\beta ({E}_{j}-{E}_{i})}.$ The energy change ${E}_{j}-{E}_{i}$ from such a step is called heat, Q_m.

A trajectory is the time-sequence of energy eigenstates occupied: $| {i}_{0},{\lambda }_{0}\rangle \to | {i}_{0}^{\prime },{\lambda }_{1}\rangle \to | {i}_{1},{\lambda }_{1}\rangle$ $\to \,\ldots \,\to \,| {i}_{f-1},{\lambda }_{f-1}\rangle \to | {i}_{f-1}^{\prime },{\lambda }_{f}\rangle \to | {i}_{f},{\lambda }_{f}\rangle .$

The probability of a given trajectory is accordingly, assuming a Markovian heat bath,

$$\begin{eqnarray}p({\rm{traj}}) & = & p(| {i}_{0},{\lambda }_{0}\rangle )\displaystyle \prod _{m=0}^{f}p(| {i}_{m},{\lambda }_{m}\rangle \to | {i}_{m}^{\prime },{\lambda }_{m+1}\rangle )\\ & & \times p(| {i}_{m}^{\prime },{\lambda }_{m+1}\rangle \to | {i}_{m+1},{\lambda }_{m+1}\rangle ).\end{eqnarray} \tag{ 1 }$$

A trajectory's inverse is the reverse of the sequence. The inverse corresponds, in the discrete-classical case, to the Hamiltonian changes running in reverse, from ${\lambda }_{f}$ to ${\lambda }_{0}$, and to the same thermalizations as in the forward protocol, with the sequence exactly inverted. This process is termed the reverse process. Beyond the discrete-classical case, the Hamiltonian changes corresponding to the reverse process are defined such that $p(| {i}_{m}^{\prime },{\lambda }_{m+1}\rangle \to | {i}_{m},{\lambda }_{m}\rangle )=p(| {i}_{m},{\lambda }_{m}\rangle \to | {i}_{m}^{\prime },{\lambda }_{m+1}\rangle )$. Our results will hold under that condition. There are at least two ways of satisfying that condition: (i) Simply let the unitary of the corresponding elementary step in the reverse process be ${U}^{-1}$, where U is that of the forward process, (ii) apply a suitable 'time-reversal' operator Θ to all states and operators involved, as in [19]. The reverse trajectory is then the reverse sequence of the time-reversed energy eigenstates: ${\rm{\Theta }}| {i}_{f},{\lambda }_{f}\rangle ...{\rm{\Theta }}| {i}_{0},{\lambda }_{0}\rangle$, with the condition $p(| {i}_{m}^{\prime },{\lambda }_{m+1}\rangle \to | {i}_{m},{\lambda }_{m}\rangle )=p({\rm{\Theta }}| {i}_{m},{\lambda }_{m}\rangle \to {\rm{\Theta }}| {i}_{m}^{\prime },{\lambda }_{m+1}\rangle )$ being satisfied, as time reversal implies taking the complex conjugate of the states, in a preferred basis, and the transpose of the time-evolution in the same basis: $U\to {U}^{T}$. The condition is thus satisfied as $\langle b| U| a\rangle ={(\langle a| {U}^{\dagger }| b\rangle )}^{* }=\langle a{| }^{* }{U}^{T}| b{\rangle }^{* }$.

A given trajectory has some work cost $w={\sum }_{m}{w}_{m}$, in line with the definition of the Hamiltonian-changing steps. The inverse trajectory has work cost $-w$. A given protocol on a given initial state induces a probability distribution over trajectories, with an associated probability distribution over work p(w). The forward and reverse protocols give rise to ${p}_{\mathrm{fwd}}(w)$ and ${p}_{\mathrm{rev}}(-w)$ respectively.

If the initial density matrices of the forwards and reverse processes are both thermal, $\exp (-\beta H({\lambda }_{0}))/{Z}_{0}$ and $\exp (-\beta H({\lambda }_{f}))/{Z}_{f}$, Crooks' theorem holds [19]:

$$\begin{eqnarray}&&\displaystyle \frac{{p}_{\mathrm{fwd}}(w)}{{p}_{\mathrm{rev}}(-w)}=\displaystyle \frac{{Z}_{f}}{{Z}_{0}}\exp (\beta w).\end{eqnarray} \tag{ 2 }$$

(To derive it take the ratio of equation (1) and the corresponding reverse trajectory expression. Apply thermal detailed balance and the equality of reverse hopping probabilities for the Hamiltonian-changing steps. Sum over trajectories with the same w, and note that the reverse of a trajectory has the same work up to a minus sign [19].)

The central object of our interest is the worst-case work

$$\begin{eqnarray*}&&{w}^{0}:= \max \{w\,:p(w)\gt 0\},\end{eqnarray*}$$

also known as the guaranteed work [7]. In physical situations this will have a finite upper bound as no battery has infinite energy. The worst-case value may be realized by a very unlikely trajectory. It is then natural to consider the worst-case work of a subset of trajectories ${ \mathcal T }$:

$$\begin{eqnarray*}&&{w}_{{ \mathcal T }}^{0}:= \max \{w\,:p(w)\gt 0\ \mathrm{and}\ \mathrm{traj}\in { \mathcal T }\}.\end{eqnarray*}$$

The standard relative entropy is $D(\rho | | \sigma ):= -\mathrm{Tr}(\rho [\mathrm{log}\rho -\mathrm{log}\sigma ])$ [20], wherein logdenotes (in this paper) the natural logarithm, or ln. This D belongs to a class of relative entropies known as the Rényi relative entropies, parameterized by $\alpha \in {\mathbb{R}}$. We shall use two other members of that family: the (classical version of the) $\infty$-relative entropy ${D}_{\infty }(P| | Q):= {\sup }_{x}\mathrm{log}\left(\tfrac{{p}_{x}}{{q}_{x}}\right)$ and the 0-relative entropy ${D}_{0}(\rho | | \sigma ):= -\mathrm{logTr}({\pi }_{\rho }\sigma )$, wherein ${\pi }_{\rho }$ projects onto the support of ρ [21]. These are called one-shot relative entropies, as they arise naturally in one-shot (or single-shot) information theory [13, 14, 21].

It was shown in [15] that one can recover an expression for the worst-case work w⁰ from Crooks' theorem equation (2). Consider the equality of Crooks' theorem (for values of w such that ${p}_{\mathrm{fwd}}(w)\gt 0$) and select the value for w which maximizes the lhs (and thus the rhs) [15]:

$$\begin{eqnarray*}&&\max \displaystyle \frac{{p}_{\mathrm{fwd}}(w)}{{p}_{\mathrm{rev}}(-w)}=\max \displaystyle \frac{{Z}_{f}}{Z}{e}^{\beta w}.\end{eqnarray*}$$

The rhs is monotonic in w, so maximizing the rhs over the support of ${p}_{\mathrm{fwd}}(w)$ leads to the maximum w-value w⁰. Taking the logarithm and recalling the ${D}_{\infty }$ definition yields:

$$\begin{eqnarray}&&\beta {w}^{0}={D}_{\infty }({p}_{\mathrm{fwd}}(w)| | {p}_{\mathrm{rev}}(-w))-\mathrm{log}({Z}_{f}/Z).\end{eqnarray} \tag{ 3 }$$

Note that this derivation assumes Crook's theorem which does not in general hold for athermal initial states.

Consider an initial state ${\rho }_{0}$, and a protocol of thermalizations and Hamiltonian changes with initial and final Hamiltonians $H({\lambda }_{0})$ and $H({\lambda }_{f})$ respectively. This induces a work probability distribution p(w) and an associated w⁰. We shall derive an equality of the form ${w}^{0}\ =\,$penalty—optimum.

We consider initial states of form ${\rho }_{0}={\sum }_{i}{p}_{i}| {i}_{0},{\lambda }_{0}\rangle \langle {i}_{0},{\lambda }_{0}|$, i.e., diagonal in the energy eigenbasis though not necessarily thermal (energy coherence may still arise during the protocol). We take ${p}_{i}\ne 0$. This is because we wish to avoid divergences from dividing by p_i. (See [22] for an alternative approach to this divergence problem.)

To apply Crooks' theorem (equation (2)) here, even though the initial state is not assumed to be thermal, our approach is as follows. For example, if one has a degenerate two-level system, the thermal state is $\gamma =1/2| 0\rangle \langle 0| +1/2| 1\rangle \langle 1|$. If one instead had ${\rho }_{0}=2/3| 0\rangle \langle 0| +1/3| 1\rangle \langle 1|$, the worst-case work would be the same as for γ. This follows because the set of trajectories with non-zero probability is the same in both cases, as can be seen from equation (1) which gives the probability of a trajectory. Given a ${\rho }_{0}$, we will find a corresponding thermal state with the same worst-case work and apply Crooks' theorem to that.

An important practical consideration which makes this more subtle is that some p_i may be negligible and even arbitrarily close to 0. It is natural to exclude trajectories starting in those states when calculating the worst-case work. We therefore divide the initial energy eigenstates into two sets. One set is the one of interest: ${{ \mathcal E }}_{\mathrm{IN}}$. The set of the other eigenstates, we call ${{ \mathcal E }}_{\mathrm{OUT}}$, corresponding to those we shall exclude when calculating the worst-case work. The probability of being in ${{ \mathcal E }}_{\mathrm{OUT}}$ is given by

$$\begin{eqnarray*}&&p(\mathrm{OUT})=\displaystyle \sum _{| {i}_{0},{\lambda }_{0}\rangle \in {{ \mathcal E }}_{\mathrm{OUT}}}\mathrm{Tr}(| {i}_{0},{\lambda }_{0}\rangle \langle {i}_{0},{\lambda }_{0}| {\rho }_{0}).\end{eqnarray*}$$

We define ${{ \mathcal T }}_{\mathrm{IN}}$ as the set of possible ($p\gt 0$) trajectories beginning in ${{ \mathcal E }}_{\mathrm{IN}}$. Similarly, we define ${{ \mathcal T }}_{\mathrm{OUT}}$ as the set of possible trajectories beginning in ${{ \mathcal E }}_{\mathrm{OUT}}$. Recall that each trajectory corresponds to some work value. We call the worst-case work of ${{ \mathcal T }}_{\mathrm{IN}}$, ${w}_{\mathrm{IN}}^{0}$. This cannot be worse than the worst-case over all trajectories: ${w}_{\mathrm{IN}}^{0}\leqslant {w}^{0}$.

Let us design an associated thermal state that yields the same worst-case work as ${\rho }_{0}$: ${w}_{\mathrm{IN}}^{0}$. Later, we show that this is indeed the case, under an additional mild assumption. The associated thermal state has the same Hamiltonian as the system apart from the OUT levels. We define the Hamiltonian as $\widetilde{H}:= {\sum }_{\mathrm{IN}}{E}_{i}| i\rangle \langle i| +{\sum }_{\mathrm{OUT}}\widetilde{{E}_{i}}| i\rangle \langle i|$, changing the energies of the states in ${{ \mathcal E }}_{\mathrm{OUT}}$ to new ones, ${\widetilde{E}}_{i}$, such that ${p}_{i}=\exp (-\beta {\widetilde{E}}_{i})/\widetilde{Z}$, and leaving the other energy levels the same. The thermal state associated with that Hamiltonian is then

$$\begin{eqnarray*}&&\widetilde{\gamma }=\displaystyle \sum _{| {i}_{0},{\lambda }_{0}\rangle \in {{ \mathcal E }}_{\mathrm{IN}}}\displaystyle \frac{{{e}}^{-\beta {E}_{i}}}{\widetilde{Z}}| {i}_{0},{\lambda }_{0}\rangle \langle {i}_{0},{\lambda }_{0}| \,+\displaystyle \sum _{| {i}_{0},{\lambda }_{0}\rangle \in {{ \mathcal E }}_{\mathrm{OUT}}}{p}_{i}| {i}_{0},{\lambda }_{0}\rangle \langle {i}_{0},{\lambda }_{0}| .\end{eqnarray*}$$

The definition implies that

$$\begin{eqnarray}&&\widetilde{Z}=\displaystyle \frac{{\displaystyle \sum }_{| {i}_{0},{\lambda }_{0}\rangle \in {{ \mathcal E }}_{\mathrm{IN}}}{e}^{-\beta {E}_{i}}}{1-p({\rm{OUT}})}.\end{eqnarray} \tag{ 4 }$$

This partition function differs from that of the actual Hamiltonian $H({\lambda }_{0})$.

In this scenario with $\widetilde{\gamma }$ as the initial state and the ${{ \mathcal E }}_{\mathrm{OUT}}$ levels lifted, the protocol is the same as in the actual scenario, except that initially the energies of the states in ${{ \mathcal E }}_{\mathrm{OUT}}$ are lowered down to the levels of the actual Hamiltonian of interest. The worst-case work of this scenario is called ${\widetilde{w}}^{0}$. Under a mild additional restriction on protocols considered, roughly speaking that the worst-case work is bounded from below—as is the case for physically realizable protocols (see Methods), we then have

$$\begin{eqnarray}&&{\widetilde{w}}^{0}={w}_{\mathrm{IN}}^{0}.\end{eqnarray} \tag{ 5 }$$

To get ${\widetilde{w}}^{0}$ from Crooks' theorem (equation (2)) we shall make use of equation (3) from [15]. This applies in the scenario with $\widetilde{\gamma }$ as the initial state, as Crook's theorem holds in that scenario (see discussion around equation (2)), and thus

$$\begin{eqnarray}&&\beta {\widetilde{w}}^{0}={D}_{\infty }({\widetilde{p}}_{\mathrm{fwd}}(w)| | {\widetilde{p}}_{\mathrm{rev}}(-w))-\mathrm{log}({Z}_{f}/\widetilde{Z}).\end{eqnarray} \tag{ 6 }$$

Combining equation (6) and equation (5) we have

$$\begin{eqnarray}&&\beta {w}_{\mathrm{IN}}^{0}={D}_{\infty }({\widetilde{p}}_{\mathrm{fwd}}(w)| | {\widetilde{p}}_{\mathrm{rev}}(-w))-\mathrm{log}({Z}_{f}/\widetilde{Z}).\end{eqnarray} \tag{ 7 }$$

Thus the worst case work of the trajectories of interest ${w}_{\mathrm{IN}}^{0}$ is this equal to (kT times) a relative entropy minus (the logarithm) of the ratio of two partition functions, one of which encodes information about how many of the initial energy eigenstates have negligible occupation probability.

Note that equation (7) has the form

$$\begin{eqnarray*}&&\beta {w}^{0}=\ \mathrm{penalty}{\mbox{---}}\mathrm{optimum}.\end{eqnarray*}$$

The penalty is essentially given by the difference between the forward and reverse distributions, quantified by D_∞. The optimum one can hope for, with a given initial state and given initial and final Hamiltonian, is to set the penalty to 0 (as relative entropies are non-negative), which leaves $-\mathrm{log}({Z}_{f}/\widetilde{Z})$. This term is more negative the smaller the support of ρ is and the lower the final energies are relative to the initial ones. This optimum is achieved by a protocol in [7].

We now consider the optimum term in two important special cases where the single-shot entropy of the initial state emerges. To simplify the considerations we here set $p(\mathrm{OUT})\approx 0$, although our full expressions do not assume that to be the case. The 'epsilonics' are dealt with explicitly in the Methods. (i) Consider firstly the special case of $H({\lambda }_{0})=H({\lambda }_{f})$ which has been studied in the single-shot statistical mechanics literature. Then $\mathrm{log}({Z}_{f}/\widetilde{Z})=-\mathrm{log}{\sum }_{i\in {\rm{IN}}}{e}^{-\beta {E}_{i}}/{Z}_{f}$. This can be rewritten, informally, using the definition of D₀ in the technical introduction and noting that the sum is only over IN levels, as ${D}_{0}({\rho }_{0}^{({\rm{IN}})}| | \gamma )$, where ${\rho }_{0}^{({\rm{IN}})}$ is ${\rho }_{0}$ with the probability tail in OUT cut off—this is made general and precise in the discussion on smooth relative entropies in appendix B. Thus in this case the equality of equation (7) has the form

$$\begin{eqnarray*}&&\beta {w}^{0}={D}_{\infty }-{D}_{0}.\end{eqnarray*}$$

(Recall D_∞ concerns work distributions and D₀ states.) (ii) Further restricting the Hamiltonian such that $H({\lambda }_{0})=H({\lambda }_{f})=0$, we have ${D}_{0}({\rho }_{0}| | \gamma )=\mathrm{log}d-{S}_{\max }({\rho }_{0})$ (noting $\gamma ={\mathbb{1}}/d$ and recalling ${S}_{\max }(\rho ):= {S}_{0}(\rho ):= \mathrm{log}(\mathrm{rank}(\rho ))$. This recovers the known results from [3, 5, 6] that these are optimal in the respective cases. The message is that it is the max entropy S_max which determines the optimal worst-case work, rather than the von Neumann entropy. If one defines thermodynamic entropy in terms of optimally extractable worst-case work, it is the max entropy which should be used.

To make the connection to physics clear, we apply the results to a recent realization of a Szilard engine with an electron box [16–18]. A great advantage of using this trajectories model from the fluctuation theorem approach is that it allows the application of single-shot results to such experiments. We described the set-up in figure 2 and in the Methods section, we analyze what controls the penalty term D_∞ in this scenario. We also describe in the Methods how the penalty term, up to vanishing probabilistic error, goes to 0 in the isothermal quasistatic limit.

Zoom In Zoom Out Reset image size

As described in the trajectories section, these results also apply if the evolution includes unitaries that create energy coherences, including sudden changes of the energy eigenbasis such that a state that was orginally an energy eigenstate is now classified as a superposition of energy eigenstates. Such coherences are normally viewed as associated with entropy production and extra work costs. We note an interesting counter-example. Working within this trajectory model, suppose $H({\lambda }_{0})=0;$ ${\rho }_{0}=1/3| 0\rangle \langle 0| +2/3| 1\rangle \langle 1|$, and $H({\lambda }_{f})=\delta E| {i}_{f}\rangle \langle {i}_{f}|$. If the energy eigenstates stay the same throughout such that $| {i}_{f}\rangle =| 1\rangle$, the worst-case work is $\delta E$, and it has probability 2/3 (even if the shift is done quickly). If instead the Hamiltonian eigenstates change such that $| 0\rangle \to | +\rangle$, and $| 1\rangle \to | {i}_{f}\rangle =| -\rangle$ then the worst-case work is still $\delta E$, corresponding to outcome $| -\rangle$ of the final energy measurement. However the probability of this can be as low as 1/2 (if H is changed suddenly $p(| -\rangle )=\mathrm{Tr}({\rho }_{0}| -\rangle \langle -| )=1/2$). This shows that the probability of the worst case can actually be improved (lowered) by coherence due to a sudden change of the Hamiltonian. This improvement comes at the cost of randomizing the work distribution.

The derivation of the main result relies very little on the specifics of the trajectory model. It would e.g. also go through with the quantum-jump type model in [23, 24]. That model uses no intermediate projective measurements on the system but rather on the heat bath, as is natural in quantum optics.

We showed that in any protocol with a time-varying Hamiltonian and thermalizations, the worst-case work takes the form of 'penalty—optimum'. The model we used could be generalized in various ways, including non-Markovian baths and baths that decohere in other bases than the energy basis. It is also important to find more bounds for the penalty term in terms of controllable parameters. Finally we note that the results of [25] suggest that the bounds obtained here also apply to what is known as thermal operations in the context of resource theories—another interesting question.

We are grateful to comments on a draft by Dario Egloff. We acknowledge funding from the EPSRC (UK), the Templeton Foundation, the Leverhulme Trust, the Oxford Martin School, the National Research Foundation (Singapore), the EU collaborative project TherMiQ (Grant agreement No. 618074), the Ministry of Education (Singapore), NSF grant PHY-0803371, an IQIM Fellowship, and a Virginia Gilloon Fellowship, the NRF Grant 2015-003689 (Korea), a BK21 Plus Project (Korea), and the Gordon and Betty Moore Foundation.

As Crooks' theorem, which we use, is for initial thermal states, we designed a thermal state $\widetilde{\gamma }$, with the aim of creating a set-up which yields the same worst-case work as for the actual initial state, which is not necessarily thermal. For a given initial state ${\rho }_{0}=\sum {p}_{i}| i\rangle \langle i|$ and initial energy eigenvalues E_i, the associated thermal state is defined as $\widetilde{\gamma }={\sum }_{i}{e}^{-\beta \widetilde{{E}_{i}}}/\widetilde{Z}| i\rangle \langle i|$, where $\widetilde{{E}_{i}}={E}_{i}$ for $| i\rangle \in {{ \mathcal E }}_{\mathrm{IN}}$, but for $| i\rangle \in {{ \mathcal E }}_{\mathrm{OUT}}$, $\widetilde{{E}_{i}}$ is chosen such that ${e}^{-\beta \widetilde{{E}_{i}}}/\widetilde{Z}={p}_{i}$. Physically, this implies replacing the energy levels with small occupation probability p_i by much higher energy levels such that their thermal occupation probability is as small as p_i. The Hamiltonian associated with $\widetilde{\gamma }$ is accordingly $\widetilde{H}:= {\sum }_{\mathrm{IN}}{E}_{i}| i\rangle \langle i| +{\sum }_{\mathrm{OUT}}\widetilde{{E}_{i}}| i\rangle \langle i|$. The normalizing factor is $\widetilde{Z}={\sum }_{| i\rangle }{e}^{-\beta \widetilde{{E}_{i}}}$. These definitions imply that

$$\begin{eqnarray}&&\widetilde{Z}=\displaystyle \frac{{\displaystyle \sum }_{| i\rangle \in {{ \mathcal E }}_{\mathrm{IN}}}{e}^{-\beta {E}_{i}}}{1-p({\rm{OUT}})}.\end{eqnarray} \tag{ A1 }$$

Apart from the given actual protocol, we then also design a ∼-protocol such that it gives the same worst-case work in the case of $\widetilde{\gamma }$ as the initial state. We define the ∼-protocol as beginning with $\widetilde{H}$, then lowering the OUT levels back to E_i, i.e., setting $\widetilde{H}\to H$. After that it is the same as the actual protocol. We call the ∼-protocol applied to $\widetilde{\gamma }$ 'the ∼-scenario.'

In the ∼-scenario we similarly have ${\widetilde{{ \mathcal T }}}_{\mathrm{IN}}$ and ${\widetilde{{ \mathcal T }}}_{\mathrm{OUT}}$, and ${\widetilde{w}}_{\mathrm{IN}}^{0}$. The following holds:

$$\begin{eqnarray}&&{\widetilde{w}}_{\mathrm{IN}}^{0}={w}_{\mathrm{IN}}^{0},\end{eqnarray} \tag{ A2 }$$

i.e., the worst-case work is the same in the ∼-scenario as in the actual scenario, for the ${{ \mathcal T }}_{\mathrm{IN}}$ subset of trajectories. This is because the protocol is defined above such that the added initial step in the ∼-scenario only involves shifting the OUT levels (without any thermal hopping). The set of possible work values is the same in ${{ \mathcal T }}_{\mathrm{IN}}$ and ${\widetilde{{ \mathcal T }}}_{\mathrm{IN}}$.

We now make the following mild restriction on protocols allowed:

$$\begin{eqnarray}&&{\widetilde{w}}_{\mathrm{IN}}^{0}={\widetilde{w}}^{0}.\end{eqnarray} \tag{ A3 }$$

We say this is mild, because the trajectories ${\widetilde{{ \mathcal T }}}_{\mathrm{OUT}}$ have an extra work gain relative to their sister trajectories in ${{ \mathcal T }}_{\mathrm{OUT}}$ following from their initial lowering. This gain tends to infinity as $p({\rm{OUT}})\to 0$. Thus the ∼-protocol will not have worse possible work values than the actual protocol. The sort of protocols that are ruled out by the mild assumption are those where there is initially a dramatic quench without thermal hopping such that one of the OUT levels is raised very high. As long as it is not raised above the initial value of the level in the ∼-protocol the mild assumption is not violated however—the shifting down and up of the level would have a net work gain. This way the mild assumption does not rule out e.g. 2-level Szilard engine protocols where the second, likely empty, level is lifted very high initially and then lowered quasistatically. This example is studied in a physical system in appendix F below. Moreover, protocols involving thermalization at the start, including quasistatic protocols, or appropriately bounded initial quenches, are not ruled out.

Combining equations (A2) and (A3) gives the desired expression used in the main body:

$$\begin{eqnarray*}&&{w}_{\mathrm{IN}}^{0}={\widetilde{w}}^{0}.\end{eqnarray*}$$

To help illustrate the notation a simple example of the ∼-protocol and how it relates to the actual protocol is in figure A1.

Zoom In Zoom Out Reset image size

As noted in the main body, the optimum term reduces to a relative entropy in a special case. If $H({\lambda }_{0})=H({\lambda }_{f})$, $\mathrm{log}(Z/\widetilde{Z})=-\mathrm{log}{\sum }_{i\in {\rm{supp}}({\rho }_{0})}{e}^{-\beta {E}_{i}}/{Z}_{f}={D}_{0}({\rho }_{0}| | \gamma )$. Moreover if $H({\lambda }_{0})=H({\lambda }_{f})=0$, ${D}_{0}({\rho }_{0}| | \gamma )=\mathrm{log}d-{S}_{\max }({\rho }_{0})$ (noting $\gamma ={\mathbb{1}}/d$). This recovers the known results from [3, 5, 6] that these are optimal in the respective cases. If p(OUT) defined above is not necessarily zero, this optimal term depends on which levels are chosen to be in ${{ \mathcal E }}_{{\rm{OUT}}}$. If one chooses the best cut between IN and OUT, in the sense of minimizing $\widetilde{Z}$ and thus the worst-case work, the optimum one can hope for becomes in those cases $-{D}_{0}^{\epsilon }({\rho }_{0}| | \gamma ):= \min -{D}_{0}(\rho ^{\prime} | | \gamma )$ such that $d({\rho }_{0},\rho ^{\prime} )\leqslant \epsilon$ where d is the trace distance (this is called the smooth relative entropy). The interpretation is that the optimal worst-case work one can hope for allowing for an error tolerance of $\epsilon =p({\rm{OUT}})$ is ${{kTD}}_{0}^{\epsilon }({\rho }_{0}| | \gamma )$, consistent with [3, 5, 6].

Certain protocols studied in single-shot statistical mechanics give a pseudo-deterministic work output, i.e. the work distribution is highly concentrated around some value. (It is standard to say that a certain amount of work A is (δ, )-deterministic if one will have $W\in [A\pm \delta ]$ except with probability $\leqslant \epsilon$.) For example in [3] one may compress all the randomness onto some bits and extract work from the others with essentially deterministic work output. In [5, 6] and and [26] the expressions given concern the optimal pseudo-deterministic work, optimized over protocols for some given boundary conditions. This is in contrast to this paper and e.g. [7] which make bounds on the optimal worst-case work. We note that from the definitions one sees that bounds on worst-case work are also bounds on deterministic work but not vice versa. This is because demanding that the work cost is bounded from above is a weaker restriction than demanding that it is bounded from both above and below.

There can actually be (sets of) trajectories which are unlikely even if the initial state of the trajectory is likely, as the hopping probability may be low. For example, if one lifts one level towards a very high value while thermalizing, there is one trajectory corresponding to staying in the rising level throughout, which would be the trajectory that gives the worst-case work. However, if such a trajectory is very unlikely, one would wish to ignore the trajectory when stating the worst-case work. In this section, we show a way to ignore such unlikely trajectories, by not only cutting off part of the initial state as previously, but also cutting off a part of the work distribution. This strategy gives a different penalty term—lower, in general—in the equality for the worst-case work.

D.1. Proof overview

D.2. The set of trajectories of interest

We wish to ignore unlikely trajectories. We identify a set of trajectories of interest, defined as excluding trajectories of two types:

1. ${\rho }_{0}$-tail trajectories: These are those which are called ${{ \mathcal T }}_{{\rm{OUT}}}$ above, i.e., trajectories which start in ${{ \mathcal E }}_{{\rm{OUT}}}$. We now call them ${\rho }_{0}$-tail trajectories as using OUT risks generating confusion because of the second type of cut we shall make on the set of trajectories.

2. Work-tail trajectories: We also ignore trajectories associated with the worst work values, if those values are sufficiently improbable. This ignoring amounts to cutting off the worst-case tail of the work probability distribution. To simplify the proof, we define this tail in terms of the work probability distribution of the fictional thermal state $\widetilde{\gamma }$. By 'the work-tail,' we mean the set of trajectories associated with the following work values w: if the initial state is $\widetilde{\gamma }$, there is an associated work probability distribution ${\widetilde{p}}_{\mathrm{fwd}}(w)$ for the given protocol, and an associated worst-case work ${\widetilde{w}}^{\epsilon }$. The work tail trajectories are by definition those with work cost $w\gt {\widetilde{w}}^{\epsilon }$. Since the actual initial state ${\rho }_{0}$ may differ from $\widetilde{\gamma }$, the probability that some trajectory begins in the work tail does not necessarily equal .

These sets are depicted in figure D1. We shall call the worst-case work in the set of interest ${w}_{{\rm{IN}},{\rm{IN}}}^{0}$.

Zoom In Zoom Out Reset image size

D.3. The worst-case work in that set

We now derive the worst-case work in the set of trajectories of interest: We maximize the work cost w over that set of interest. We shall, for the first part, draw inspiration from an argument, in [15], concerning scenarios governed by Crooks' Theorem. Take the initial state of the forwards process to be ${\rho }_{0}=\widetilde{\gamma };$ and the initial state of the reverse process as $\gamma ={{e}}^{-\beta H({\lambda }_{f})}/{Z}_{f}$.

Maximize Crooks' theorem over the support of ${\tilde{p}}_{\mathrm{fwd}}(w)$ [15]:

$$\begin{eqnarray*}&&\max \displaystyle \frac{{\widetilde{p}}_{\mathrm{fwd}}(w)}{{p}_{\mathrm{rev}}(-w)}=\max \displaystyle \frac{Z}{\widetilde{Z}}{e}^{\beta w}.\end{eqnarray*}$$

The rhs is monotonic in w, so maximizing the rhs over the support of ${\tilde{p}}_{\mathrm{fwd}}(w)$ leads to the maximum w-value w⁰. Taking the logarithm and recalling the ${D}_{\infty }$ definition yields [15],

$$\begin{eqnarray*}&&\beta {w}^{0}={D}_{\infty }({\widetilde{p}}_{\mathrm{fwd}}(w)| | {p}_{\mathrm{rev}}(-w))-\mathrm{log}\left(\displaystyle \frac{Z}{\widetilde{Z}}\right).\end{eqnarray*}$$

Now, we cut off the work tail by defining a cut-off probability distribution ${\widetilde{p}}_{\mathrm{fwd}}^{\epsilon }(w):= 0$, if $w\geqslant {\widetilde{w}}^{\epsilon }$ and $\tfrac{{\widetilde{p}}_{\mathrm{fwd}}(w)}{1-\epsilon }$, otherwise, wherein ${\widetilde{w}}^{\epsilon }$ denotes the work guaranteed up to probability if $\widetilde{\gamma }$ is the initial state. (Dividing by ($1-\epsilon$) normalizes the distribution.) For work values outside the work tail, Crooks' theorem can be reformulated as

$$\begin{eqnarray*}&&\displaystyle \frac{{\widetilde{p}}_{\mathrm{fwd}}^{\epsilon }(w)}{{p}_{\mathrm{rev}}(-w)}(1-\epsilon )=\displaystyle \frac{Z}{\widetilde{Z}}{e}^{\beta w}.\end{eqnarray*}$$

Since the rhs is monotonic,

$$\begin{eqnarray*}&&\max \displaystyle \frac{{\widetilde{p}}_{\mathrm{fwd}}^{\epsilon }(w)}{{\widetilde{p}}_{\mathrm{rev}}(-w)}(1-\epsilon )=\left(\displaystyle \frac{Z}{\widetilde{Z}}\right){e}^{\beta {\widetilde{w}}^{\epsilon }},\end{eqnarray*}$$

wherein the maximization is over the support of ${\widetilde{p}}_{\mathrm{fwd}}^{\epsilon }$. Taking the logarithm and rearranging yields

$$\begin{eqnarray*}&&\beta {\widetilde{w}}^{\epsilon }={D}_{\infty }({\widetilde{p}}_{\mathrm{fwd}}^{\epsilon }(w)\,| | {\widetilde{p}}_{\mathrm{rev}}(-w))+\mathrm{log}(1-\epsilon )-\mathrm{log}\left(\displaystyle \frac{Z}{\widetilde{Z}}\right).\end{eqnarray*}$$

The lhs is the worst-case work in the set of trajectories of interest.

D.4. Probability that a trajectory is in the set of interest

The trajectories of interest are effectively the possible trajectories. To make precise what is meant by 'effective,' we bound the probability that any particular trajectory lies outside that set.

Consider a trajectory followed by a system initialized to ${\rho }_{0}$. The probability that the trajectory lies outside the set of interest is bounded by $p({\rho }_{0}-\mathrm{tail})+p(\mathrm{work}-\mathrm{tail})$, as shown in figure D1. $p({\rho }_{0}-\mathrm{tail})$, defined via ${\rho }_{0}$ and the choice of effective support, is specified by input parameters. $p(\mathrm{work}-\mathrm{tail})$ denotes the probability that the trajectory is in the set associated with a worse work cost than ${\widetilde{w}}^{\epsilon }$ (the work guaranteed up to probability not to be exceeded, if the initial state is $\widetilde{\gamma }$). $p(\mathrm{work}-\mathrm{tail})$ does not necessarily equal for an arbitrary ${\rho }_{0}$. As $p(\mathrm{work}-\mathrm{tail})$ is not an input parameter, we wish to bound it with input parameters.

Let us drop the subscript 'fwd' and refer simply to p(w). The weight $p(w\gt x)$ in the actual work tail associated with ${\rho }_{0}$ cannot differ arbitrarily from the weight $\widetilde{p}(w\gt x)$ in the work tail associated with $\widetilde{\gamma }$:

$$\begin{eqnarray*}&&| p(w\gt x)-\widetilde{p}(w\gt x)| \leqslant d(p(w),\widetilde{p}(w)).\end{eqnarray*}$$

This bound follows from the definition of the variational distance d, which equals the trace distance between diagonal states⁸ .

The variational distance d is contractive under stochastic matrices, because the trace distance is contractive under completely positive trace-preserving maps. The work distribution is the result of a stochastic matrix's acting on the probability distribution over initial energy eigenstates. Let us now in this paragraph for convenience use Dirac notation for classical probability vectors, representing a probability distribution p(w) as $\langle w| p\rangle$. The work distribution comes from the stochastic matrix ${\sum }_{j}| {p}_{j}\rangle \langle j|$ mapping a state $| {\rho }_{0}\rangle$ to a work distribution, wherein j labels projectors onto $H({\lambda }_{0})$ eigenstates, $| {p}_{j}\rangle$ labels the work distribution when starting with an initial state $| j\rangle$ (i.e., ${p}_{j}(w)=\langle w| {p}_{j}\rangle$), and $| {\rho }_{0}\rangle ={\sum }_{j}{q}_{j}| j\rangle$. For example, if there are two possible eigenstates, we can write $| {\rho }_{0}\rangle ={q}_{1}| 1\rangle +{q}_{2}| 2\rangle ={({q}_{1}{q}_{2})}^{T}$, and the resulting work distribution $p(w)=(\langle w| {p}_{1}\rangle \langle 1| +\langle w| {p}_{2}\rangle \langle 2| )| {\rho }_{0}\rangle ={q}_{1}{p}_{1}(w)+{q}_{2}{p}_{2}(w)$.

Thus,

$$\begin{eqnarray*}&&| p(w\gt x)-\widetilde{p}(w\gt x)| \leqslant d(p(w),\widetilde{p}(w))\leqslant d({\rho }_{0},\widetilde{\gamma })\,\forall x.\end{eqnarray*}$$

For some $x=x^{\prime}$, by definition, $\widetilde{p}(w\gt x^{\prime} )=\widetilde{p}(\mathrm{work}-\mathrm{tail})=\epsilon$, and $p(\mathrm{work}-\mathrm{tail}):= p(w\gt x^{\prime} )$. Thus

$$\begin{eqnarray*}&&p(\mathrm{work}-\mathrm{tail})\leqslant d({\rho }_{0},\widetilde{\gamma })+\epsilon .\end{eqnarray*}$$

D.5. Main result, also cutting work tail

We conclude that the worst-case work from the trajectories of interest, ${{w}}_{{\rm{IN}},{\rm{IN}}}^{0}$ respects

$$\begin{eqnarray}&&\beta {w}_{{\rm{IN}},{\rm{IN}}}^{0}={D}_{\infty }({\widetilde{p}}_{\mathrm{fwd}}^{\epsilon }(w)\,| | {p}_{\mathrm{rev}}(-w))+\mathrm{log}(1-\epsilon )-\mathrm{log}Z/\widetilde{Z}.\end{eqnarray} \tag{ D1 }$$

The probability that the trajectory is not in the set of interest is upper-bounded by $p({\rho }_{0}-{\rm{tail}})+p({\rm{work}}-{\rm{tail}})\leqslant p({\rho }_{0}-{\rm{tail}})+d({\rho }_{0},\widetilde{\gamma })+\epsilon$.

We have mainly focused on the discrete-time protocol. Experimental realizations of thermodynamic protocols are often described by a continuous master equation. Here, we show that the discrete protocol leads to a master equation in the continuum model and vice versa. In this section we restrict ourselves to scenarios without energy coherences, i.e., the discrete-classical case.

E.1. From discrete to continuous

We consider a discrete sequence of times, ${t}_{m}={t}_{0}+m\ {\rm{d}}{t}$ ($m=0,1,2\ldots$), and the sequence ${\lambda }_{m}\equiv \lambda ({t}_{m})$ of values of the external parameter. As the waiting time decreases (${\rm{d}}{t}\to 0$), the transition probability $p(| i,\lambda (t),t\rangle \to | j,\lambda (t+{\rm{d}}{t}),t+{\rm{d}}{t}\rangle )$ due to thermalization should vanish. To first order, it behaves as

$$\begin{eqnarray}&&p(| i,\lambda (t),t\rangle \to | j,\lambda (t+{\rm{d}}{t}),t+{\rm{d}}{t}\rangle )\approx {\delta }_{{ij}}+{{\rm{\Gamma }}}_{i\to j}(t){\rm{d}}{t}+{ \mathcal O }({{\rm{d}}{t}}^{2}).\end{eqnarray} \tag{ E1 }$$

The transition rate ${{\rm{\Gamma }}}_{i\to j}(t)$ is a possibly complicated function of instantaneous energy levels $E(| i,\lambda (t),t\rangle )$. However, the transition rates inherit the condition

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Gamma }}}_{i\to j}(t)}{{{\rm{\Gamma }}}_{j\to i}(t)}={e}^{-\beta [E(| j,\lambda (t),t\rangle )-E(| i,\lambda (t),t\rangle )]}\end{eqnarray} \tag{ E2 }$$

from detailed balance and the condition

$$\begin{eqnarray}&&\displaystyle \sum _{j}{{\rm{\Gamma }}}_{i\to j}(t)=0\end{eqnarray} \tag{ E3 }$$

from probability conservation. The occupation probability is

$$\begin{eqnarray*}p(| j,\lambda (t+{\rm{d}}{t}),t+{\rm{d}}{t}\rangle ) & = & \displaystyle \sum _{i}p(i,| \lambda (t),t\rangle )p(| i,\lambda (t),t\rangle \to | j,\lambda (t+{\rm{d}}{t}),t+{\rm{d}}{t}\rangle )\\ & \approx & p(| j,\lambda (t),t\rangle )+\displaystyle \sum _{i}p(i,| \lambda (t),t\rangle ){{\rm{\Gamma }}}_{i\to j}(t){\rm{d}}{t}-\displaystyle \sum _{i}p(j,| \lambda (t),t\rangle ){{\rm{\Gamma }}}_{j\to i}(t){\rm{d}}{t}.\end{eqnarray*}$$

If the occupation probability is a smooth function of time, the master equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}p(| j,\lambda (t+{\rm{d}}{t}),t+{\rm{d}}{t}\rangle )=\displaystyle \sum _{i}p(i,| \lambda (t),t\rangle ){{\rm{\Gamma }}}_{i\to j}(t)-\displaystyle \sum _{i}p(j,| \lambda (t),t\rangle ){{\rm{\Gamma }}}_{j\to i}(t)\end{eqnarray} \tag{ E4 }$$

follows. The equivalence is further illustrated in appendix F in the example of an electron box.

E.2. From continuous to discrete

Going in the other direction, we now show explicitly how the discrete-time model can be derived from a physical master equation. Consider a two-level system that has a state $| 0\rangle$, kept at zero energy, and a state $| 1\rangle$ whose energy ${\hslash }\omega (t)$ changes. The Hamiltonian is $H(t)={\hslash }\omega (t)| 1\rangle \langle 1|$, and the system interacts with a temperature-T heat bath. In [27], a master equation for the density matrix $\rho (t)$ was derived for a such system. In the present case, the master equation is

$$\begin{eqnarray}&&\dot{\rho }(t)=-{\rm{i}}[H(t),\rho (t)]+{ \mathcal L }(t)\rho (t)\end{eqnarray} \tag{ E5 }$$

$$\begin{eqnarray}&&{ \mathcal L }(t)\rho ={\rm{\Gamma }}{\rm{d}}(\omega (t))([{n}_{{\rm{th}}}(\omega (t))+1]\{[{\sigma }_{-},\rho (t){\sigma }_{+}]+{\rm{h}}.{\rm{c}}.\}+{n}_{{\rm{th}}}(\omega (t))\{[{\sigma }_{+},\rho (t){\sigma }_{-}]+{\rm{h}}.{\rm{c}}.\}).\end{eqnarray} \tag{ E6 }$$

The heat bath, modeled as as set of harmonic oscillators, has a thermal occupation number ${n}_{\mathrm{th}}{(\omega )=({e}^{\beta {\hslash }\omega }-1)}^{-1}$ that depends on time because the upper level shifts. ${\rm{d}}(\omega )$ is the dimensionless heat-bath density of states; Γ denotes a rate assumed to be constant; ${\sigma }_{-}=| 0\rangle \langle 1|$ denotes the usual lowering operator; and ${\sigma }_{+}={\sigma }_{-}^{\dagger }$. Equation (E5) has the form of the usual Lindblad master equation, but the Lindblad operator depends on time. The dependence arises only from the level spacing's time dependence. The Hamiltonian part contains the Lamb shift.

In the derivation of equation (E5) one assumes, as usual, weak coupling to the heat bath, the Markovian approximation, and the rotating-wave approximation. One also assumes that the adiabatic approximation holds, i.e., the system always remains in its time-local energy eigenstates when the interaction with the heat bath is ignored. This condition is always fullfilled under the assumption of vanishing energy coherences at all times that we made in this section. Indeed, the part of (E5) pertaining to the diagonal elements of $\rho (t)$ can be derived without the adiabatic assumption [28].

We now consider discrete times ${t}_{n}:= n{\rm{\Delta }}t$, $n=0,\ldots ,N$, with $\omega (t)$ constant during the time intervals ${\rm{\Delta }}t$, ${\omega }_{n}:= \omega ({t}_{n})$. Restricting ourselves to changes of the Hamiltonian that only involve its spectrum, H(t) and ${ \mathcal L }(t)$ are constant during a given time interval.

Consider first the Hamiltonian changes. Heisenberg's equation of motion for the system-and-bath composite implies that $\dot{\rho }(t)$ has a finite jump when the Hamiltonian has a finite jump. Therefore, $\rho (t)$ is continuous when the Hamiltonian has a finite jump. Hence for finite Hamiltonian changes during a time $\delta t$, the system-and-bath composite's density matrix is unchanged in the limit as δt → 0. Hence the system's reduced density matrix is unchanged during the instantaneous shift of energy levels. As for the relaxation process, the initial thermal state is described in terms of occupation probabilities p_n for the nth level. The evolution during the relaxation process is given by $p(t)={e}^{{Tt}}p(0)$, where T is a matrix that connects the diagonal matrix elements of ρ in the master equation (E5): ${\dot{\rho }}_{{nn}}={\sum }_{m}{T}_{{nm}}{\rho }_{{mm}}$. The transition rates T_nm inherit detailed balance from the rates appearing in the master equation, i.e., ${T}_{{ij}}={e}^{-\beta ({\epsilon }_{i}-{\epsilon }_{j})}{T}_{{ji}}$. Detailed balance holds for each power T^k of T, ${({T}^{k})}_{{ij}}={e}^{-\beta ({\epsilon }_{i}-{\epsilon }_{j})}{({T}^{k})}_{{ji}}$, for all $k\in {\mathbb{N}}$. Therefore, by the power-series expansion of e^Tt, detailed balance holds also for e^Tt. We thus have derived, from a physical model of a system that is coupled to a heat bath and whose energy levels are piecewise-constant, the discrete-time model considered in the paper.

To illustrate this let us consider a two level system: expressing $\rho (t)={p}_{0}(t)| 0\rangle \langle 0| +[1-{p}_{0}(t)]| 1\rangle \langle 1|$, we obtain a differential equation for ${p}_{0}(t)$,

$$\begin{eqnarray}&&{\dot{p}}_{0}(t)+g(\omega (t)){p}_{0}(t)=\displaystyle \frac{g(\omega (t))}{2}+{\rm{\Gamma }}{\rm{d}}(\omega (t)),\end{eqnarray} \tag{ E7 }$$

wherein $g(\omega (t)):= 2{\rm{\Gamma }}{\rm{d}}(\omega (t))[2{n}_{\mathrm{th}}(\omega (t))+1]$. This equation has the general solution

$$\begin{eqnarray}&&{p}_{0}(t)=\left({p}_{0}(0)-\displaystyle \frac{1}{2}+{\rm{\Gamma }}{\int }_{0}^{t}{{\rm{d}}{t}}_{1}{\rm{d}}(\omega ({t}_{1}))G({t}_{1})\right)/G(t)+\displaystyle \frac{1}{2},\end{eqnarray} \tag{ E8 }$$

wherein $G(t):= {e}^{{\int }_{0}^{t}g(\omega ({t}_{1})){{\rm{d}}{t}}_{1}}$. The integrals in equation (E8) can be calculated analytically:

$$\begin{eqnarray}&&{p}_{0}(t)={p}_{0}(0){{\rm{e}}}^{-2{\rm{d}}(\omega )\coth (\tfrac{\beta {\hslash }\omega }{2}){\rm{\Gamma }}t}+{p}_{0,\mathrm{th}}\left[1-{{\rm{e}}}^{-2{\rm{d}}(\omega )\coth (\tfrac{\beta {\hslash }\omega }{2}){\rm{\Gamma }}t}\right],\end{eqnarray} \tag{ E9 }$$

wherein ${p}_{0,\mathrm{th}}:= 1/({e}^{-\beta {\hslash }\omega }+1)$ denotes the ground state's thermal occupation. For large times, the memory of the initial state is lost, and the system relaxes towards thermal equilibrium. From equation (E9), we obtain the transition probabilities during relaxation over the time interval between $n{\rm{\Delta }}t$ and $(n+1){\rm{\Delta }}t$: $p(| {0}_{n},{\omega }_{n}\rangle \to | {0}_{n+1},{\omega }_{n}\rangle )={p}_{0}({\rm{\Delta }}t){| }_{{p}_{0}(0)=1}$, $p(| {0}_{n},{\omega }_{n}\rangle \to | {1}_{n+1},{\omega }_{n}\rangle )=1-p(| {0}_{n},{\omega }_{n}\rangle \to | {0}_{n+1},{\omega }_{n}\rangle )$, $p(| {1}_{n},{\omega }_{n}\rangle \to | {0}_{n+1},{\omega }_{n}\rangle )={p}_{0}({\rm{\Delta }}t){| }_{{p}_{0}(0)=0}$, and $p(| {1}_{n},{\omega }_{n}\rangle \to | {1}_{n+1},{\omega }_{n}\rangle )=1-p(| {1}_{n},{\omega }_{n}\rangle \to | {0}_{n+1},{\omega }_{n}\rangle )$. These transition probabilities obey detailed balance. As they remain unchanged by the inclusion of an instantaneous Hamiltonian change at the end of each time interval, $p(| {i}_{n},{\omega }_{n}\rangle \to | {j}_{n+1},{\omega }_{n}\rangle )=p(| {i}_{n},{\omega }_{n}\rangle \to | {j}_{n+1},{\omega }_{n+1}\rangle )$ for $i,j\in \{0,1\}$.

To demonstrate the physical relevance of our results, we take a realistic example, the electron box [16–18, 29, 30], and apply our results to it. We first derive a time-local master equation for the level-occupation probabilities in appendix F.1. As shown in appendix E, the master equation is equivalent to the discrete-time trajectory model discussed in the main text. The work distribution functions are analyzed numerically in appendix F.2 and analytically in appendix F.3. Finally, we upper-bound the penalty term ${D}_{\infty }$, which reveals the direct physical relevance of our results.

F.1. Theoretical model and its justification

We consider the type of system in [16–18]. Following a semiclassical theory (known as 'the orthodox theory') such as in [31], we derive a master equation and illustrate the work fluctuations. A more complete quantum description is possible [29, 30]. Yet the semiclassical approach is useful for interpreting and identifying work and heat, which are often ambiguous.

The system (figure F1) consists of a large metallic electrode R that serves as a charge reservoir, a small metallic island (or quantum dot) D, and a gate electrode. The island D is coupled only capacitively to the gate electrode but couples to the reservoir R capacitively and via tunneling. The Hamiltonian has four parts: $H={H}_{R}+{H}_{D}+{H}_{C}+{H}_{T}.$ The first two terms,

$$\begin{eqnarray}&&{H}_{R}=\displaystyle \sum _{k}{\varepsilon }_{k}{c}_{k}^{\dagger }{c}_{k}\quad \mathrm{and}\quad {H}_{D}=\displaystyle \sum _{q}{\varepsilon }_{q}{d}_{q}^{\dagger }{d}_{q},\end{eqnarray} \tag{ F1 }$$

describe the non-interacting parts of the electrode R and the island D. Here, ${c}_{k}^{\dagger }$ (${d}_{q}^{\dagger }$) creates an electron with momentum ${\hslash }k$ (${\hslash }q$) and energy ${\varepsilon }_{k}$ (${\varepsilon }_{q}$). The single-particle dispersions ${\varepsilon }_{k}$ and ${\varepsilon }_{q}$ form continua of energy levels. H_C is responsible for the electron–electron interaction on the island. We adopt a capacitive model as recounted below. Finally, the tunneling of electrons between R and D is described by

$$\begin{eqnarray}&&{H}_{T}=\eta \displaystyle \sum _{{kq}}{c}_{k}^{\dagger }{d}_{q}+{\rm{h}}.{\rm{c}}.,\end{eqnarray} \tag{ F2 }$$

wherein η is the tunneling amplitude. η is assumed not to depend on momenta (or on energy), as in common metals that have wide conduction bands.

Zoom In Zoom Out Reset image size

The effective semiclassical model describes equilibrium: suppose that an electron tunnels between the island D and the reservoir R. The tunneling jolts the system out of equilibrium, but the system equilibrates quickly: Coulomb repulsions redistribute the electrons throughout the circuit. After the redistribution, the junction carries the equilibrium charge Q_J, and the gate carries the equilibrium charge Q_g. These charges are regarded as being 'on' the island D, due to the island's capacitive couplings to the reservoir R and to the gate. (The island carries also excess electrons, discussed below.) The electrons continue to repel each other, imbuing the system with the equilibrium Coulomb energy

$$\begin{eqnarray}&&{H}_{C}=\displaystyle \frac{{Q}_{J}^{2}}{2{C}_{J}}+\displaystyle \frac{{Q}_{g}^{2}}{2{C}_{g}},\end{eqnarray} \tag{ F3 }$$

wherein C_J and C_g denote the junction and gate capacitances. One can find that

$$\begin{eqnarray}&&{Q}_{J}=C({V}_{g}-{Ne}/{C}_{g}),\end{eqnarray} \tag{ F4a }$$

$$\begin{eqnarray}&&{Q}_{g}=C({V}_{g}+{Ne}/{C}_{J}),\end{eqnarray} \tag{ F4b }$$

wherein $C:= {C}_{g}{C}_{J}/({C}_{g}+{C}_{J})$ is the system's effective capacitance. $N={\sum }_{k}{d}_{k}^{\dagger }{d}_{k}$ denotes the number of excess electrons, relative to the charge-neutral state, on the island D. When N = 0, the island has zero net charge. H_C can thus be rewritten as

$$\begin{eqnarray}&&{H}_{C}={E}_{C}{N}^{2}+\displaystyle \frac{1}{2}{{CV}}_{g}^{2},\end{eqnarray} \tag{ F5 }$$

wherein ${E}_{C}:= {e}^{2}/2({C}_{g}+{C}_{J})$ is the single-electron charging energy, one of the largest energy scales of the system.

We are primarily interested in the macroscopic variable N but not in the microscopic degrees of freedom c_k and d_q, whose dynamics is typically much faster. One can thus integrate out c_k and d_q to get the effective Hamiltonian expressed only in terms of N. In the semiclassical approach, this can be achieved by considering the energy that an electron gains by tunneling.

Suppose that an electron tunnels into the island D from the reservoir R. This will change the charge ${Q}_{J}\to {Q}_{J}-e$ and the excess number of electrons $N\to N+1$. This new charge configuration, right after the tunneling, is redistributed quickly to a new equilibrium configuration,

$$\begin{eqnarray}&&{Q}_{J}^{\prime }=C[{V}_{g}-(N+1)e/{C}_{g}],\end{eqnarray} \tag{ F6a }$$

$$\begin{eqnarray}&&{Q}_{g}^{\prime }=C[{V}_{g}+(N+1)e/{C}_{J}],\end{eqnarray} \tag{ F6b }$$

by the gate voltage source. The voltage source has moved the amount

$$\begin{eqnarray}&&{\rm{\Delta }}Q:= {Q}_{J}^{\prime }-({Q}_{J}-e)={{eC}}_{g}/({C}_{g}+{C}_{J})\end{eqnarray} \tag{ F7 }$$

of charge through the transmission line from the junction interface to the gate capacitor by doing the amount

$$\begin{eqnarray}&&W={V}_{g}{\rm{\Delta }}Q={{eV}}_{g}{C}_{g}/({C}_{g}+{C}_{J})\end{eqnarray} \tag{ F8 }$$

of work on the system. Therefore, the electron's overall energy gain ${\rm{\Delta }}E$ equals the work W minus the change in the electrostatic energy:

$$\begin{eqnarray}&&{\rm{\Delta }}E={E}_{C}[2{C}_{g}{V}_{g}/e-(2N+1)].\end{eqnarray} \tag{ F9 }$$

As this energy gain comes from the transition $N\to N+1$, the effective Hamiltonian for the macroscopic variable N can be regarded as

$$\begin{eqnarray}&&{H}_{\mathrm{eff}}={E}_{C}({N}^{2}-2{{NN}}_{g}),\end{eqnarray} \tag{ F10 }$$

wherein ${N}_{g}:= {C}_{g}{V}_{g}/e$. Recall that the second term comes from the work done on the system by the voltage source.

The microscopic degrees of freedom removed from the effective macroscopic model cause N to fluctuate randomly. The transition $N\to N\pm 1$ is associated with tunneling of an electron into or from the island. Hence the transition rate follows from Fermi's Golden Rule:

$$\begin{eqnarray}&&{\rm{\Gamma }}({\rm{\Delta }}E)\approx \displaystyle \frac{2\pi | \eta {| }^{2}{\rho }_{R}{\rho }_{D}}{{\hslash }}\displaystyle \frac{{\rm{\Delta }}E}{{e}^{\beta {\rm{\Delta }}E}+1},\end{eqnarray} \tag{ F11 }$$

wherein ${\rho }_{R}$ and ${\rho }_{D}$ are the densities of states of R and D, respectively, and

$$\begin{eqnarray}&&{\rm{\Delta }}E={H}_{\mathrm{eff}}(N\pm 1)-{H}_{\mathrm{eff}}(N).\end{eqnarray} \tag{ F12 }$$

Finally, at sufficiently low temperatures ($\beta {E}_{C}\gg 1$), higher energy levels play no role. Considering the two lowest levels N = 0 and N = 1 suffices for ${N}_{g}\in [0,1]$ ⁹ . Let p₀ denote the probability that N = 0, and let p₁ denote the probability that N = 1. With equations (F10) and (F11), this two-level approximation leads to the master equation

$$\begin{eqnarray}&&{\dot{p}}_{0}=-{{\rm{\Gamma }}}_{+}{p}_{0}+{{\rm{\Gamma }}}_{-}{p}_{1},\end{eqnarray} \tag{ F13a }$$

$$\begin{eqnarray}&&{\dot{p}}_{1}=-{{\rm{\Gamma }}}_{-}{p}_{1}+{{\rm{\Gamma }}}_{+}{p}_{0}.\end{eqnarray} \tag{ F13b }$$

The transition rates are [29, 30]

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\pm }(t):= {\rm{\Gamma }}(\pm \epsilon (t))\quad \mathrm{and}\quad {\rm{\Gamma }}(\epsilon ):= \displaystyle \frac{{{\rm{\Gamma }}}_{0}\epsilon (t)/{\varepsilon }_{c}}{{e}^{\beta \epsilon (t)}-1}.\end{eqnarray} \tag{ F14 }$$

Here, ${\varepsilon }_{c}$ is the bath's high-frequency cutoff (i.e., ${\hslash }/{\varepsilon }_{c}$ is the correlation time), and ${{\rm{\Gamma }}}_{0}$ is a constant that characterizes the strength of the coupling to the bath. ${{\rm{\Gamma }}}_{0}/{\varepsilon }_{c}$ is related to the material properties by ${{\rm{\Gamma }}}_{0}/{\varepsilon }_{c}=2\pi | \eta {| }^{2}{\rho }_{R}{\rho }_{D}/{\hslash }$. Note that the transition rates satisfy the detailed-valance relation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Gamma }}(+\epsilon )}{{\rm{\Gamma }}(-\epsilon )}={e}^{-\beta \epsilon }.\end{eqnarray} \tag{ F15 }$$

The time-local master equation (F13) is equivalent to the discrete-time trajectory model (see appendix E). Therefore, the electron box is a realistic prototype system to which our results can apply.

F.2. Monte Carlo simulation of the electron box

We performed a Monte Carlo simulation of an erasure protocol in the electron box set-up. Our simulation discretizes the protocol into time steps $\delta t$ small enough to justify the linear approximation that the population of level i evolves from time step t to $t+\delta t$ according to ${p}_{i}(t+\delta t)={p}_{i}(t)+\delta t{\dot{p}}_{i}(t)$. Using equations (F13), we can write a stochastic matrix acting on the probabilities:

$$\begin{eqnarray}&&\left[\begin{array}{c}{p}_{0}(t+\delta t)\\ {p}_{1}(t+\delta t)\end{array}\right]=\left[\begin{array}{cc}1-{{\rm{\Gamma }}}_{+}\delta t & {{\rm{\Gamma }}}_{-}\delta t\\ \delta t{{\rm{\Gamma }}}_{+} & 1-{{\rm{\Gamma }}}_{-}\delta t\end{array}\right]\left[\begin{array}{c}{p}_{0}(t)\\ {p}_{1}(t)\end{array}\right].\end{eqnarray} \tag{ F16 }$$

For a two-level system which does not build up quantum coherences, a stochastic thermalizing matrix (which by its definition evolves all states towards the Gibbs state) has only one degree of freedom remaining once the Gibbs state has been chosen: the thermalization speed. All models of two-level thermalizations for a given Gibbs state are equivalent. We pick the conceptually straightforward partial swap: with some probability p_sw, the system's current state is exchanged with the Gibbs state. With probability $1-{p}_{\mathrm{sw}}$, the state remains unchanged: ${M}_{\mathrm{swap}}=(1-{p}_{\mathrm{sw}}){\mathbb{1}}+{p}_{\mathrm{sw}}| \mathrm{Gibbs}\rangle \langle \mathrm{ones}|$, where $| \mathrm{ones}\rangle$ means the vector of 1's. For a Gibbs state associated with an energy-level splitting ,

$$\begin{eqnarray}{M}_{\mathrm{swap}}=\left[\begin{array}{cc}1-\displaystyle \frac{{p}_{\mathrm{sw}}\exp (-\beta \epsilon )}{1+\exp (-\beta \epsilon )} & \displaystyle \frac{{p}_{\mathrm{sw}}}{1+\exp (-\beta \epsilon )}\\ \displaystyle \frac{{p}_{\mathrm{sw}}\exp (-\beta \epsilon )}{1+\exp (-\beta \epsilon )} & 1-{p}_{\mathrm{sw}}\displaystyle \frac{1-\exp (-\beta \epsilon )}{1+\exp (-\beta \epsilon )}\end{array}\right].\end{eqnarray} \tag{ F17 }$$

Equating equation (F17) with the matrix in equation (F16), we can find the partial-swap probability in terms of the electron box's physical parameters:

$$\begin{eqnarray}&&{p}_{\mathrm{sw}}(t)=\displaystyle \frac{{{\rm{\Gamma }}}_{0}\delta t}{{\varepsilon }_{c}}\epsilon (t)\coth [\beta \epsilon (t)/2].\end{eqnarray} \tag{ F18 }$$

The swap probability ${p}_{\mathrm{sw}}(t)$ and the energy level splitting $\epsilon (t)$ appear as functions of time, as the swap probability changes as the protocol evolves. The probability changes only as a function of an external parameter, the splitting (as opposed to e.g., the current state). Hence Crooks' theorem is still applicable to thermalizations of this type.

Figure F2 depicts our Monte Carlo simulation. We randomly generate trajectories by picking a random initial microstate according to the initial-state probability distribution. Then, we evolve the system by small steps, testing at each step if a swap should occur (with probability p_sw). If a swap occurs, we replace the state with a new microstate randomly chosen according to the Gibbs state associated with the current Hamiltonian. By recording which microstate is occupied when the energy level is raised, we calculate the work cost associated with a particular trajectory. Repeated runs of the simulation allow us to build up a work distribution, to which the results in this paper apply.

Zoom In Zoom Out Reset image size

F.3. Analytic expression for the work distribution

The work distribution function for an electron box can also be obtained explicitly from the master equation.

Consider an arbitrary work protocol that runs from t = 0 to $t=\tau$. The gap is tuned as a function $\epsilon (t)$. The trajectory $\sigma (t)\in \{0,1\}$ of the system is piece-wise constant, jumping discontinuously from one energy level to another at some random instants t_j ($j=1,2,\ldots$). Therefore, the trajectory is specified uniquely by the initial condition ${\sigma }_{0}$, the number J of jumps, and the corresponding instants t_j ($j=1,2,\ldots ,J$). The probability distribution function for the trajectory is

$$\begin{eqnarray}&&{P}_{J}({t}_{1},\ldots ,{t}_{J};{\sigma }_{0})=\displaystyle \prod _{j=1}^{J}{\rm{\Gamma }}({(-1)}^{{\sigma }_{0}+j+1}\epsilon ({t}_{j}))\exp [-{S}_{J}({t}_{1},\ldots ,{t}_{J};{\sigma }_{0})],\end{eqnarray} \tag{ F19 }$$

where the effective action associated with a given trajectory has been defined as

$$\begin{eqnarray}&&{S}_{J}({t}_{1},\ldots ,{t}_{J};{\sigma }_{0})=\displaystyle \sum _{j=1}^{J+1}{\int }_{{t}_{j-1}}^{{t}_{j}}{\rm{d}}{s}\,{\rm{\Gamma }}({(-1)}^{{\sigma }_{0}+j+1}\epsilon (s))\end{eqnarray} \tag{ F20 }$$

and ${t}_{0}=0$ and ${t}_{J+1}=\tau$ are implied. Checking the normalization is straightforward:

$$\begin{eqnarray}&&{P}_{0}({\sigma }_{0})+\displaystyle \sum _{J=1}^{\infty }\displaystyle \prod _{j=1}^{J}{\int }_{{t}_{j-1}}^{\tau }{{\rm{d}}{t}}_{j}\,{P}_{J}({t}_{1},\ldots ,{t}_{J};{\sigma }_{0})=1,\end{eqnarray} \tag{ F21 }$$

wherein ${t}_{0}=0$ is again implied.

The work is done only while the system is in the state $\sigma =1$. Hence the contribution to the work along the trajectory is

$$\begin{eqnarray}&&{W}_{J}({t}_{1},\ldots ,{t}_{J};{\sigma }_{0})=\displaystyle \sum _{j=1}^{J}{(-1)}^{{\sigma }_{0}+j}\epsilon ({t}_{j})+({\sigma }_{0}+J\mathrm{mod}2){\epsilon }_{f}-{\sigma }_{0}{\epsilon }_{0}.\end{eqnarray} \tag{ F22 }$$

The work distribution function along a trajectory with J jumps is

$$\begin{eqnarray}&&{P}_{J}(W;{\sigma }_{0})=\displaystyle \prod _{j=1}^{J}{\int }_{{t}_{j-1}}^{\tau }{{\rm{d}}{t}}_{j}\,{P}_{J}({t}_{1},\ldots ,{t}_{J};{\sigma }_{0})\delta (W-{W}_{J}({t}_{1},\ldots ,{t}_{J};{\sigma }_{0})).\end{eqnarray} \tag{ F23 }$$

The total work distribution function can be written in a series

$$\begin{eqnarray}&&P(W)={p}_{0}{e}^{-{S}_{0}(0)}\delta (W)+{p}_{1}{e}^{-{S}_{0}(1)}\delta (W-{W}_{c})+\displaystyle \sum _{J=1}^{\infty }\displaystyle \sum _{{\sigma }_{0}}{p}_{{\sigma }_{0}}{P}_{J}(W;{\sigma }_{0}).\end{eqnarray} \tag{ F24 }$$

P_J(W) has a factor of ${({{\rm{\Gamma }}}_{0}^{2}{e}^{-\beta \epsilon })}^{J}$. At low temperatures, P_J is rapidly suppressed as J increases.

The expression (F24) for the work distribution is essentially a perturbative expansion in ${{\rm{\Gamma }}}_{0}^{2}$ and converges very quickly for small ${{\rm{\Gamma }}}_{0}$. For large ${{\rm{\Gamma }}}_{0}$, however, it becomes impractical to use it for actual calculation because of its slow convergence. Therefore, it will be useful to devise a more general method. We examine the characteristic function $Z(\xi )=\langle {e}^{\xi W}\rangle$ of the work distribution function P(W). We first consider the characteristic function ${Z}_{\sigma }(\xi )=\langle {e}^{\xi W}{\rangle }_{\sigma }$ conditioned on all trajectories' starting from a definite initial state ${\sigma }_{0}$. Regarded as a function of the operation time τ, ${Z}_{\sigma }(\xi ;\tau )$ satisfies the master equation [32, 33]

$$\begin{eqnarray}&&{\partial }_{\tau }{Z}_{\sigma }(\lambda ;\tau )=\displaystyle \sum _{\sigma }[{{\rm{\Gamma }}}_{\sigma \sigma ^{\prime} }(\tau )+\lambda {\partial }_{\tau }{\epsilon }_{\sigma }(\tau ){\delta }_{\sigma \sigma ^{\prime} }]{Z}_{\sigma ^{\prime} }(\lambda ;\tau )\end{eqnarray} \tag{ F25 }$$

and the initial condition

$$\begin{eqnarray}&&{Z}_{\sigma }(\xi ;0)={e}^{\xi {\epsilon }_{\sigma }(0)}.\end{eqnarray} \tag{ F26 }$$

Compared with the original master equation (F13) for the level-occupation probabilities, the new master equation (F25) for the characteristic function contains additional diagonal terms. The full characteristic function is

$$\begin{eqnarray}&&Z(\xi )=\displaystyle \sum _{{\sigma }_{0}}{p}_{{\sigma }_{0}}{Z}_{{\sigma }_{0}}(\xi ).\end{eqnarray} \tag{ F27 }$$

Recall that $Z(\xi )$ contains the same information as P(W). From $Z(\xi )$, one can calculate P(W) itself and, as shown in section F.4 below, a bound for ${D}_{\infty }({P}_{\mathrm{fwd}}(W)\parallel {P}_{\mathrm{rev}}(-W))$.

Let us show that the work distribution in equation (F24) satisfies Crooks' fluctuation theorem,

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{\mathrm{fwd}}(W)}{{P}_{\mathrm{rev}}(-W)}=\displaystyle \frac{{Z}_{f}}{{Z}_{0}}{e}^{\beta W},\end{eqnarray} \tag{ F28 }$$

wherein Z₀ and Z_f are the partition functions for the initial and final Hamiltonians in the forward protocol. Given a forward ramping $\epsilon (t)$, the reverse ramping ${\epsilon }^{\mathrm{rev}}(t)$ is defined by

$$\begin{eqnarray}&&{\epsilon }^{\mathrm{rev}}(t)=\epsilon (\tau -t).\end{eqnarray} \tag{ F29 }$$

In the forward protocol, consider a trajectory $\sigma (t)$ characterized by the initial condition ${\sigma }_{0}$, the number J of energy-level jumps and the jump instants t_j ($j=1,2,\ldots ,J$). One can find a unique trajectory ${\sigma }^{\mathrm{rev}}(t)$ in the reverse protocol, which is defined by the initial condition

$$\begin{eqnarray}&&{\sigma }_{0}^{\mathrm{rev}}={\sigma }_{0}+J(\mathrm{mod}\,\,2)\end{eqnarray} \tag{ F30 }$$

and the flip instants

$$\begin{eqnarray}&&{t}_{j}^{\mathrm{rev}}=\tau -{t}_{J-j+1}.\end{eqnarray} \tag{ F31 }$$

Note that

$$\begin{eqnarray}&&{\epsilon }^{\mathrm{rev}}({t}_{j}^{\mathrm{rev}})=\epsilon ({t}_{J-j+1}).\end{eqnarray} \tag{ F32 }$$

The effective action along the reverse trajectory is the same as that along the forward trajectory (see (F20)):

$$\begin{eqnarray}&&{S}_{J}^{\mathrm{rev}}({t}_{1}^{\mathrm{rev}},\ldots ,{t}_{J}^{\mathrm{rev}};{\sigma }_{0}^{\mathrm{rev}})={S}_{J}({t}_{1},\ldots ,{t}_{J};{\sigma }_{0}).\end{eqnarray} \tag{ F33 }$$

Further, the work contribution along the reverse trajectory is the negative of that along the forward trajectory (see (F22)):

$$\begin{eqnarray}&&{W}_{J}^{\mathrm{rev}}({t}_{1}^{\mathrm{rev}},\ldots ,{t}_{J}^{\mathrm{rev}};{\sigma }_{0}^{\mathrm{rev}})=-{W}_{J}({t}_{1},\ldots ,{t}_{J};{\sigma }_{0}).\end{eqnarray} \tag{ F34 }$$

These observations lead to

$$\begin{eqnarray}&&{P}_{J}^{\mathrm{rev}}({t}_{1}^{\mathrm{rev}},\ldots ,{t}_{J}^{\mathrm{rev}};{\sigma }_{0}^{\mathrm{rev}})={P}_{J}({t}_{1},\ldots ,{t}_{J};{\sigma }_{0}){e}^{-\beta {W}_{J}({t}_{1},\ldots ,{t}_{J};{\sigma }_{0})}\exp [({\sigma }_{0}+J\mathrm{mod}2)\beta {\epsilon }_{f}-{\sigma }_{0}\beta {\epsilon }_{0}]\end{eqnarray} \tag{ F35 }$$

and

$$\begin{eqnarray}&&{P}_{J}^{\mathrm{rev}}(-W;{\sigma }_{0}^{\mathrm{rev}})={P}_{J}(W;{\sigma }_{0}){e}^{-\beta W}{e}^{({\sigma }_{0}+J\mathrm{mod}2)\beta {\epsilon }_{f}-{\sigma }_{0}\beta {\epsilon }_{0}}\end{eqnarray} \tag{ F36 }$$

It is then straightforward to prove Crooks' theorem:

$$\begin{eqnarray*}{P}_{\mathrm{rev}}(-W) & = & \displaystyle \frac{1}{1+{e}^{-\beta {\epsilon }_{0}^{\mathrm{rev}}}}\displaystyle \sum _{J=0}^{\infty }\displaystyle \sum _{{\sigma }_{0}^{\mathrm{rev}}}{e}^{-\beta {\sigma }_{0}^{\mathrm{rev}}{\epsilon }_{f}^{\mathrm{rev}}}{P}_{J}^{\mathrm{rev}}(-W;{\sigma }_{0}^{\mathrm{rev}})\\ & = & \displaystyle \frac{1}{1+{e}^{-\beta {\epsilon }_{f}}}\displaystyle \sum _{J=0}^{\infty }\displaystyle \sum _{{\sigma }_{0}}{e}^{-\beta ({\sigma }_{0}+J\mathrm{mod}2){\epsilon }_{0}}\\ & & \times {e}^{-\beta W}{e}^{\beta ({\sigma }_{0}+J\mathrm{mod}2){\epsilon }_{f}-{\sigma }_{0}\beta {\epsilon }_{0}}{P}_{J}(W;{\sigma }_{0})\\ & = & \displaystyle \frac{{e}^{-\beta W}}{1+{e}^{-\beta {\epsilon }_{f}}}\displaystyle \sum _{J=0}^{\infty }\displaystyle \sum _{{\sigma }_{0}}{e}^{-{\sigma }_{0}\beta {\epsilon }_{0}}{P}_{J}(W;{\sigma }_{0})\\ & = & {e}^{-\beta W}\displaystyle \frac{1+{e}^{-\beta {\epsilon }_{0}}}{1+{e}^{-\beta {\epsilon }_{f}}}{P}_{\mathrm{fwd}}(W).\end{eqnarray*}$$

For illustration, examples of the forwards and reverse distributions appear in figure F3.

Zoom In Zoom Out Reset image size

F.4. Upper bound on ${D}_{\infty }$

Recall the Markov inequality for a non-negative random variable X:

$$\begin{eqnarray*}&&p(X\geqslant a)\leqslant \langle X\rangle /a.\end{eqnarray*}$$

This is derived by noting that there cannot be too much probability of having a value much greater than the average, or else the average would have to be greater. In our case, it reads

$$\begin{eqnarray*}&&p(w\geqslant {\widetilde{w}}^{\epsilon })=: \epsilon \leqslant \langle w\rangle /{\widetilde{w}}^{\epsilon }.\end{eqnarray*}$$

Thus,

$$\begin{eqnarray}&&{\widetilde{w}}^{\epsilon }\leqslant \langle w\rangle /\epsilon .\end{eqnarray} \tag{ F37 }$$

We rearrange the main result, equation (D1):

$$\begin{eqnarray}&&{D}_{\infty }({\widetilde{p}}_{\mathrm{fwd}}^{\epsilon }(w)| | {p}_{\mathrm{rev}}(-w))=\beta {\widetilde{w}}^{\epsilon }-\mathrm{log}(1-\epsilon )+\mathrm{log}Z/\widetilde{Z}.\end{eqnarray} \tag{ F38 }$$

Substituting in equation (F37) yields

$$\begin{eqnarray}&&{D}_{\infty }({\widetilde{p}}_{\mathrm{fwd}}^{\epsilon }(w)| | {p}_{\mathrm{rev}}(-w))\leqslant \beta \langle w\rangle /\epsilon -\mathrm{log}(1-\epsilon )+\mathrm{log}Z/\widetilde{Z}.\end{eqnarray} \tag{ F39 }$$

(here $\mathrm{log}Z/\widetilde{Z}=\mathrm{log}2$). One has only to upper-bound $\langle w\rangle$. $\langle w\rangle$ can be upper-bounded most easily with the characteristic function $\langle {e}^{\lambda w}\rangle$, which bounds $\langle w\rangle$ due to convexity. This has been illustrated in figure F4.

Zoom In Zoom Out Reset image size

We finally remark that, as shown in [6, 7], in the isothermal limit, the penalty (meaning again the lhs of equation (F38)), goes to zero. The isothermal limit means that the hopping probabilities multiplying together to give a trajectory's probability as in equation (1) take the form of thermal state occupation probabilities ${\gamma }_{j}^{{i}_{j}}:= \exp -\beta E(| {i}_{j},{\lambda }_{j}\rangle )/Z({\lambda }_{j})$. The probability of a forwards trajectory becomes $p({\rm{traj}})={\gamma }_{0}^{{i}_{0}}{\gamma }_{1}^{{i}_{1}}...{\gamma }_{f}^{{i}_{f}}$ whereas the reverse trajectory has the probability $p({\rm{traj}}-{\rm{inv}})={\gamma }_{f}^{{i}_{f}}{\gamma }_{f}^{{i}_{f-1}}...{\gamma }_{1}^{{i}_{0}}$. The probability of a given time sequence of work values from the elementary steps is then a product of individual distributions $p({w}_{1},{w}_{2},\ldots )={p}_{1}({w}_{1}){p}_{2}({w}_{2})...$. This allows one to use the McDiarmid inequality for independent random variables as in [7] to show that there is concentration around the average in the limit of breaking up an isothermal time evolution into infinitely many substeps of energy shifts. If we write ${\widetilde{w}}^{\epsilon }=\langle w\rangle +\varepsilon$, both and ε tend towards zero in this limit. Combining that with the also known fact that $\langle w\rangle =-{kT}\mathrm{log}Z/\widetilde{Z}$ in the isothermal case, we see that the rhs of equation (F38) tends to zero in this limit; thus the lhs also tends to zero. To find, for a given initial state and initial and final Hamiltonians, a protocol such that this penalty tends to zero, one can accordingly begin the protocol with lifting the OUT levels to the corresponding thermal levels, and then perform isothermal quasistatic extraction as described above. The lifting of the OUT levels then undoes their initial lowering in the ∼-protocol, undoing any work cost and returning the state to being a thermal state. This limit also illustrates why the -versions of the worst-case work and associated penalty are physically natural to introduce. Strictly speaking the worst case work w⁰ could be much larger than the average even in the isothermal case, but the probability of this happening can be arbitrarily small.