Optimal protocols and universal time-energy bound in Brownian thermodynamics

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.

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 [1][2][3][4][5], 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 [6][7][8]. From a mathematical viewpoint, this is an interesting optimal control problem, which can be studied using the Pontryagin's principle [9,10].
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 [6]. Earlier proposals have discussed the possibility to minimize the work expended through a transfer whose duration is initially fixed [1]. 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 [11,12].
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 (κ, s), where κ is the stiffness of the optical trap and s is the variance of the position of the particle. Second, we construct a functional J[κ, s] that is composed of two terms, corresponding respectively to the total work and total duration, with a Lagrange multiplier λ regulating the trade-off between these two quantities. Minimizing the above functional with fixed boundary conditions, leads to the desired optimal protocol κ λ (s). For instance, λ 1 yields a protocol that has a low energetic cost but long duration; conversely, λ 1 leads to a fast protocol that requires a large amount of work.
Remarkably, this approach leads to a universal relation ∆t ∆W ≥ (∆t ∆W ) opt between the transfer duration ∆t and the expended work ∆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 [13,14].
Our Brownian particle is a polystyrene microsphere optically trapped in water at room temperature -see Appendix A for a detailed description of the setup [15][16][17]. 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 κ proportional to the trapping laser intensity I, it is possible to define, by modulating I(t), a given protocol κ(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 where γ = 6πRη is the Stokes drag coefficient, which depends on the radius of the particle R = 500 nm and the dynamic viscosity of the fluid η ∼ 10 −3 Pa s, and D = k B T /γ ∼ 0.4 µm 2 /s the Brownian diffusion coefficient fixed by the temperature T of the fluid and the Boltzmann constant k B . arXiv:1906.04171v1 [cond-mat.stat-mech] 9 Jun 2019 Equation (1) fully determines the statistical properties of the system where the initial and final equilibria correspond to the stationary solutions s i κ i = s f κ 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) = 1 2 t t i dts(t)κ(t) and dissipated heat Q(t) = − 1 2 t t i dtṡ(t)κ(t) (following the convention of a positive flow when heat is transferred from the trapped microsphere to the bath) [14,18].
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 κ i to κ f > κ i (i.e., closing the trap in a step-like way), the typical relaxation time to the new equilibrium is given by τ relax = 2γ/κ 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 κ(t) =κ (s(t)). The advantage of this approach is that we can easily write down, as functionals ofκ(s), both the transfer duration: and the expended work: 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 ∆t: where λ 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, κ)-space that minimizes J[κ(s)] while keeping the extrema fixed at equilibrium, i.e.
whereκ ≡ dκ/ds, yielding the following solution sκ(s) = Dγ + γs/λ for a protocol that eventually closes the trap with κ f > κ i [19]. Equation (6) encapsulates the main result obtained so far. For instance, the quasi-static solution (ṡ ≈ 0) is obtained by taking λ → ∞, which yields an infinite duration but the smallest possible expended work W QS = 1 2 Dγ ln(κ f /κ i ), equal to the free energy difference between the two equilibria, as expected for a quasi-static process. For finite λ and making use of Eq. (6), Eq. (1) can be rewritten asṡ = −2 s/γλ , which possesses the general solution s(t) = ( √ s i −t/ γλ ) 2 . Inserting this expression into Eq. (6) yields the optimal evolution of the trap stiffness 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 κ i, f = Dγ, we need to add to Eq. (7) two discontinuities (as already noticed in [1,20]). The optimal protocol thus consists of three successive sequences: 1. At time t = t i , the trap stiffness is suddenly changed while keeping the variance equal to s i ; 2. Between t + i and t − f , the stiffness varies according to Eq. (7), reaching κ(t − f ) ≡ κ − f = Dγ/s f + γ/(λ s f ); 3. At time t = t f , the stiffness is suddenly changed from κ − f to κ f = Dγ/s f (final equilibrium), while keeping the variance equal to s f . Experimentally, we have performed N = 2 × 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 ∆t chosen to be shorter than τ relax . The time-evolution of s between two thermal equilibrium configurations is displayed in Fig. 1 for a shortened duration ∆t ∼ τ relax /10, where τ relax = 2γ/κ f = 3.22 ± 0.09 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 κ f /κ i ∼ 1.85. The time evolution of s calculated using these values for s i and λ 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κ(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 (κ + i − κ i )/2 is injected instantaneously into the system at the time t i as the trap is suddenly stiffened from κ i to κ + i . During the second sequence, the injection of work continues as the trapping volume is progressively reduced, reaching W (2) 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 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: Similarly, the total duration of the optimal protocol is obtained by inserting Eq. (6) into Eq. (2), yielding: 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 λ . Indeed, one can choose λ 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 λ 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 ∆W opt ≡ W opt −W QS and to note that the product is independent of λ 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 ∆t ∆W ≥ γ/2( √ s i − √ s f ) 2 that bounds from below all protocols that are not optimal.
The frontier can be drawn experimentally by changing the transfer duration ∆t within the conditions of optimal control, i.e. changing the Lagrange multiplier. To do so, we have measured ∆W opt for a series of eight optimal protocols with different ∆t. By normalizing each measured value of ∆W opt to the associated value of γ/2( Extracted excess works ∆W opt for a series (square) of optimal protocols defined by transfer durations ∆t = τ relax /n, with successively n ∼ 34, 30, 22, 16, 10, 6, 3 and 2, revealing the mutually exclusive relation between ∆W opt and ∆t. For each n− protocol, we normalize ∆W opt by the corresponding value γ/2( √ s i − √ s f ) 2 , considering that the precise values for (s i , s f ) slightly vary from protocol to protocol. The universality of the bound is clearly verified experimentally by observing that all optimized coordinates {∆t, } precisely fall (within error bars) on the 1/∆t curve. The excess work measured for an "engineered swift equilibration" protocol [6] (black diamond) defined for ∆t = 3.47 × 10 −4 s ∼ τ relax /10, and the excess work measured for a steplike protocol (black star) at τ relax ∼ 3.22 ms clearly fall above the optimal bound -see Appendix B for details. Insets: (∆t, ∆W ) coordinates measured for smooth (thus suboptimal) protocols for n ∼ 22 and n ∼ 10, and smoothness parameters ε = 5 × 10 −6 , ε = 10 −6 , and ε = 0, expressed in units of s 2 i /(Dκ 2 i ). Such smooth protocols are defined using the same Lagrange multiplier λ as their associated optimal protocols. For each case, the product ∆t∆W converges to the optimal lower bound (solid blue line) as ε → 0. For n ∼ 10, the excess work for the ESE protocol is plotted again (black diamond). universal nature of the bound. This is clearly confirmed in Fig. 3, with all the optimal solutions implemented experimentally falling precisely on the 1/∆t curve. To further prove that the frontier corresponds to a lower-bound, we have verified experimentally that the (∆t, ∆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 [10,21]. 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κ(s), which can be done by adding the gradient term s f s i |κ (s)| 2 ds to the functional J[κ(s)] in Eq. (4), with a second Lagrange multiplier ε. Hence, we arrive at the modified Euler-Lagrange equation: which can be solved numerically as a boundary value problem, with initial and final conditions at thermal equilibrium κ(s i, f ) = Dγ/s i, f . Once the solutionκ(s) is known, the timeevolution of the variance s(t) is found by integrating Eq. (1) (more details are given in Appendix C).

FIG. 4. Comparison of the cumulative energetics [expended work
W (t) (blue curves) and dissipated heat Q(t) (red curves)] between an optimal protocol and two smooth protocols with ε = 5 × 10 −6 and ε = 10 −6 , expressed in units of s 2 i /(Dκ 2 i ), and identical value of λ = (2.97 ± 0.12) × 10 16 s/J. As seen on the insets of Fig. 3, although the smooth protocols involve slightly less work than the optimal one, they correspond to longer transfer durations. Inset: superimposed evolutions s(t) vs. κ(t) for the three protocols, showing the continuous nature of the smooth protocol and illustrating the progressive convergence to the optimal protocol in the ε → 0 limit. For each protocol, the curves are normalized to the corresponding κ i for κ(t) and s f for s(t).
Using the same values of λ that defined the optimal protocols with, respectively, ∆t ∼ τ relax /22 and ∆t ∼ τ relax /10 (Fig. 3, inset), we implemented two smooth protocols for two different values of the Lagrange multiplier ε = 5 × 10 −6 and ε = 10 −6 (here and in the following, ε is expressed in units of s 2 i /(Dκ 2 i )). 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 λ , smooth protocols give slightly longer transfer durations (2.48 × 10 −4 s for ε = 5 × 10 −6 and 2.14 × 10 −4 s for ε = 10 −6 ) than the optimal protocol (∆t = 1.72 × 10 −4 s) but, as expected, the expended work is slightly smaller (1.36±0.06 k B T for ε = 5×10 −6 and 1.65 ± 0.06 k B T for ε = 10 −6 ) in the smooth case than in the optimized limit (1.69 ± 0.06 k B T ). The non-optimal character of the smooth protocols is clearly seen in the insets of Fig. 3, where all (∆t, ∆W ) coordinates lie above the universal bound, and only converge towards it in the ε → 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 ∆t ∆W is bounded from below, in a way reminiscent of energy-time uncertainty relations. Similar bounds were noticed in earlier works [13,14], 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 [22][23][24][25]. All experiments are performed on single optically trapped polystyrene spheres (radius R = 500 nm) taken from a monodisperse (δ R/R = 0.028) solution (ThermoFisher, Flu-oSpheres) and enclosed inside a fluidic cell filled with dionized water. The microfluidic cell is made with a microscope slide and a 170 µm thick glass coverslip, sealed with a 120 µm thick spacer.
The optical trap, described in details in Fig. 5, is an evolution of the setup described in our previous work [15][16][17]. It uses a CW near-infrared (λ 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.
FIG. 5. The trapping laser (λ T = 785 nm, 100 mW, TEM 00 , CW, Coherent, OBIS LX785) is modulated externally using a waveform generator (Agilent, 33220A). Linearly polarized along the z−axis, the beam is sent to a water-immersion objective (O1, 100×, 1.2 numerical aperture (NA)) through a polarizing beam splitter (PBS) and a quarter-wave plate (λ /4). The intensity I(t) partially reflected by the end-surface of the fluidic cell varies linearly with the displacement x(t) of the polystyrene microsphere inside the trap. This intensity I(t) is collected and recorded by a p-i-n photodiode (Thorlabs, DET10A), while a CCD camera is used in the other port of the non polarizing beam splitter (NPBS) for imaging. The probe beam consists of a second laser (639 nm, 70 mW Thorlabs laser diode, linearly polarized) of low power (400 µW). It is injected inside the trap collinearly with the trapping beam but from behind the fluidic cell using a dry objective (O2, 60×, NA 0.7). This second beam is separated from the trapping beam using a dichroic mirror (DM) and the interference between the transmitted beam and the diffracted light by the bead is recorded using a second p-i-n photodiode (Thorlabs, DET10A) placed in a plane conjugated to the back focal plane of the trapping objective. In order to ensure that a single bead is trapped without other beads in its vicinity, potentially perturbing the dynamics, the optical trap is equipped with an interferometric scattering microscope not shown here but described in details in our previous work [17].
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 = 2 18 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).

Stiffness modulation calibration
The trapping laser is modulated according to a given protocol κ(t), defined and calculated with chosen transition parameters (κ i , κ f , ∆t). In order to convert this protocol κ(t) into a modulating voltage V 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.
FIG. 6. Evolution of the trap stiffnesses as a function of DC waveform generator voltages. The red dots represent the stiffness values extracted from the motional PSD, with error bars for each point combining the uncertainties of the Lorentzian fit of each PSD and the error made on the Stokes drag γ = 6πRη due to the polystyrene sphere radius dispersion δ R/R. The solid line is the linear fit and the shaded area represents a 95 % confidence interval for the estimated linear regression parameters taking into account the weights of the data points.

Monitoring Brownian dynamics
The time evolution of the Brownian system is monitored by recording the stochastic trajectory of the trapped bead over 2 × 10 4 cycles of the protocol κ(t). Each cycle lasts 50 ms, where the first 30 ms correspond to the initial thermal equilibrium with κ i and the remaining (20 − ∆t) ms correspond to the final thermal equilibrium at κ f . Each stationary region of the full trajectory, i.e. corresponding to a constant κ (κ i or κ 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 κ. Figs. 7 (a) and (b) respectively show the PSD of the concatenated trajectories for the equilibria κ i and κ f for the case ∆t ∼ τ relax /10 described in the main text. FIG. 7. The power spectral density of the concatenated trajectories corresponding to the sections of the cycles for which κ is fixed to κ i is displayed in greeen. The best-fitted roll-off frequency f c = 52.63 ± 0.01 Hz (vertical red line) yields κ i = 2.78 ± 0.08 pN/µm, and the position sensitivity parameter is β = k B T /γD f it = 1.21 ± 0.02 µm/V -see below. The blue curve is the power spectral density of the concatenated trajectories corresponding to the sections of the cycles for which κ is fixed to κ f . The best-fitted rolloff frequency is f c = 98.98 ± 0.02 Hz (vertical purple line) gives κ f = 5.22 ± 0.15 pN/µmfor this case. Here, the positional calibration factor is β = 1.31 ± 0.02 µm/V. Lorentzian fits (continuous red and purple lines superimposed to the PSDs) are calculated by implementing a MATLAB Levenberg-Marquardt algorithm for non-linear leasts squares.
Implementing the same procedure, the full temporal trace of the particle positions undergoing 2 × 10 4 cycles is chopped into trajectories that correspond to a single cycle of the protocol κ(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 κ i , as displayed in Fig. 8 below.
The instantaneous ensemble variance s(t j ) at a time t = t j ( j = 1, · · · , T × 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 ρ(x,t j ) is a Gaussian of zero mean µ x (t j ) and variance s(t j ). Fig. 9 displays the position distribution functions (PDF) before (equilibrium at κ i ) and after (equilibrium at κ f ) the change in trapping stiffness imposed by the protocol κ(t). The corresponding trapping potentials calculated as U(x,t j ) = −k B T log(ρ(x,t j )) + cst are also shown and compared to the expected harmonic profiles U = 1 2 κx 2 evaluated from the stiffnesses κ i , κ 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 κ(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.

Statistical uncertainties
The uncertainties for the instantaneous ensemble variances are obtained following a χ 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 κ(t).

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 κ). The fit parameter D fit is compared to the diffusion coefficient D = k B T /γ expected from the Fluctuation-Dissipation Theorem, assuming known temperature and viscosity. This gives a conversion factor β = D/D 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 δ R/R = 2.8% size dispersion deviation of the trapped beads.
Instantaneous positions are thus given from the conversion factor as x(t) = (β ± δ β )V (t), and therefore the variance, up to first-order in uncertainty, x 2 (t) = (β 2 ± 2β δ β )V 2 (t), (since µ x (t) = 0). The total error of the variance writes as: where σ 2 is the estimator of the instantaneous ensemble variance over N cycles, δ σ 2 χ 2 corresponds to the statistical uncertainty in the motional variance determination (see above) and δ β σ 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 ∆t as the interval over which κ(t) remains constant (either at κ i or κ f ), the temporal average of the corresponding variance is: taking ∆t as the interval over which κ(t) remains constant (either at κ i or κ f ) and n = ∆t · f s with f s = 2 18 Hz, the sampling frequency. The standard deviation of the temporal average is simply evaluated as: The stationary variances s i and s f and their uncertainties are thus simply given by: where δ s t = 1/∆t ∑ n j=1 δ s(t j ).

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.

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 κ i to a final one κ f with, for all protocols, fixed and identical κ f , κ 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. [6]. We calculate κ ESE (t) following [6] for a transfer duration of ∆t = 3.47 × 10 −4 s. Over the same transfer duration, we also implement our optimal protocol κ opt (t). All protocols are displayed in Fig. 10. FIG. 10. Calibrated signal of the function generator, for a step-like (green), ESE (pink), and optimal (blue) protocols. The stiffness κ(t) is normalized to the initial stiffness κ i . The jump for the transition κ i → κ f starts at t 0 = 0 s and, for the case of ESE and optimal ends at ∆t = 3.47 × 10 −4 s, with κ i = 2.77 ± 0.08, κ f = 5.22 ± 0.15 pNµm. The ESE protocol κ ESE was computed based on Eq. (8) in [6]. Fig. 11 gathers the time evolutions of the motional variances associated with each protocol. As expected, the steplike 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. FIG. 11. Temporal evolution of the variance s(t), after t 0 = 0 s for the step-like protocol (in green), and the ESE (in purple) and optimal (in blue) protocols. The variances are normalized to the final equilibrated variance s f . The data points represent ensemble mean values of the variance s(t) for each protocol. The shaded areas show the respective 95% confidence intervals. Both ESE and optimal protocols reach an equilibrium regime s f at ∆t = 3.47 × 10 −4 s ∼ τ relax /10 by construction. Inset: The control parameterκ(s) as a function of the variance s, with the same color codes as in the main figure.

APPENDIX C: SMOOTH PROTOCOLS
The optimal protocol obtained in this work [Eq. (6) in the main text] was derived using the Lagrangian density A peculiar feature of L[s,κ(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κ(s) (i.e. s i κ i = s f κ f = Dγ) 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κ(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κ(s) by adding a further term to the Lagrangian density (16), which becomes: FIG. 12. Temporal evolution of the mean cumulative energetics of the different protocols, step-like (lower inset), ESE and optimal. The mean cumulative work for the optimal protocol is displayed in blue, with total work W (t) opt = 0.981 ± 0.059 k B T . The mean cumulative work for the ESE protocol is displayed in pink, with total work W (t) ESE = 1.142 ± 0.075 k B T . The mean cumulative heat generated through the optimal protocol is displayed in orange and the ESE protocol in yellow. Both are superimposed to the work, with total heat Q(t) opt = 0.983 ± 0.060 k B T and Q(t) ESE = 1.142 ± 0.076 k B T . Shaded areas represent 95% confidence levels. Lower inset: Energetics for the step-like protocol. As expected, the mean cumulative work (in green) reaches immediately W step = 0.45 ± 0.04 k B T . In brown, the heat, in contrast, achieves the equilibrium value W = Q with Q step = 0.45 ± 0.04 k B T only after τ relax . Upper inset: Comparison between the excess work values of the ESE protocol (pink) and the optimal one (blue) for the transfer duration of duration ∆t = t f = 3.47 × 10 −4 s. The non-optimal character of the ESE protocol is directly measured with ∆W ESE = 0.81 ± 0.08 k B T larger than the optimal value ∆W opt = 0.65 ± 0.07 k B T . The universal bound ∆W = γ( √ s i − √ s f ) 2 /∆t discussed in the main text is shown by the continuous line.
where ε is an additional Lagrange multiplier. The above Lagrangian density yields the Euler-Lagrange equation (11) in the main text, which we reproduce here: As a second-order differential equation, Eq. (18) needs two independent boundary conditions, thus enabling us to set s i κ i = s f κ f = Dγ, as requested for our protocols. When ε → 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κ(s) at s i and s f . Through the Lagrange multiplier ε, one can limit the value of such derivative, so that the protocol becomes smoother and smoother as ε increases. Equation (18) can be solved numerically by successive iterations. We used the following scheme: where the superscript n denotes the n-th iteration, while the subscript i refers to the discrete grid s i = i δ s, with spacing equal to ds. The second derivative is then approximated with the standard finite-difference formula: The parameter α > 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κ n+1 i =κ n i ). As an example, we have solved Eq. (18) with physical parameters D = γ = 1 and λ = 0.81, corresponding to a total duration for the optimal protocol ∆t opt ∼ τ relax /6 according to Eq. (9) in the main text. The boundary values are s i = 1 and s f = 0.5, κ i = 1 and κ f = 2. The smoothness parameter is ε = 10 −5 . The numerical convergence parameter is set to α = 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 × τ relax , longer than that of the optimal one. But the total work is smaller W smooth = 1.32 < W opt = 1.38. The time-energy product is (∆t ∆W ) smooth = 0.356 > (∆t ∆W ) opt = 0.343, in agreement with the theoretical considerations detailed in the main text.