Optimal protocols and universal time-energy bound in Brownian thermodynamics

TL;DR

This paper introduces an optimal control strategy for stochastic thermodynamic systems, experimentally demonstrating a universal time-energy bound in Brownian thermodynamics that is achieved by optimized protocols.

Contribution

It develops a variational approach to optimize work and transfer time in Brownian systems, revealing a universal time-energy bound and experimentally validating it.

Findings

Universal relation between transfer duration and work expended

Optimized protocols reach the theoretical lower bound

Experimental verification with Brownian particles in optical traps

Abstract

We propose an optimization strategy to control the dynamics of a stochastic system transferred from one thermal equilibrium to another and apply it experimentally to a Brownian particle in an optical trap under compression. Based on a variational principle that treats the transfer duration and the expended work on an equal footing, our strategy leads to a family of protocols that are either optimally cheap for a given duration or optimally fast for a given energetic cost. This approach unveils a universal relation $\Delta t\,\Delta W \ge (\Delta t\,\Delta W)_{\rm opt}$ between the transfer duration and the expended work. We verify experimentally that the lower bound is reached only with the optimized protocols.

Equations

$$\gamma \frac{ds}{dt}=-2\kappa (t)s+2D\gamma ,$$

$$\Delta t[\hat{\kappa }(s)]\equiv t_f-t_i=\frac{1}{2}{\int }_{s_i}^{s_f}\frac{\gamma ds}{D\gamma -s\hat{\kappa }(s)},$$

$$W[\hat{\kappa }(s)]=-\frac{1}{2}{\int }_{s_i}^{s_f}\hat{\kappa }(s)ds+\frac{1}{2}\left({\kappa }_f{s}_f-{\kappa }_i{s}_i\right),$$

$$J[\hat{\kappa }(s)]={\int }_{s_i}^{s_f}\frac{\gamma ds}{D\gamma -s\hat{\kappa }(s)}-\lambda {\int }_{s_i}^{s_f}\hat{\kappa }(s)ds,$$

$$\frac{d}{ds}\frac{\partial L}{\partial {\hat{\kappa }}^{\prime }}-\frac{\partial L}{\partial \hat{\kappa }}=0,$$

$$s\hat{\kappa }(s)=D\gamma +\sqrt{\gamma s/\lambda }$$

$$\kappa (t)=\frac{D\gamma +\sqrt{\gamma s_i/\lambda }-t/\lambda }{{\left(\sqrt{s_i}-t/\sqrt{\gamma \lambda }\right)}^2},$$

$$W_{\mathrm{opt}}=W_{\mathrm{QS}}+\sqrt{\gamma /\lambda }\left(s_i^{1/2}-s_f^{1/2}\right).$$

$$\Delta t_{\mathrm{opt}}=\sqrt{\gamma \lambda }\left(s_i^{1/2}-s_f^{1/2}\right).$$

$$\Delta t_{\mathrm{opt}}\Delta W_{\mathrm{opt}}=\frac{\gamma }{2}{\left(s_i^{1/2}-s_f^{1/2}\right)}^2$$

$$2\epsilon \frac{d^2\hat{\kappa }}{d{s}^2}=\frac{\gamma s}{(D\gamma -s\hat{\kappa }{)}^2}-\lambda ,$$

$$s(t_j)={\sigma }_x^2(t_j)\pm \underset{\delta s(t_j)}{\underbrace{(\delta {\sigma }_{{\chi }^2}^2(t_j)+\beta \delta \beta {\sigma }_x^2(t_j))}},$$

$$⟨s{⟩}_t=\frac{1}{\Delta t}\sum _{j=1}^{n}s(t_j),$$

$${\delta }_t⟨s⟩=\sqrt{\frac{1}{\Delta t}\sum _{j=1}^n|s(t_j)-⟨s{⟩}_t{|}^2}$$

$$s_{i,f}=⟨s{⟩}_t\pm \underset{{\delta }_t{s}_{i,f}}{\underbrace{\left({\delta }_t⟨s⟩+⟨\delta s{⟩}_t+{\delta }_t⟨\delta s⟩\right)}},$$

$$L[s,\hat{\kappa }(s)]=\frac{\gamma }{D\gamma -s\hat{\kappa }(s)}-\lambda \hat{\kappa }(s).$$

$$L[s,\hat{\kappa }(s)]=\frac{\gamma }{D\gamma -s\hat{\kappa }(s)}-\lambda \hat{\kappa }(s)+\epsilon |{\hat{\kappa }}^{\prime }(s){|}^2,$$

$$2\epsilon \frac{d^2\hat{\kappa }}{d{s}^2}=\frac{\gamma s}{(D\gamma -s\hat{\kappa }{)}^2}-\lambda .$$

$$-\alpha {\hat{\kappa }}_i^{n+1}+2\epsilon {\left(\frac{d^2\hat{\kappa }}{d{s}^2}\right)}_i^{n+1}=\frac{\gamma s_i}{(D\gamma -s_i{\hat{\kappa }}_i^n{)}^2}-\lambda -\alpha {\hat{\kappa }}_i^n,$$

$${\left(\frac{d^2\hat{\kappa }}{d{s}^2}\right)}_i\approx \frac{{\hat{\kappa }}_{i-1}-2{\hat{\kappa }}_i+{\hat{\kappa }}_{i+1}}{\delta s^2}.$$

Full text

Controlling the transformation of equilibrium states (and related quantities) is a major concern in stochastic energetics. Though still in its infancy, this is an important research area with promising applications for nanotechnologies. Recently, many experimental and theoretical perspectives, both in the classical and quantum regimes Schmiedl and Seifert (2007); Chen et al. (2010); Aurell et al. (2011); Seifert (2012); Weber et al. (2014), have demonstrated the possibility to control the evolution of a small system while constraining a set of thermodynamic variables using appropriate protocols. For instance, recent work proposed protocols that can force a nano- or micro-system to evolve from one equilibrium state to another much faster than the relaxation time expected from the energy difference between the two equilibria Martínez et al. (2016); Chupeau et al. (2018); Le Cunuder et al. (2016). From a mathematical viewpoint, this is an interesting optimal control problem, which can be studied using the Pontryagin’s principle Pontryagin et al. (1962); Plata et al. (2019).

Accelerated equilibration protocols have direct thermodynamic consequences. A protocol that reduces the transfer duration is necessarily more expensive energetically, so that the requirements of being fast and cheap cannot be satisfied simultaneously Martínez et al. (2016). Earlier proposals have discussed the possibility to minimize the work expended through a transfer whose duration is initially fixed Schmiedl and Seifert (2007). But this approach does not treat duration and work on an equal footing: while the former is arbitrarily fixed by the experimentalist, only the latter is minimized. This strategy prevents one from deriving and exploiting the mutually exclusive relation between a protocol’s duration and its energetic cost, which is of paramount importance to design protocols that are optimized from both points of view. The possibility for optimal control turns out to be particularly relevant in the field of stochastic engines, where it is necessarily related to the global figure of merit of the system Blickle and Bechinger (2011); Martínez et al. (2015).

In order to derive such optimal protocols, we adopt an original approach that treats both the duration of the transfer and the expended work in a completely symmetric way. Our strategy, implemented on an optically trapped Brownian particle, is based on two novel ingredients. First, each protocol is defined by a path in the phase space $(\kappa,s)$, where $\kappa$ is the stiffness of the optical trap and $s$ is the variance of the position of the particle. Second, we construct a functional $J[\kappa,s]$ that is composed of two terms, corresponding respectively to the total work and total duration, with a Lagrange multiplier $\lambda$ regulating the trade-off between these two quantities. Minimizing the above functional with fixed boundary conditions, leads to the desired optimal protocol $\kappa_{\lambda}(s)$. For instance, $\lambda\gg 1$ yields a protocol that has a low energetic cost but long duration; conversely, $\lambda\ll 1$ leads to a fast protocol that requires a large amount of work.

Remarkably, this approach leads to a universal relation $\Delta t\,\Delta W\geq(\Delta t\,\Delta W)_{\rm opt}$ between the transfer duration $\Delta t$ and the expended work $\Delta W$ (in excess of the free energy difference), where the lower bound depends exclusively on the initial and final states and can be reached only under optimal control conditions. This result unveils a fundamental feature that underpins all optimization procedures in stochastic thermodynamics Sekimoto (2010); Ciliberto (2017).

Our Brownian particle is a polystyrene microsphere optically trapped in water at room temperature – see Appendix A for a detailed description of the setup Schnoering and Genet (2015); Schnoering et al. (2019, 2018). In this overdamped regime, the conditions are carefully set so that the trapping potential is harmonic. We record the instantaneous motion $x(t)$ of the microsphere along the optical axis of the trap. With the trap stiffness $\kappa$ proportional to the trapping laser intensity $I$, it is possible to define, by modulating $I(t)$, a given protocol $\kappa(t)$ for the transfer from an initial thermal equilibrium at time $t_{i}$ to a final equilibrium at time $t_{f}$.

By performing a series of $N$ identical protocols on the trapped microsphere, we build a statistical ensemble of trajectories that yields a probability density function (PDF) of positions $x$. The dynamics of the system is described by the variance $s(t)$ extracted from the PDF that evolves according to

[TABLE]

where $\gamma=6\pi R\eta$ is the Stokes drag coefficient, which depends on the radius of the particle $R=500$ nm and the dynamic viscosity of the fluid $\eta\sim 10^{-3}\ {\rm Pa}\,{\rm s}$, and $D=k_{B}T/\gamma\sim 0.4\ \mu{\rm m}^{2}/{\rm s}$ the Brownian diffusion coefficient fixed by the temperature $T$ of the fluid and the Boltzmann constant $k_{B}$.

Equation (1) fully determines the statistical properties of the system where the initial and final equilibria correspond to the stationary solutions $s_{i}\kappa_{i}=s_{f}\kappa_{f}=k_{B}T$ with Gaussian PDF. It is also clear that the PDF will remain Gaussian for all intermediate times between $t_{i}$ and $t_{f}$. The cumulative energetics involved during the protocol is directly related to the time evolution of $s(t)$, giving the ensemble-averaged expended work $W(t)=\frac{1}{2}\int_{t_{i}}^{t}{\rm d}ts(t)\dot{\kappa}(t)$ and dissipated heat $Q(t)=-\frac{1}{2}\int_{t_{i}}^{t}{\rm d}t\dot{s}(t)\kappa(t)$ (following the convention of a positive flow when heat is transferred from the trapped microsphere to the bath) Sekimoto (1998); Ciliberto (2017).

Our purpose is to control the dynamics of the microsphere so that the transfer between the two equilibria is optimal with respect to both duration and energetics. If one switches instantaneously the trap stiffness from $\kappa_{i}$ to $\kappa_{f}>\kappa_{i}$ (i.e., closing the trap in a step-like way), the typical relaxation time to the new equilibrium is given by $\tau_{\rm relax}=2\gamma/\kappa_{f}$. This time will be taken as the reference value with respect to which the reduction in transfer duration is measured.

Our optimization strategy starts with the idea of using the variance $s$ as the independent variable of the problem, instead of the time $t$. This is possible whenever the function $s(t)$ is monotonic and it enables us to express the control parameter as $\kappa(t)=\hat{\kappa}\left(s(t)\right)$. The advantage of this approach is that we can easily write down, as functionals of $\hat{\kappa}(s)$, both the transfer duration:

[TABLE]

and the expended work:

[TABLE]

where the second term on the right-hand side vanishes because the initial and final configurations are thermally equilibrated.

Next, we define the functional to be minimized as twice the sum of $W$ and $\Delta t$:

[TABLE]

where $\lambda$ is a Lagrange multiplier that regulates the trade-off between transfer duration and work. Within this framework, the optimization strategy can be interpreted as the search for the trajectory in the $(s,\kappa)$-space that minimizes $J[\hat{\kappa}(s)]$ while keeping the extrema fixed at equilibrium, i.e. $s_{i}\kappa_{i}=s_{f}\kappa_{f}=k_{\rm B}T=D\gamma$. Once written as $J=\int_{s_{i}}^{s_{f}}L[s,\hat{\kappa}(s)]ds$, this functional can be minimized using the standard Euler-Lagrange equation

[TABLE]

where $\hat{\kappa}^{\prime}\equiv d\hat{\kappa}/ds$, yielding the following solution

[TABLE]

for a protocol that eventually closes the trap with $\kappa_{f}>\kappa_{i}$ 111In order to satisfy $\Delta t>0$, a negative sign should be used on the r.h.s. of Eq. (6) when $s_{f}>s_{i}$ (opening trap, $\kappa_{f}<\kappa_{i}$), whereas a positive sign should be used when $s_{f}<s_{i}$ (closing trap, $\kappa_{f}>\kappa_{i}$)..

Equation (6) encapsulates the main result obtained so far. For instance, the quasi-static solution ($\dot{s}\approx 0$) is obtained by taking $\lambda\to\infty$, which yields an infinite duration but the smallest possible expended work $W_{\rm QS}={1\over 2}D\gamma\ln(\kappa_{f}/\kappa_{i})$, equal to the free energy difference between the two equilibria, as expected for a quasi-static process. For finite $\lambda$ and making use of Eq. (6), Eq. (1) can be rewritten as $\dot{s}=-2\sqrt{s/\gamma\lambda}$, which possesses the general solution $s(t)=(\sqrt{s_{i}}-t/\sqrt{\gamma\lambda})^{2}.$ Inserting this expression into Eq. (6) yields the optimal evolution of the trap stiffness

[TABLE]

which defines our protocol for the optimized transfer.

It is important to stress that, in our case, the Euler-Lagrange equation (5) is purely algebraic, so that one cannot enforce the initial and final conditions. Thus, except for the quasi-static limit, Eq. (7) does not satisfy the conditions for which the system is at thermal equilibrium in the initial and final states. In order to ensure that $s_{i,f}\kappa_{i,f}=D\gamma$, we need to add to Eq. (7) two discontinuities (as already noticed in Schmiedl and Seifert (2007, 2008)). The optimal protocol thus consists of three successive sequences:

At time $t=t_{i}$, the trap stiffness is suddenly changed from $\kappa_{i}=D\gamma/s_{i}$ (initial equilibrium) to $\kappa(t_{i}^{+})=\kappa_{i}^{+}$, such that: $\kappa_{i}^{+}-\kappa_{i}=\sqrt{\gamma/(\lambda s_{i})}$, while keeping the variance equal to $s_{i}$; 2. 2.

Between $t_{i}^{+}$ and $t_{f}^{-}$, the stiffness varies according to Eq. (7), reaching $\kappa(t_{f}^{-})\equiv\kappa_{f}^{-}=D\gamma/s_{f}+\sqrt{\gamma/(\lambda s_{f})}$; 3. 3.

At time $t=t_{f}$, the stiffness is suddenly changed from $\kappa_{f}^{-}$ to $\kappa_{f}=D\gamma/s_{f}$ (final equilibrium), while keeping the variance equal to $s_{f}$.

Experimentally, we have performed $N=2\times 10^{4}$ successive and identical optimal protocols on the trapped microsphere, each built on the above three sequences, forcing the system to relax to thermal equilibrium within a time $\Delta t$ chosen to be shorter than $\tau_{\rm relax}$. The time-evolution of $s$ between two thermal equilibrium configurations is displayed in Fig. 1 for a shortened duration $\Delta t\sim\tau_{\rm relax}/10$, where $\tau_{\rm relax}=2\gamma/\kappa_{f}=3.22\pm 0.09\ {\rm ms}$. The reduction in $s(t)$ from its initial value corresponds to the fact that the trap is stiffer at the end of the protocol with $\kappa_{f}/\kappa_{i}\sim 1.85$. The time evolution of $s$ calculated using these values for $s_{i}$ and $\lambda$ is in excellent agreement with the experiment.

This reduction in the transfer duration has an energetic cost that can be evaluated for each sequence using Eq. (3). Such cost is measured experimentally by evaluating the cumulative work $W(t)$ from the recorded evolution of $\hat{\kappa}(s)$. As seen in Fig. 2, the time-evolution of $W(t)$ can also be split into three sequences. First, the quantity of work $W^{(1)}=s_{i}(\kappa_{i}^{+}-\kappa_{i})/2$ is injected instantaneously into the system at the time $t_{i}$ as the trap is suddenly stiffened from $\kappa_{i}$ to $\kappa_{i}^{+}$. During the second sequence, the injection of work continues as the trapping volume is progressively reduced, reaching $W^{(2)}=W_{\rm QS}+\sqrt{\gamma/\lambda}\left[(s_{i}^{3/2}-s_{f}^{3/2})/3-(s_{i}^{1/2}-s_{f}^{1/2})/2\right]$ at $t=t_{f}^{-}$, when $\kappa=\kappa_{f}^{-}$. Finally, the trap is suddenly expanded at $t=t_{f}$, and the system instantly reaches its final equilibrium state, delivering to the thermal bath a quantity of work equal to $W^{(3)}=s_{i}(\kappa_{f}-\kappa_{f}^{-})/2$. For $t>t_{f}$, the thermal steady state is characterized by $W=Q$, with $W(t>t_{f})=(0.981\pm 0.059)\ k_{\rm B}T$ and $Q(t>t_{f})=(0.983\pm 0.060)\ k_{\rm B}T$.

In contrast, heat is continuously dissipated from the microsphere to the thermal bath, as seen in Fig. 2 from the monotonic increase of the ensemble-averaged cumulative heat $Q(t)$ throughout the protocol. The evolution of the dissipated heat between the two equilibrium states is almost exactly linear in time, which corresponds to a constant production of entropy. We stress again that the experimentally measured values of the heat and work injected in and extracted from the system are in excellent agreement with the theoretical predictions.

The total work expended throughout the optimal protocol is evaluated by adding the contributions from each of the three steps described above. One obtains:

[TABLE]

Similarly, the total duration of the optimal protocol is obtained by inserting Eq. (6) into Eq. (2), yielding:

[TABLE]

The above expressions clearly show that our optimization procedure is perfectly symmetric as far as duration and work are concerned and that the trade-off between these two quantities is governed by the Lagrange multiplier $\lambda$. Indeed, one can choose $\lambda$ using Eq. (9) to fix the total duration and then the minimum expended work will be given by Eq. (8); or, alternatively, one can determine $\lambda$ through Eq. (8) to fix the total work and then the minimum duration of the process will be given by Eq. (9).

This leads us to define the “excess work” of the optimal protocol as $\Delta W_{\rm opt}\equiv W_{\rm opt}-W_{\rm QS}$ and to note that the product

[TABLE]

is independent of $\lambda$ and only depends on the initial and final states. This equality fixes the mutually exclusive relation between transfer duration and expended work under optimal control. It corresponds to the frontier value of a universal exclusion region $\Delta t\,\Delta W\geq\gamma/2(\sqrt{s_{i}}-\sqrt{s_{f}})^{2}$ that bounds from below all protocols that are not optimal.

The frontier can be drawn experimentally by changing the transfer duration $\Delta t$ within the conditions of optimal control, i.e. changing the Lagrange multiplier. To do so, we have measured $\Delta W_{\rm opt}$ for a series of eight optimal protocols with different $\Delta t$. By normalizing each measured value of $\Delta W_{\rm opt}$ to the associated value of $\gamma/2(\sqrt{s_{i}}-\sqrt{s_{f}})^{2}$, one can test the universal nature of the bound. This is clearly confirmed in Fig. 3, with all the optimal solutions implemented experimentally falling precisely on the $1/\Delta t$ curve. To further prove that the frontier corresponds to a lower-bound, we have verified experimentally that the ($\Delta t,\Delta W)$ coordinates of typical non-optimal protocols – continuous (see below), step-like, and “engineered swift equilibration” protocols (see Appendix B) – all fall above the expected bound, as displayed in Fig. 3.

A salient feature of our optimal control procedure is represented by the sudden jumps in stiffness that have to augment the solution of Eq. (7) in order to comply with thermally equilibrated initial and final configurations. From an experimental point of view, such discontinuities do not constitute a weakness of the procedure, as they correspond to finite and measurable quantities of work exchanged between the bath and the system Aurell et al. (2012); Plata et al. (2019). But it is interesting to stress that one asset of our variational strategy is its capacity to construct smooth protocols that are as close as desired to the optimal ones. For this, we need to control the derivatives of the function $\hat{\kappa}(s)$, which can be done by adding the gradient term $\int_{s_{i}}^{s_{f}}|\hat{\kappa}^{\prime}(s)|^{2}ds$ to the functional $J[\hat{\kappa}(s)]$ in Eq. (4), with a second Lagrange multiplier $\varepsilon$. Hence, we arrive at the modified Euler-Lagrange equation:

[TABLE]

which can be solved numerically as a boundary value problem, with initial and final conditions at thermal equilibrium $\hat{\kappa}(s_{i,f})=D\gamma/s_{i,f}$. Once the solution $\hat{\kappa}(s)$ is known, the time-evolution of the variance $s(t)$ is found by integrating Eq. (1) (more details are given in Appendix C).

Using the same values of $\lambda$ that defined the optimal protocols with, respectively, $\Delta t\sim\tau_{\rm relax}/22$ and $\Delta t\sim\tau_{\rm relax}/10$ (Fig. 3, inset), we implemented two smooth protocols for two different values of the Lagrange multiplier $\varepsilon=5\times 10^{-6}$ and $\varepsilon=10^{-6}$ (here and in the following, $\varepsilon$ is expressed in units of $s_{i}^{2}/(D\kappa_{i}^{2})$). As shown in Fig. 4, the smooth protocols follow closely the optimal ones, except near the beginning and the end of the process, where they approach the equilibrium states in a continuous way. For the same value of $\lambda$, smooth protocols give slightly longer transfer durations ($2.48\times 10^{-4}\ {\rm s}$ for $\varepsilon=5\times 10^{-6}$ and $2.14\times 10^{-4}\ {\rm s}$ for $\varepsilon=10^{-6}$) than the optimal protocol ($\Delta t=1.72\times 10^{-4}\ {\rm s}$) but, as expected, the expended work is slightly smaller ($1.36\pm 0.06~{}k_{\rm B}T$ for $\varepsilon=5\times 10^{-6}$ and $1.65\pm 0.06~{}k_{\rm B}T$ for $\varepsilon=10^{-6}$) in the smooth case than in the optimized limit ($1.69\pm 0.06~{}k_{\rm B}T$). The non-optimal character of the smooth protocols is clearly seen in the insets of Fig. 3, where all ($\Delta t,\Delta W)$ coordinates lie above the universal bound, and only converge towards it in the $\varepsilon\rightarrow 0$ limit.

In conclusion, we have devised a family of optimal protocols that transfer an optically trapped microsphere between two equilibrium positions, minimizing both the transfer duration and the associated energetic cost. Within such protocols, the trade-off between duration and work can be modulated at will by tuning a single Lagrange multiplier given by our variational approach. A key result of our work is to show that the product $\Delta t\,\Delta W$ is bounded from below, in a way reminiscent of energy-time uncertainty relations. Similar bounds were noticed in earlier works Sekimoto (2010); Ciliberto (2017), but only for some special cases. Here, our bound is universal (it depends exclusively on the initial and final states) and is only reached for the optimal protocol, as we demonstrated both theoretically and experimentally. Further extending the present results to quantum systems may open new interesting perspectives in the burgeoning field of quantum stochastic thermodynamics Elouard et al. (2017); Vinjanampathy and Anders (2016); Roßnagel et al. (2016); Cavina et al. (2018).

I Acknowledgments

This work was supported in part by Agence Nationale de la Recherche (ANR), France, ANR Equipex Union (Grant No. ANR-10-EQPX-52-01), the Labex NIE projects (Grant No. ANR-11-LABX-0058-NIE), and USIAS within the Investissements d’Avenir program (Grant No. ANR-10-IDEX-0002-02). Y. R.-C. is a member of the International Doctoral Program of the Initiative d’Excellence of the University of Strasbourg, whose support is acknowledged. L. M. is supported by the National Natural Science Foundation of China, Research Fund for International Young Scientists under the project No. 1161101053 and the Young Scientist Program under the project No. 11601335.

II Appendix A: Setup, calibration, uncertainties

II.1 Optical trap setup

All experiments are performed on single optically trapped polystyrene spheres (radius $R=500$ nm) taken from a monodisperse ($\delta R/R=0.028$) solution (ThermoFisher, FluoSpheres) and enclosed inside a fluidic cell filled with dionized water. The microfluidic cell is made with a microscope slide and a 170 $\mu$m thick glass coverslip, sealed with a 120 $\mu$m thick spacer.

The optical trap, described in details in Fig. 5, is an evolution of the setup described in our previous work Schnoering and Genet (2015); Schnoering et al. (2019, 2018). It uses a CW near-infrared ($\lambda_{T}=$785 nm) laser whose intensity – hence the trap stiffness – can be modulated externally using a waveform generator. Any trapping protocol can then be implemented by computer-programming the waveform generator so that the time-evolution of the trap stiffness follows the desired profile.

Under such trapping laser modulation, the instantaneous axial motion $x(t)$ of the bead is monitored using an auxiliary laser propagating in the opposite direction of the trapping laser (see Fig. 5). We checked that this low-power probe beam, injected in the fluidic cell from its back-side, does not exerts any spurious optical force of the trapped bead. The signal collected by the photodiode and the output voltage of the waveform generator are simultaneously registered by a multichannel acquisition card (National Instruments, NI-6251) with a sampling rate $f_{s}=$ 218 Hz. In order to span the signal in the full dynamic range of the acquisition card, the generator output voltage was re-scaled using a scaling amplifier (Stanford Research Systems, SIM983) and the voltage time series of the photodiode was amplified and filtered using low-noise pre-amplifiers (Stanford Research Systems, SR560).

II.2 Stiffness modulation calibration

The trapping laser is modulated according to a given protocol $\kappa(t)$, defined and calculated with chosen transition parameters ($\kappa_{i},\kappa_{f},\Delta t$). In order to convert this protocol $\kappa(t)$ into a modulating voltage $V_{\rm mod}(t)$ for the waveform generator, a calibration procedure is performed. This procedure consists in measuring the trap stiffnesses associated with a series of consecutive values of DC voltages, i.e. consecutive trapping laser intensities. Each stiffness is extracted from a Lorentzian fit of the corresponding motional power spectral density (PSD) of the trapped bead. Associated error bars are obtained from the uncertainties of the Lorentzian fits (MATLAB Levenberg-Marquardt algorithm). The calibration curve shown in Fig. 6 corresponds to a linear fit of the evolution of such measured stiffnesses (including their error bars) as a function of the DC voltages.

II.3 Monitoring Brownian dynamics

The time evolution of the Brownian system is monitored by recording the stochastic trajectory of the trapped bead over $2\times 10^{4}$ cycles of the protocol $\kappa(t)$. Each cycle lasts 50 ms, where the first 30 ms correspond to the initial thermal equilibrium with $\kappa_{i}$ and the remaining ($20-\Delta t$) ms correspond to the final thermal equilibrium at $\kappa_{f}$. Each stationary region of the full trajectory, i.e. corresponding to a constant $\kappa$ ($\kappa_{i}$ or $\kappa_{f}$ ), is sectioned and concatenated with all the other sliced trajectories under the same stiffness. The PSD of this concatenated trajectory is computed and a Lorentzian fit yields the ensemble average $\kappa$. Figs. 7 (a) and (b) respectively show the PSD of the concatenated trajectories for the equilibria $\kappa_{i}$ and $\kappa_{f}$ for the case $\Delta t\sim\tau_{\rm relax}/10$ described in the main text.

Implementing the same procedure, the full temporal trace of the particle positions undergoing $2\times 10^{4}$ cycles is chopped into trajectories that correspond to a single cycle of the protocol $\kappa(t)$. The ensemble of traces then consists of all the sub-trajectories superimposed within the same time interval, in such a way that they all start $t=-30$ ms with $\kappa_{i}$, as displayed in Fig. 8 below.

The instantaneous ensemble variance $s(t_{j})$ at a time $t=t_{j}$ ($j=1,\cdots,T\times f_{s}$), with $T=50$ ms and $f_{s}=2^{18}$ Hz) is obtained by a vertical cross-cut of the ensemble of trajectories plotted in Fig. 8. The resulting distribution of positions $\rho(x,t_{j})$ is a Gaussian of zero mean $\mu_{x}(t_{j})$ and variance $s(t_{j})$. Fig. 9 displays the position distribution functions (PDF) before (equilibrium at $\kappa_{i}$) and after (equilibrium at $\kappa_{f}$) the change in trapping stiffness imposed by the protocol $\kappa(t)$. The corresponding trapping potentials calculated as $U(x,t_{j})=-k_{\rm B}T\log(\rho(x,t_{j}))+{\rm cst}$ are also shown and compared to the expected harmonic profiles $U=\frac{1}{2}\kappa x^{2}$ evaluated from the stiffnesses $\kappa_{i},\kappa_{f}$ that were extracted from the measured PSD shown in Fig 7.

Proceeding in the same manner but for all times $t_{j}$, we can obtain the temporal evolution of the ensemble variance $s(t)$ over the full protocol $\kappa(t)$. To confirm that all PDF remain Gaussian for all times, we calculate their kurtosis and verify -see Fig. 9, bottom panel- that all-time kurtosis remain very close to $3$ throughout the entire protocol.

II.4 Statistical uncertainties

The uncertainties for the instantaneous ensemble variances are obtained following a $\chi^{2}$ law with $N-1$ degrees of freedom where $N=N_{cycles}$ is the number of independent trajectories $x_{i}(t)$ undergoing one cycle of the protocol $\kappa(t)$.

II.5 PSD calibration uncertainties

Under a trapping laser intensity, the registered p-i-n voltage values $V(t)$ that correspond to the position fluctuations of the trapped bead are converted into displacement units using the best-fit parameter of the Lorentzian fit of the PSD of the trajectory (at constant $\kappa$). The fit parameter $D_{\rm fit}$ is compared to the diffusion coefficient $D=k_{B}T/\gamma$ expected from the Fluctuation-Dissipation Theorem, assuming known temperature and viscosity. This gives a conversion factor $\beta=\sqrt{D/D_{\rm fit}}$ from p-i-n voltages to meters. The uncertainty on the position sensitivity is obtained from standard error propagation including the uncertainty on the viscosity resulting from the $\delta R/R=2.8\%$ size dispersion deviation of the trapped beads.

Instantaneous positions are thus given from the conversion factor as $x(t)=(\beta\pm\delta\beta)V(t)$, and therefore the variance, up to first-order in uncertainty, $x^{2}(t)=(\beta^{2}\pm 2\beta\delta\beta)V^{2}(t)$, (since $\mu_{x}(t)=0$). The total error of the variance writes as:

[TABLE]

where $\sigma^{2}_{x}(t_{j})=\sum^{N}_{i=1}|x_{i}(t_{j})-\mu(t_{j})|^{2}/(N-1)$ is the estimator of the instantaneous ensemble variance over $N$ cycles, $\delta\sigma^{2}_{\chi^{2}}$ corresponds to the statistical uncertainty in the motional variance determination (see above) and $\delta\beta\sigma^{2}_{x}$ the PSD calibration uncertainty just discussed.

The temporal average variances related to the initial an final stiffness $s_{i}$ and $s_{f}$ are obtained from temporal average. Assuming $\Delta t$ as the interval over which $\kappa(t)$ remains constant (either at $\kappa_{i}$ or $\kappa_{f}$), the temporal average of the corresponding variance is:

[TABLE]

taking $\Delta t$ as the interval over which $\kappa(t)$ remains constant (either at $\kappa_{i}$ or $\kappa_{f}$) and $n=\Delta t\cdot f_{s}$ with $f_{s}=2^{18}$ Hz, the sampling frequency. The standard deviation of the temporal average is simply evaluated as:

[TABLE]

The stationary variances $s_{i}$ and $s_{f}$ and their uncertainties are thus simply given by:

[TABLE]

where $\langle\delta s\rangle_{t}=1/\Delta t\sum^{n}_{j=1}\delta s(t_{j})$.

II.6 Energetics uncertainties

The confidence interval of the mean cumulative work are computed taking into account the uncertainties related to both variances and stiffnesses. They are displayed on all energetic figures at a $95\%$ confidence level.

III Appendix B: Comparing optimal, step-like and ESE protocols

We compare here three protocols that transfer the bead between two equilibria, going from an initial stiffness $\kappa_{i}$ to a final one $\kappa_{f}$ with, for all protocols, fixed and identical $\kappa_{f},\kappa_{i}$ values given in the main text.

The first protocol consists of a sudden step-like change of the optical trap stiffness – see Fig. 10, green trace. The second protocol is the “engineered swift equilibriation” (ESE) protocol recently proposed and implemented by Martinez, et al. Martínez et al. (2016). We calculate $\kappa_{\rm ESE}(t)$ following Martínez et al. (2016) for a transfer duration of $\Delta t=3.47\times 10^{-4}$ s. Over the same transfer duration, we also implement our optimal protocol $\kappa_{\rm opt}(t)$. All protocols are displayed in Fig. 10.

Fig. 11 gathers the time evolutions of the motional variances associated with each protocol. As expected, the step-like protocol displays the longest equilibration time when compared to the ESE and optimal protocols. From an energetic point of view, the comparison between the two latter protocols, shown in Fig. 12, clearly reveals the non-optimal character of the ESE protocol with a cumulated work expense larger than for the the optimal protocol. This can also be seen in the inset of Fig. 12 where the excess work expended during the ESE protocol lies clearly above the optimal lower bound discussed in the main text.

IV Appendix C: Smooth protocols

The optimal protocol obtained in this work [Eq. (6) in the main text] was derived using the Lagrangian density

[TABLE]

A peculiar feature of $L[s,\hat{\kappa}(s)]$ is that the corresponding Euler-Lagrange equation is purely algebraic (as opposed to a differential equation). Hence, it is not possible to impose the desired boundary conditions on the control parameter $\hat{\kappa}(s)$ (i.e. $s_{i}\kappa_{i}=s_{f}\kappa_{f}=D\gamma$) and two jumps have to be added “by hand” at the beginning and the end of the protocol, as explained in the main text.

Although these jumps can be realized without much trouble in the experiments, it is interesting to develop a theoretical procedure capable of furnishing a suboptimal protocol $\hat{\kappa}(s)$ that is continuous in the variable $s$ and converges towards the optimal protocol as some parameter tends to zero. To do this, we need to limit the gradient of $\hat{\kappa}(s)$ by adding a further term to the Lagrangian density (16), which becomes:

[TABLE]

where $\varepsilon$ is an additional Lagrange multiplier. The above Lagrangian density yields the Euler-Lagrange equation (11) in the main text, which we reproduce here:

[TABLE]

As a second-order differential equation, Eq. (18) needs two independent boundary conditions, thus enabling us to set $s_{i}\kappa_{i}=s_{f}\kappa_{f}=D\gamma$, as requested for our protocols. When $\varepsilon\to 0$, we obtain the correct limit case of Eq. (6) in the main text, i.e., the optimal protocol containing two points of infinite derivative (jumps) for the function $\hat{\kappa}(s)$ at $s_{i}$ and $s_{f}$. Through the Lagrange multiplier $\varepsilon$, one can limit the value of such derivative, so that the protocol becomes smoother and smoother as $\varepsilon$ increases.

Equation (18) can be solved numerically by successive iterations. We used the following scheme:

[TABLE]

where the superscript $n$ denotes the $n$-th iteration, while the subscript $i$ refers to the discrete grid $s_{i}=i\,\delta s$, with spacing equal to $ds$. The second derivative is then approximated with the standard finite-difference formula:

[TABLE]

The parameter $\alpha>0$ is needed to ensure the convergence of the iterative procedure, but does not affect the final result (indeed it disappears from Eq. (19) when $\hat{\kappa}_{i}^{n+1}=\hat{\kappa}_{i}^{n}$).

As an example, we have solved Eq. (18) with physical parameters $D=\gamma=1$ and $\lambda=0.81$, corresponding to a total duration for the optimal protocol $\Delta t_{\rm opt}\sim\tau_{\rm relax}/6$ according to Eq. (9) in the main text. The boundary values are $s_{i}=1$ and $s_{f}=0.5$, $\kappa_{i}=1$ and $\kappa_{f}=2$. The smoothness parameter is $\varepsilon=10^{-5}$. The numerical convergence parameter is set to $\alpha=0.3$. The result of the numerical integration is given in Figs. 13 and 14, for both the optimal (black lines) and smooth (red lines) protocols. As expected, the smoothed protocol follows closely the optimal one, except near the extremities where it reaches its boundary values smoothly and without jumps. The total time of the smoothed protocol is $0.182\times\tau_{\rm relax}$, longer than that of the optimal one. But the total work is smaller $W_{\rm smooth}=1.32<W_{\rm opt}=1.38$. The time-energy product is $(\Delta t\,\Delta W)_{\rm smooth}=0.356>(\Delta t\,\Delta W)_{\rm opt}=0.343$, in agreement with the theoretical considerations detailed in the main text.