Optimal work in a harmonic trap with bounded stiffness

TL;DR

This paper uses optimal control theory to determine the minimal work required to transition a Brownian particle between states in a harmonic trap with bounded stiffness, revealing new solution types and practical implications.

Contribution

It extends previous models by incorporating bounded stiffness constraints, showing how they affect optimal protocols and work minimization in thermodynamic processes.

Findings

Bounded stiffness can prevent certain state transitions.

Different solution regimes depend on operation time and compression ratio.

Work minimization is significantly impacted by stiffness bounds.

Abstract

We apply Pontryagin's principle to drive rapidly a trapped overdamped Brownian particle in contact with a thermal bath between two equilibrium states corresponding to different trap stiffness $\kappa$. We work out the optimal time dependence $\kappa(t)$ by minimising the work performed on the particle under the non-holonomic constraint $0\leq\kappa\leq\kappa_{\max}$, an experimentally relevant situation. Several important differences arise, as compared with the case of unbounded stiffness that has been analysed in the literature. First, two arbitrary equilibrium states may not always be connected. Second, depending on the operating time $t_{\text{f}}$ and the desired compression ratio $\kappa_{\text{f}}/\kappa_{\text{\i}}$, different types of solutions emerge. Finally, the differences in the minimum value of the work brought about by the bounds may become quite large, which may have a relevant impact on the optimisation of heat engines.

Equations

$$\frac{d⟨x^2(t)⟩}{dt}=\phi (⟨x^2(t)⟩,\kappa (t)),$$

$$W=\frac{1}{2}{\int }_0^{t_{\text{f}}}dt⟨x^2⟩\dot{\kappa }(t)=-\frac{1}{2}{\int }_0^{t_{\text{f}}}dt\kappa (t)\phi (⟨x^2(t)⟩,\kappa (t)),$$

$$L(⟨x^2⟩,\kappa )=-\frac{\kappa }{2}\phi (⟨x^2⟩,\kappa ),$$

$$W={\int }_0^{t_{\text{f}}}dtL(⟨x^2(t)⟩,\kappa (t)).$$

$$0\le \kappa (t)\le {\kappa }_{\max }.$$

$$\gamma \frac{dx(t)}{dt}=-\kappa (t)x(t)+\xi (t),$$

$$⟨\xi (t)⟩=0,⟨\xi (t)\xi (t^{\prime })⟩=2D\delta (t-t^{\prime }),$$

$$\hat{\kappa }=\frac{\kappa }{{\kappa }_{\text{i}}},\hat{t}=\frac{{\kappa }_{\text{i}}}{\gamma }t,{\hat{y}}_1=\sqrt{\frac{⟨x^2⟩}{⟨x^2{⟩}_{\text{i}}}},$$

$$\frac{d{y}_1(t)}{dt}=f_1(y_1(t),\kappa (t)),$$

$$f_1(y_1,\kappa )\equiv \frac{1}{y_1}-\kappa y_1,$$

$$W=\frac{1}{2}\ln {\kappa }_{\text{f}}+{\int }_0^{t_{\text{f}}}dt{\left[f_1(y_1(t),\kappa (t))\right]}^2,$$

$$\kappa (0)={\kappa }_{\text{i}}=1,\kappa (t_{\text{f}})={\kappa }_{\text{f}},$$

$$y_1(0)\equiv y_{1,\text{i}}=1,y_1(t_{\text{f}})\equiv y_{1,\text{f}}=\frac{1}{\sqrt{{\kappa }_{\text{f}}}}.$$

$$W_{\text{irr}}[y_1,\kappa ]={\int }_0^{t_{\text{f}}}dt{f}_0(y_1(t),\kappa (t)),$$

$$f_0(y_1,\kappa )\equiv {\left[f_1(y_1,\kappa )\right]}^2,$$

$$t_{\text{f}}<t_{\text{eq}},t_{\text{eq}}\simeq \frac{3}{{\kappa }_{\text{f}}},$$

$$\frac{d{y}_0}{dt}=f_0(y_1(t),\kappa (t))={\left[f_1(y_1(t),\kappa (t))\right]}^2.$$

$$\Pi (y_1,{\psi }_0,{\psi }_1,\kappa )$$

$$\mathcal{H}(y_1,{\psi }_0,{\psi }_1)$$

$$\frac{d{y}_i}{dt}=\frac{\partial \Pi }{\partial {\psi }_i},\frac{d{\psi }_i}{dt}=-\frac{\partial \Pi }{\partial y_i},i=0,1,$$

$$\frac{d{\psi }_0}{dt}=0\Rightarrow {\psi }_0(t)={\psi }_0\text{ (constant),}$$

$$\frac{d{\psi }_1}{dt}=-{\psi }_0\frac{\partial f_0}{\partial y_1}-{\psi }_1\frac{\partial f_1}{\partial y_1}=-\frac{\partial f_1}{\partial y_1}\left(2{\psi }_0{f}_1+{\psi }_1\right).$$

$$\Pi (y_1^{\ast }(t),{\psi }_0^{\ast }(t),{\psi }_1^{\ast }(t),{\kappa }^{\ast }(t))=\mathcal{H}(y_1^{\ast }(t),{\psi }_0^{\ast }(t),{\psi }_1^{\ast }(t)).$$

$$W_{\text{irr}}^{\min }=y_0^{\ast }(t_{\text{f}}).$$

$$0={\frac{\partial \Pi }{\partial \kappa }|}_{\tilde{\kappa }}$$

$$\frac{d{\psi }_1}{dt}=\left(\frac{1}{y_1^2}+\kappa \right)\left[2{\psi }_0\left(\frac{1}{y_1}-\kappa y_1\right)+{\psi }_1\right],$$

$$0={\frac{\partial \Pi }{\partial \kappa }|}_{\tilde{\kappa }}=-y_1\left[2{\psi }_0\left(\frac{1}{y_1}-\tilde{\kappa }{y}_1\right)+{\psi }_1\right].$$

$$\tilde{\kappa }=\frac{{\psi }_1}{2{\psi }_0{y}_1}+\frac{1}{y_1^2}.$$

$$\frac{d{y}_0}{dt}={\left(\frac{{\psi }_1}{2{\psi }_0}\right)}^2,$$

$$\frac{d{\psi }_0}{dt}=0,$$

$$\frac{{\psi }_{1,\text{u}}}{2{\psi }_{0,\text{u}}}=-\frac{y_{1,\text{f}}-1}{t_{\text{f}}},$$

$$y_{1,\text{u}}(t)=1+\frac{y_{1,\text{f}}-1}{t_{\text{f}}}t.$$

$$y_{0,\text{u}}(t)=\frac{(y_{1,\text{f}}-1{)}^2}{t_{\text{f}}^2}t.$$

Full text

I Introduction

One of the key parameters in non-equilibrium transformations is the characteristic relaxation time of the system under study. In general, equilibrium states of a system depend on the values of certain physical properties $\lambda$ that can be externally controlled, such as the available volume for a gas or the spring constant of the harmonic potential that confines a colloidal particle. When one relevant external parameter is abruptly changed from $\lambda_{\text{i}}$ to $\lambda_{\text{f}}$, a system that was at the equilibrium state corresponding to $\lambda_{i}$ begins to evolve and, as time increases, approaches the new equilibrium state corresponding to $\lambda_{\text{f}}$. The system’s equilibration time $t_{\text{eq}}$ can be loosely defined as the time that the system needs to reach the new equilibrium configuration, and it is an intrinsic property for each physical system that depends on the underlying interactions, encoded in the transport coefficients, the external parameters $\lambda$, and the temperature.

Recently, there has been a growing interest in the development of engineered techniques capable of beating the natural time scale for relaxation between equilibrium states. Inspired by the so-called shortcut to adiabaticity processes Chen et al. (2010a, b), specific procedures that make it possible to connect equilibrium states using linking times much shorter than the natural equilibration time have been devised. The term Engineered Swift Equilibration (ESE) has been coined to describe these kind of procedures. The general idea of an ESE process is to design a tailor-made time dependent protocol $\lambda(t)$ for the externally controlled parameter, such that the system is driven from the equilibrium state corresponding to $\lambda_{\text{i}}$ to the equilibrium state corresponding to $\lambda_{\text{f}}$ in a finite time $t_{\text{f}}$, ideally much shorter than the equilibration time $t_{\text{eq}}$. Such protocols have been established for an isolated dilute gas confined in a 3D isotropic harmonic trap Guéry-Odelin et al. (2014), and for nano systems in contact with a thermostat both in the over or underdamped regime Martínez et al. (2016a); Le Cunuder et al. (2016).

Here, we focus on a colloidal particle confined by a harmonic trap of stiffness $\kappa$ Martínez et al. (2016a); Le Cunuder et al. (2016); Rondin et al. (2017); Chupeau et al. (2018a, b), the relevant physical quantity is the variance $\langle x^{2}\rangle$ of its position. Initially, the stiffness of the trap is $\kappa_{\text{i}}$ and the particle is equilibrated at the temperature $T$ of the fluid in which it is immersed, $\langle x^{2}\rangle_{\text{i}}=k_{B}T/\kappa_{\text{i}}$, $k_{B}$ being Boltzmann’s constant. Throughout this work, we consider processes in which the temperature of the bath is kept constant, at difference with the approach in Chupeau et al. (2018b). In a STEP process, the stiffness of the trap is suddenly changed to a different value $\kappa_{\text{f}}$ at $t=0^{+}$, and the relaxation of the colloidal particle to the new equilibrium state is tracked. Basically, its variance $\langle x^{2}\rangle$ relaxes exponentially to its new equilibrium value $\langle x^{2}\rangle_{\text{f}}=k_{B}T/\kappa_{\text{f}}$ after a characteristic time $t_{\text{eq}}\simeq 3k_{B}T/(\kappa D)$, where $D$ is the diffusion coefficient. Alternatively, the system can be compressed/decompressed isothermally by introducing a suitable time protocol $\kappa(t)$ for the stiffness that drives the system from the initial equilibrium state with $\kappa_{\text{i}}$ to the final equilibrium state with $\kappa_{\text{f}}$ in a finite time $t_{\text{f}}$. The ESE procedure consists of choosing in a smart way the stiffness protocol $\kappa(t)$, so that $t_{\text{f}}\ll t_{\text{eq}}$, thus beating the system’s natural rate of equilibration. For example, the protocol employed in Ref. Martínez et al. (2016a) beats the natural relaxation time by two orders of magnitude, $t_{\text{f}}/t_{\text{eq}}\simeq 0.01$.

Once it has been shown that ESE processes are indeed possible, an optimisation problem arises. There is a wide class of functions $\lambda(t)$ that connect the initial and final equilibrium states in a given time $t_{\text{f}}$. For each of the possible functions $\lambda(t)$, one can calculate the work performed in the process $W=\int_{0}^{t_{\text{f}}}\frac{\partial{H}}{\partial{\lambda}}\dot{\lambda}\,dt$, where $H$ is the Hamiltonian of the system Sekimoto (2010); mathematically, $W$ is a functional of $\lambda$. Hence the question, for a given connection time $t_{\text{f}}$: what is the optimal time evolution $\lambda^{*}(t)$ that minimises (on average) the work $W$?

For the colloidal particle in a harmonic trap, the optimal time evolution for the stiffness $\kappa^{*}(t)$ has been obtained for different boundary conditions Schmiedl and Seifert (2007, 2008). The specific boundary conditions that are adequate for the ESE process were considered in Schmiedl and Seifert (2008), in the context of building a stochastic heat engine. This result has also been rederived in later works, see for example Aurell et al. (2011). The optimal protocol for the stiffness has finite jumps both at the initial and final times, $\kappa^{*}(t=0^{+})\neq\kappa_{\text{i}}$ and $\kappa^{*}(t_{\text{f}}^{-})\neq\kappa_{\text{f}}$. This kind of discontinuity at the endpoints of the time interval is usual in stochastic thermodynamics and stem from the “Lagrangian” of the considered variational problem being linear in the “velocities” Band et al. (1982), which is sometimes known as the Miele problem Tolle (2012). This discontinuities can be regularised by introducing an additional small term in the Lagrangian, which introduces two boundary layers of finite width at the endpoints of the time interval that eliminate the finite jumps Aurell et al. (2012); Muratore-Ginanneschi and Schwieger (2017).

The main shortcoming of previous protocols, be they optimal or not, comes about in decompression processes. Any protocol involving a short enough time $t_{\text{f}}$ entails that the stiffness has to be transiently negative inside a certain time window Chupeau et al. (2018a, b), similarly to the situation found in other systems Torrontegui et al. (2013); Guéry-Odelin et al. (2014). The arising of negative values for the stiffness is challenging from an experimental point of view, since the potential should change from confining to repulsive. In the usual experimental setups, the stiffness $\kappa$ of the harmonic potential is always positive and, in addition, has a certain upper bound $\kappa_{\max}$ depending on the technique employed—mainly atomic force microscopy (AFM) or laser optical tweezers (LOT) Ritort (2006); Wen et al. (2007); Manosas et al. (2007); Hoffmann and Dougan (2012); Marszalek and Dufrêne (2012); Rondin et al. (2017); Ciliberto (2017)—to implement the harmonic trap. The existence of this upper limit is related to the validity of the harmonic approximation. The intrinsic limit of ESE protocols is dictated by the accuracy of the mathematical model that describes the physical system.

In light of the above remarks, it is relevant to investigate the optimisation problem of the work described above when the stiffness of the trap is restricted to a certain interval, $\kappa\in[0,\kappa_{\max}]$. The existence of an upper bound also changes the problem, since very high compression ratios $\kappa(t)/\kappa_{\text{i}}$ have to be applied to accelerate the equilibration in the compression case. For example, in Ref. Martínez et al. (2016a), transient compression ratios of the order of $40$ were applied in order to speed up the equilibration of the particle, even when $\kappa_{\text{f}}$ only doubled $\kappa_{\text{i}}$.

These drawbacks are important for the optimisation of irreversible heat engines, a field of research that has become quite active in the last few years Esposito et al. (2010); Roßnagel et al. (2014); Martínez et al. (2016b, 2017); Taye (2017); Apertet et al. (2017). In fact, Brownian particles trapped by optical tweezers have been recently employed to build stochastic heat engines, both theoretically and experimentally Schmiedl and Seifert (2008); Blickle and Bechinger (2012); Martínez et al. (2016b), for a review see Ciliberto (2017). In these studies, the stiffness of the trap is changed as a function of time by tuning the laser power, and decreasing (resp. increasing) the stiffness is equivalent to decompressing (resp. compressing) the system. Cyclic engines are thus built by connecting isothermal compression/decompression branches with either isochoric Blickle and Bechinger (2012) or isoentropic branches Schmiedl and Seifert (2008); Blickle and Bechinger (2012). In the decompression (resp. compression) branch the corresponding work $W_{d}$ (resp. $W_{c}$) is negative (resp. positive), and the total work $W=W_{c}+W_{d}$ must be negative to build a heat engine.

In this work, we focus on the analysis of isothermal compression/decompression processes, i.e. the isothermal branches of the heat engines described in the previous paragraph. Note that the optimisation of the work considered here is relevant in the context of heat engines, since the extracted work $-W$ has to be a maximum, i.e $W$ must be minimum Schmiedl and Seifert (2008). In addition, the stiffness is restricted in experiments to a certain interval as explained above, and thus the externally controlled function $\kappa(t)$ obeys the non-holonomic constraint $0\leq\kappa\leq\kappa_{\max}$. Therefore, the currently available “unconstrained” results Schmiedl and Seifert (2008); Aurell et al. (2011) are not useful for short enough times $t_{\text{f}}$, because the optimal $\kappa(t)$ becomes negative (resp. larger than $\kappa_{\max}$) in decompression (resp. compression) processes.

The time evolution of the colloidal particle is governed by a first-order differential equation,

[TABLE]

where $\varphi$ is a smooth function of both $\langle x^{2}\rangle$ and $\kappa$, see for example Schmiedl and Seifert (2008); Martínez et al. (2016a). Then, $\kappa(t)$ is a control function, in the sense used in control theory. The mean work in a finite time isothermal process can be written as

[TABLE]

where we have made use of the relation $\kappa_{\text{i}}\langle x^{2}\rangle_{\text{i}}=\kappa_{\text{f}}\langle x^{2}\rangle_{\text{f}}=k_{B}T$. By defining

[TABLE]

we can write

[TABLE]

We then have a well-posed problem in control theory Pontryagin (1987); Liberzon (2012). We seek the minimum of $W$, taking into account that the evolution of $\langle x^{2}\rangle$ is controlled by $\kappa$, as described by (1), where $\kappa(t)$ satisfies the non-holonomic constraint

[TABLE]

This kind of optimisation problem cannot be tackled with the usual tools of variational calculus, i.e. the Euler-Lagrange equations; they must be addressed by applying more sophisticated tools from control theory, such as Pontryagin’s maximum principle Pontryagin (1987); Liberzon (2012).

The plan of the paper is as follows. Section II is devoted to the statement of the minimisation of the work as a control problem. Therein, we explain how Pontryagin’s principle can be applied to this particular situation. In Sec. III, we address the minimisation problem when the stiffness is not bounded and can thus have any value, including negative ones. Next, we look into the minimisation problem with bounds in Sec. IV, first for the decompression case in IV.1 and afterwards for the compression case in IV.2. Section V discusses the different phases that appear in the minimisation problem and a detailed comparison between the values of the optimal work for the unbounded and the bounded cases is carried out. The main conclusions are presented in Sec. VI. Finally, the appendices deal with some technicalities that are omitted in the main text.

II The control problem

II.1 Statement

We consider a colloidal particle immersed in a fluid at temperature $T$. The particle is in a harmonic trap of stiffness $\kappa(t)$, the time dependence of which is externally controlled, and we are interested in time scales such that the overdamped limit holds. Thus, the dynamics of the particle position $x$ is governed by the Langevin equation

[TABLE]

where $\gamma$ is the friction coefficient and $\xi(t)$ is a Gaussian white noise force,

[TABLE]

in which $D$ is the diffusion coefficient that is connected to $\gamma$ by the fluctuation-dissipation relation $D=k_{B}T/\gamma$. Implicitly, our modelling assumes that the relaxation of the surrounding fluid to equilibrium can be regarded as instantaneous on the time scale over which the stiffness varies.

The Fokker-Planck equation associated to the Langevin equation (6) is linear. Therefore, in the class of ESE processes described in the introduction, the probability distribution function $\rho(x,t)$ is Gaussian for all times, since it is so initially, and we can characterise the stochastic process completely by its variance $\langle x^{2}(t)\rangle$. To do the calculations, it is convenient to introduce dimensionless variables

[TABLE]

where the initial value of the variance is $\langle x^{2}\rangle_{\text{i}}=k_{B}T/\kappa_{\text{i}}$. Therefore, $\hat{y}_{1}(t)$ is the non-dimensional standard deviation. In order not to clutter our formulas, we omit the hats in the dimensionless variables henceforth.

The time evolution of the standard deviation $y_{1}$ is governed by the first-order differential equation

[TABLE]

for each given time-dependent stiffness $\kappa(t)$.

The mean work performed on the system is defined at the average level as $\,{\mathchar 22\relax\mkern-12.0mud}W=\frac{1}{2}\langle x^{2}\rangle\,d\kappa$ Sekimoto (2010), which is positive when energy is transferred from the environment to the particle and negative otherwise. The unit of energy is $k_{B}T$, then the dimensionless work for a finite transformation from $t=0$ to $t=t_{\text{f}}$ is, after using integration by parts Schmiedl and Seifert (2007)

[TABLE]

where we have made use of the boundary conditions for our ESE problem

[TABLE]

The first term on the rhs of (10) is the free energy difference between the initial and final states. Then, the second term on the rhs, which is non-negative, is the irreversible work and vanishes only in the quasi-static limit, when $t_{\text{f}}\to\infty$ Schmiedl and Seifert (2007, 2008).

Here, we are interested in minimising $W$ (i.e. maximising the “extracted” work $-W$) for a fixed time interval $t_{\text{f}}$, starting from the equilibrium state corresponding to $\kappa_{\text{i}}$, equal to unity in dimensionless variables, and ending up in the equilibrium state corresponding to $\kappa_{\text{f}}$. Therefore, we have to minimise the irreversible work as given by the functional

[TABLE]

where the stiffness of the trap $\kappa(t)$ is an externally controlled function and the time evolution of $y_{1}(t)$ is linked thereto by (9a). For the ESE processes, we are especially interested in the regime

[TABLE]

where $t_{\text{eq}}$ is the equilibration time when the system relaxes to equilibrium with time-independent stiffness, $\kappa(t)=\kappa_{\text{f}}$ for all times Martínez et al. (2016a).

Let us be more specific. For each time-dependent control function $\kappa(t)$, we obtain a certain time evolution for $y_{1}(t)$ by integrating (9a), and therefore a certain value for our functional $W_{\text{irr}}$. What we are interested in is finding out whether there is an optimal control function $\kappa^{*}(t)$, for which the corresponding time evolution of the standard deviation is $y_{1}^{*}(t)$, such that $W_{\text{irr}}[y_{1},\kappa]\geq W_{\text{irr}}^{*}\equiv W_{\text{irr}}[y_{1}^{*},\kappa^{*}]$ within a certain class $\mathcal{K}$ of admissible control functions. From a physical point of view, it is reasonable to admit functions $\kappa(t)$ with finite instantaneous jumps at certain times $t\in[0,t_{\text{f}}]$; therefore we assume that $\kappa(t)$ is piecewise continuous in $[0,t_{\text{f}}]$. Note that this entails that $y_{1}(t)$ must be continuous in $[0,t_{\text{f}}]$ since Eq. (9a) implies that $\dot{y}_{1}$ has at most finite jump discontinuities.

The boundary conditions for our minimisation problem stem from the ESE process we are interested in, and are given by (11a). At this point, we have a well-posed optimal control problem Pontryagin (1987); Gelfand and Fomin (2000); Liberzon (2012). We want to minimise the functional (12), in which the time evolution of $y_{1}(t)$ is controlled by the imposed program $\kappa(t)$ by means of the evolution equation (9a), with the boundary conditions for $y_{1}$ given by (11b). This minimisation is done over the class of admissible controls: piecewise continuous functions $\kappa(t)$ that verify the prescribed boundary conditions for $\kappa$, as given by (11a). In addition, we may have more restrictions on $\kappa$, which we summarise here by saying that the possible values of the control $\kappa(t)\in{U}$. The so-called control set $U$ is a certain subset (interval) of the real numbers, ${U}\subseteq\mathbb{R}$. Although our notation does not make it explicit, the control set ${U}$ can vary in time, see for example section $3.3$ of Liberzon (2012).

II.2 Pontryagin’s procedure

The solution to this control problem is obtained by applying Pontryagin’s principle, see section $1.8$ of Pontryagin (1987) or section $4.3.1$ of Liberzon (2012) for its general formulation. Below, we explain how Pontryagin’s maximum principle is applied to our particular physical situation.

First, we define a variable $y_{0}$ such that $y_{0}(0)=0$ and

[TABLE]

It is clear that, for each choice of the control function $\kappa(t)$, $y_{0}(t_{\text{f}})$ equals the value of the functional $W_{\text{irr}}[y,\kappa]$. Next, we introduce variables $\psi_{i}$ conjugate to each $y_{i}$, $i=0,1$, and define a function

[TABLE]

Note that, by construction, $\Pi$ does not depend on $y_{0}$. For fixed $(y_{1},\psi_{0},\psi_{1})$, the function $\Pi$ becomes a function of $\kappa$, which belongs to the control set, $\kappa\in{U}$. We denote the supremum of this function by $\mathcal{H}$,

[TABLE]

In conjunction with (15), the following system of equations hold for the variables $(y_{0},y_{1},\psi_{0},\psi_{1})$

[TABLE]

i.e. we recover (14) and (9a) for the evolution of $(y_{0},y_{1})$ and obtain the evolution equations for the conjugate variables $(\psi_{0},\psi_{1})$

[TABLE]

For any control function $\kappa(t)$ linking $y_{1,\text{i}}$ and $y_{1,\text{f}}$ in a time $t_{\text{f}}$, we have a solution $y_{1}(t)$ of (9a). Inserting both $\kappa(t)$ and the associated $y_{1}(t)$ into (18), we also obtain the solutions for the conjugate variables $(\psi_{0},\psi_{1}(t))$ associated to the considered control. This construction defines the conjugate variables, and consequently the function $\Pi$.

Pontryagin’s extremum principle states a necessary condition for having an optimal control $\kappa^{*}(t)$ that minimises the functional $W$, within the considered class of admissible controls. Let $\kappa^{*}(t)$ be an admissible control and $y_{1}^{*}(t)$ the associated solution of (9a). In order that $\kappa^{*}(t)$ yield a solution of the minimisation problem, there must exist a solution of (18) $(\psi_{0}^{*},\psi_{1}^{*}(t))\neq(0,0)$ for all $t\in[0,t_{\text{f}}]$ such that

for all $t\in[0,t_{\text{f}}]$, it is at the point $\kappa=\kappa^{*}(t)$ that the function $\Pi(y_{1}^{*}(t),\psi_{0}^{*}(t),\psi_{1}^{*}(t),\kappa)$ attains its maximum, i.e.

[TABLE] 2. 2.

The constant $\psi_{0}^{*}\leq 0$.

The latter condition assures that $\Pi$ has a maximum at $\kappa^{*}$ 111The main point is that the sign of all the momenta $\psi_{k}$ and thus the sign of $\Pi$ can be reversed, which gives a “mirrored” solution of the canonical equations. Over this “mirrored” solution, with $(-\psi_{0})>0$, the corresponding $(-\Pi)$ would reach an infimum at $(-\mathcal{H})$, instead of a supremum. It is to fix this ambiguity in Pontryagin’s procedure and formulate a maximum principle that the choice $\psi_{0}<0$ is made Liberzon (2012).. The idea behind Pontryagin’s principle is to rewrite the functional to be extremalised as $\int_{0}^{t_{\text{f}}}dt\,\psi_{0}f_{0}=\int_{0}^{t_{\text{f}}}dt\,(\Pi-\psi_{1}f_{1})$. Taking advantage of the Hamiltonian structure behind (17) yields the formalism in question.

From the optimal control, one deduces the corresponding $y_{0}^{*}(t)$ and the minimum irreversible work is

[TABLE]

Finally, it is straightforward to show that $\mathcal{H}(y_{1}^{*}(t),\psi_{0}^{*}(t),\psi_{1}^{*}(t))$ does not depend on time, i.e. it is a constant of motion.

At this point, the issue is finding the supremum of the function $\Pi(y_{1},\psi_{0},\psi_{1},\kappa)$ that leads to the optimal control $\kappa^{*}(t)$. The basic idea is that, for any time $t$, the value of the optimal control $\kappa$ can lie either inside $U$ or along its boundary $\partial U$. This is completely analogous to the situation found when seeking an extremum of a function of several variables $g(x_{1},x_{2},\ldots,x_{N})$ in a certain closed subset ${U}\subset\mathbb{R}^{N}$, which may lie inside ${U}$ or on its boundary $\partial{U}$. To find it, first we look for the extremum $(x_{1}^{*},x_{2}^{*},\ldots,x_{N}^{*})$ by imposing $(\partial g/\partial x_{i})^{*}=0$; if this equation does not have a solution inside ${U}$, the extremum must lie on the boundary $\partial{U}$. Therefore, to obtain the supremum of $\Pi$, at first $\kappa^{*}$ is sought by writing

[TABLE]

We have introduced the notation $\tilde{\kappa}$ to make it clear that $\tilde{\kappa}$ may be the “right” solution, i.e. $\tilde{\kappa}=\kappa^{*}$, or not. Being more concrete, there appear two possibilities:

The specific $\tilde{\kappa}$ found from (20) belongs to the class of admissible controls for all times $t$, then we have found the solution of the minimisation problem, $\kappa^{*}=\tilde{\kappa}$. 2. 2.

$\tilde{\kappa}$ does not belong to the class of admissible controls because at a certain time $t_{0}<t_{\text{f}}$ we have that $\tilde{\kappa}(t_{0})$ lies outside the control set $U$. Then, the optimal $\kappa^{*}(t)$ comprises in general several branches: some branches stem from (20) and lie inside $U$ whereas other branches lie over its boundary $\partial{U}$.

Now we derive some specific expressions for our system. First, we write the particular evolution equation for the conjugate variable $\psi_{1}$,

[TABLE]

where we have taken into account the definition of $f_{1}(y,\kappa)$ in (9b). Second, we derive the particular equation for $\tilde{\kappa}$. Making use of (20) and the definition of $f_{1}(y,\kappa)$,

[TABLE]

and thus

[TABLE]

The insertion of (23) into the set of differential equations (17) yields

[TABLE]

In the following sections, we analyse in depth two particular cases: (i) when the stiffness may have any value including negative ones, see Section III, and (ii) when the stiffness is bounded and lies within a certain interval $[0,\kappa_{\max}]$, see Section IV. Note that the latter is the relevant problem at the experimental level, as explained in the introduction.

III Unbounded stiffness

First, we consider the simplest situation: we have no other restrictions on the control function $\kappa(t)$ aside from the boundary conditions (11a). Therefore, the class of admissible control functions $\mathcal{K}$ comprises all piecewise continuous functions lying inside the vertical strip $S_{\text{u}}\equiv[0,t_{\text{f}}]\times(-\infty,+\infty)$ in the $(t,\kappa)$ plane that go from the point $(0,\kappa_{\text{i}})$ to $(t_{\text{f}},\kappa_{\text{f}})$.

Our starting point is the system of equations (24). We add subscripts u to all the variables to mark that we are studying the unbounded case. Both $\psi_{0,\text{u}}$ and $\psi_{1,\text{u}}$ are constants of motion and thus $y_{1,\text{u}}$ has a linear shape. The boundary conditions for $y_{1}$, as given by (11b), entail that the constant slope equals $(y_{1,\text{f}}-y_{1,\text{i}})/t_{\text{f}}$, i.e.

[TABLE]

and

[TABLE]

In addition,

[TABLE]

Within the theoretical framework of Pontryagin’s maximum principle, the above solution is valid as long as $\tilde{\kappa}$ stemming from (23),

[TABLE]

belongs to the class of admissible controls. It can be easily shown that $\tilde{\kappa}_{\text{u}}(t)\leq 1$ (resp. $\geq 1$) for decompression (resp. compression). Note that, however, $\tilde{\kappa}_{\text{u}}$ may become negative (resp. arbitrarily large) for decompression (resp. compression) as $t_{\text{f}}$ is reduced.

As already stated at the beginning of this section, the class of admissible controls for the unbounded case contains all piecewise functions in the closed interval $[0,t_{\text{f}}]$ that verify the boundary conditions (11a). Therefore, the obtained expression $\tilde{\kappa}_{\text{u}}(t)$ gives the optimal control $\kappa^{*}(t)$ in the open interval $(0,t_{\text{f}})$ but not at the initial and final times. Therein, $\kappa$ is restricted to only one value, $\kappa_{\text{i}}$ for $t=0$ and $\kappa_{\text{f}}$ for $t=t_{\text{f}}$, so it is straightforward that the respective maximums of $\Pi$ are attained at $\kappa^{*}(0)=\kappa_{\text{i}}$ and $\kappa^{*}(t_{\text{f}})=\kappa_{\text{f}}$ 222This can also be understood as having a time dependent control set ${U}$, ${U}(0)=\kappa_{\text{i}}$, ${U}(t)=[0,\kappa_{\max}]$ for $t\in(0,t_{\text{f}})$ and ${U}(t_{\text{f}})=\kappa_{\text{f}}$.. However, this poses no problem because the controls have been assumed to be piecewise continuous in our theory. Therefore, the final result for the optimal control in the unbounded case is

[TABLE]

The optimal profiles for the variables are $y_{0}^{*}(t)=y_{0,\text{u}}(t)$ and $y_{1}^{*}(t)=y_{1,\text{u}}(t)$, with $[0,t_{\text{f}}]$. Neither of them is affected by the finite jumps in $\kappa_{\text{u}}^{*}(t)$, since they are continuous functions of time. Then, we have that

[TABLE]

The above results for the optimal standard deviation and the minimum irreversible work have already been obtained Schmiedl and Seifert (2008); Aurell et al. (2011).

We would like to emphasise the important role played by the boundary conditions to write the relevant variational problem for the physical situation at hand. In the context of ESE processes, one wants to connect the equilibrium states corresponding to $\kappa_{\text{i}}$ and $\kappa_{\text{f}}$ in a finite time $t_{\text{f}}$ and, therefore, the right boundary conditions are those given by (11). Indeed, this is an important issue that affects the result of the variational problem. For example, the boundary conditions considered in Ref. Schmiedl and Seifert (2007) do not connect equilibrium states because the system is not equilibrated at the final time, $\dot{y}(t_{\text{f}})\neq 0$. In fact, this shortcoming was corrected in Ref. Schmiedl and Seifert (2008).

As already stated in the introduction, discontinuities of the optimal stiffness at the initial and final times often appear in stochastic thermodynamics Band et al. (1982); Schmiedl and Seifert (2007, 2008); Aurell et al. (2011). They are usually rationalised in a mathematical way Band et al. (1982), referring to the so-called Miele problem in which the “Lagrangian” is linear in the highest derivative Tolle (2012). We put forward an alternative, physically appealing, picture to understand the emergence of these discontinuities in Appendix A.

IV Bounded stiffness

In experiments, the stiffness of the harmonic trap cannot have an arbitrary value. As stated in the introduction, the type of device employed to design the harmonic potential (AFM, LOT,…) constrains the stiffness values to a certain interval

[TABLE]

For the sake of concreteness and simplicity, we have taken the minimum stiffness as [math] throughout this work. A more general situation with a non-zero $\kappa_{\text{min}}$ can be addressed along similar lines as here. However, note that the most important restriction from a physical point of view is the positiveness of $\kappa$ which, in addition, leads to simpler calculations.

We now turn our attention to the problem of minimising the irreversible work with the non-holonomic constraint (31). In this case, the class of admissible control functions $\mathcal{K}$ comprises all the piecewise continuous functions lying inside the rectangle $S_{\text{b}}\equiv[0,t_{\text{f}}]\times[0,\kappa_{\max}]$ in the $(t,\kappa)$ plane that go from the point $(0,\kappa_{\text{i}})$ to $(t_{\text{f}},\kappa_{\text{f}})$. Evidently, both $\kappa_{\text{i}}$ and $\kappa_{\text{f}}$ must lie in the interval $[0,\kappa_{\max}]$. The maximum value of the stiffness $\kappa_{\max}$ leads to a minimum equilibrium value for the standard deviation, namely

[TABLE]

Pontryagin’s maximum principle is specially adequate to analyse problems with the kind of non-holonomic constraint in (31). The condition $(\partial\Pi/\partial\kappa)_{\tilde{\kappa}}=0$ gives results that are identical to the unbounded case as long as the protocol $\kappa^{*}_{\text{u}}(t)$ lies inside the rectangle $[0,t_{\text{f}}]\times[0,\kappa_{\max}]$. When the optimal protocol for the unbounded case $\kappa^{*}_{\text{u}}(t)$ crosses the boundary of this rectangle at a certain time $t_{0}<t_{\text{f}}$, it can no longer be the solution of the minimisation problem.

Taking into account (23), we have three different regions, $A$, $B$ and $C$, for the optimal stiffness in the bounded case $\kappa_{\text{b}}^{*}$:

[TABLE]

Along the same lines followed in the unbounded case, it is readily shown that $y_{1,\text{b}}^{*}$ is linear in $t$ in region B with slope $-\psi_{1}/(2\psi_{0})$, as predicted by (24). We denote this behaviour by $y_{1,\text{lin}}(t)$. In appendix B, we show that if the system enters region A or region C, it remains there. In other words, once the optimal solution in region B “touches” the boundary at a certain time $t^{J}$, i.e. $\kappa_{\text{b}}^{*}(t^{J})$ equals either [math] or $\kappa_{\max}$, it moves over the boundary from then on, $\kappa_{\text{b}}^{*}(t)=\kappa_{\text{b}}^{*}(t^{J})$ for all $t>t^{J}$. In regions A and C, $y_{1}(t)$ is given by the solutions of (9a) corresponding to constant $\kappa=0$ and $\kappa=\kappa_{\max}$, which we denote by $y_{1}(t)|_{\kappa=0}$ and $y_{1}(t)|_{\kappa=\kappa_{\max}}$, respectively.

On physical grounds, we have three different cases depending on the values of the parameters: namely $\{t_{\text{f}},\kappa_{\text{f}},\kappa_{\max}\}$.

The initial and final states cannot be linked in the given time $t_{\text{f}}$, which is too short. This is due to the impossibility of compressing (resp. decompressing) the system faster than with a STEP protocol with $\kappa(t)=\kappa_{\max}$ (resp. $\kappa(t)=0$). 2. 2.

The time interval $t_{\text{f}}$ is such that the connection is possible but not with the linear solution for the unbounded case $y_{1,\text{u}}^{*}(t)$, since the associated $\kappa^{*}_{\text{u}}(t)\notin[0,\kappa_{\max}]$ for a certain range of times inside $[0,t_{\text{f}}]$. In that case, we show below that the optimal protocol is built as a linear evolution of $y_{1}$ that matches continuously and smoothly (continuous first derivative) the solution of (9a) with $\kappa=\kappa_{\max}$ (resp. $\kappa=0$) in a compression (resp. decompression) process. 3. 3.

The given time $t_{\text{f}}$ is long enough to make the connection possible with the unbounded solution $y^{*}_{1,\text{u}}(t)$ because $\kappa_{\text{u}}^{*}(t)\in[0,\kappa_{\max}]$ for all times. In this case, the bounds do not affect the minimisation problem.

IV.1 Decompression

Let us look into the decompression case, in which $0\leq\kappa_{\text{f}}<1$. To begin with, we would like to discern when $\kappa_{\text{u}}^{*}(t)$ becomes negative. Looking at (28), it is readily seen that the first term on its rhs becomes smaller than the second one for large enough $y_{1}$, and $y_{1,\text{u}}(t)$ increases linearly in time. Therefore, the value of the final time $t_{\text{f}}$ below which the unbounded solution ceases to be valid is determined by the condition $\tilde{\kappa}_{\text{u}}(t_{\text{f}})=0$, i.e. $y_{1,\text{f}}(y_{1,\text{f}}-1)=t_{\text{f}}$. Taking into account (11b), this is equivalent to $t_{\text{f}}=(1-\sqrt{\kappa_{\text{f}}})/\kappa_{\text{f}}$.

The condition $\kappa\geq 0$ implies that there are states that are impossible to connect. The fastest decompression (shortest possible $t_{\text{f}}$) corresponds to a STEP process, in which the stiffness is instantaneously changed to $\kappa=0$ at $t=0^{+}$. In that case, we have that $y_{1}(t)=\sqrt{1+2t}$ and thus $y_{1,\text{f}}=\sqrt{1+2t_{\text{f}}}$. Recalling once more (11b), we conclude that the fastest decompression occurs for $\kappa_{\text{f}}(1+2t_{\text{f}})=1$ or $t_{\text{f}}=(1-\kappa_{\text{f}})/(2\kappa_{\text{f}})$.

Therefore, cases 1, 2 and 3 above correspond here to:

Impossible to connect.

[TABLE] 2. 2.

Matched solution, i.e. a first linear branch $y_{1,\text{lin}}(t)$ and a second branch moving over the line $\kappa=0$ of $S_{\text{b}}$, $y_{1}(t)|_{\kappa=0}$.

[TABLE] 3. 3.

Linear profile for the unbounded case $y_{1,\text{u}}(t)$.

[TABLE]

In case 1, there is no solution and we already know the solution of case 3. Then, we move on to solve case 2, for which the solution comprises two branches. First, a branch corresponding to region B in (33), i.e. a linear profile $y_{1,\text{lin}}$ that verifies only the boundary condition at $t=0$ and thus has one free parameter. This solution is valid in some subinterval $[0,t_{d}^{J}]$, the free parameter can be considered to be its constant slope $m_{d}=-\psi_{1}/(2\psi_{0})$, i.e.

[TABLE]

Second, a branch corresponding to region A in (33), i.e. obtained by putting $\kappa=0$ in (9a), $y_{1}(t)|_{\kappa=0}$. This branch verifies the boundary condition at $t=t_{\text{f}}$ and is valid in the complementary subinterval $[t_{d}^{J},t_{\text{f}}]$. Its specific form is given by

[TABLE]

Note that this branch does not contain any free parameter. The two branches are matched at the joining time $t_{d}^{J}$ by imposing the continuity of both $y_{1}(t)$ and $\dot{y}_{1}(t)$, i.e.

[TABLE]

Note that this is consistent, any solution $y_{1}(t)$ of (9a) must be continuous for piecewise continuous $\kappa(t)$. Moreover, since $\kappa(t)$ is continuous for the matched solution at $t=t_{d}^{J}$, $\kappa({t_{d}^{J}}^{-})=\kappa({t_{d}^{J}}^{+})=0$, $\dot{y}_{1}(t)$ must also be continuous there. We show in appendix B that this simplest approach is the correct one for our problem.

The continuity equations (39) give rise to the conditions

[TABLE]

which can be explicitly solved for $m_{d}$ and $t_{d}^{J}$, with the result

[TABLE]

Note that the matching time $t_{d}^{J}$ is an increasing function of $t_{\text{f}}$, vanishing in the limit as $t_{\text{f}}\to t_{d}^{\min}$ and approaching $t_{\text{f}}$ in the limit as $t_{\text{f}}\to t_{d}^{\text{u}}$. In fact, this solution only makes sense in case 2: in case 1, the argument of the square root is negative whereas in case 3 we have that $t_{d}^{J}>t_{\text{f}}$. Recall that $y_{1,\text{f}}$ is given as a function of $\kappa_{\text{f}}$ by (11b).

Then, the optimal protocol for the stiffness is

[TABLE]

The finite jumps of the stiffness at the initial and final times have the same reason as in the unbounded case and thus we will not repeat the discussion here. The initial jump in the stiffness decreases it to a positive value, $\kappa_{d}^{*}(t=0^{+})=1-m_{d}$, and $0<1-m_{d}<1$. In addition, note that $\kappa_{d}^{*}(t)$ is continuous at $t=t_{d}^{J}$, since the condition $m_{d}=1/y_{1,\text{lin}}(t_{d}^{J})$ holds as a consequence of the continuity of the derivative of $\dot{y}_{1}(t)$ at $t=t_{d}^{J}$, as expressed by (40b). Consistently with our discussion below (41), the expression in (42) only makes sense for $t_{d}^{\min}\leq t_{\text{f}}\leq t_{d}^{\text{u}}$.

Optimal protocols for the decompression case are plotted in Fig. 1. We have chosen $\kappa_{\text{f}}=0.5$ and several values of the connection time $t_{\text{f}}$. The unbounded solution $\kappa_{\text{u}}^{*}(t)$ (dashed lines) only works for the longest time $t_{\text{f}}=t_{d}^{\text{u}}+0.25(t_{d}^{\text{u}}-t_{d}^{\min})$, when it remains positive over the whole time interval. For the remainder of shorter connecting times, $\kappa_{\text{u}}^{*}(t)$ becomes negative as observed in the figure and the optimal protocol equals $\kappa_{d}^{*}(t)$, as given by (42) (thick solid lines). There is no solid line for the longest time, since (42) is well-defined only for $t_{d}^{\min}\leq t_{\text{f}}\leq t_{d}^{\text{u}}$.

IV.2 Compression

When the colloidal particle is compressed, $\kappa_{\max}\geq\kappa_{\text{f}}>1$, the unbounded $\kappa_{u}^{*}(t)$ may become greater than $\kappa_{\max}$. When this is the case, the solution to the minimisation problem is built in a manner completely analogous to the decompression case, but the second branch is obtained by substituting $\kappa_{\max}$ into (9a), i.e. $y_{1}(t)|_{\kappa=\kappa_{\max}}$. Again, the two branches are smoothly joined at a certain $t=t_{c}^{J}$, i.e. with $y_{1}$ and $\dot{y}_{1}$ being continuous.

Since the scenario is analogous to that for decompression, we do not repeat the complete analysis here.

Impossible to connect:

[TABLE] 2. 2.

Matched solution, i.e. at first a linear branch $y_{1,\text{lin}}(t)$ and afterwards moving over the line $\kappa=\kappa_{\max}$ of $S_{\text{b}}$, $y_{1}(t)|_{\kappa=\kappa_{\max}}$,

[TABLE] 3. 3.

Linear profile for the unbounded case $y_{1,\text{u}}(t)$.

[TABLE]

Again, we consider case 2, for which the solution comprises two branches. First, the linear branch $y_{1,\text{lin}}$ valid in $[0,t_{c}^{J}]$, which has the slope $m_{c}=-\psi_{1}/(2\psi_{0})$,

[TABLE]

Second, the branch obtained by substituting $\kappa=\kappa_{\max}$ in (9a), $y_{1}(t)|_{\kappa=\kappa_{\max}}$, which verifies the boundary condition at $t=t_{\text{f}}$ and is valid in $[t_{c}^{J},t_{\text{f}}]$,

[TABLE]

At the joining time $t_{c}^{J}$, $y_{1}(t)$ and $\dot{y}_{1}(t)$ are continuous, which yields

[TABLE]

At variance with the decompression case, this system cannot be explicitly solved for $m_{c}$ and $t_{c}^{J}$ but we can obtain their values for any set of the parameters $\{t_{\text{f}},\kappa_{\text{f}},\kappa_{\max}\}$ numerically. Once more, $y_{1,\text{f}}$ is given by (11b) as a function of $\kappa_{\text{f}}$. Note that $m_{c}<0$ because the standard deviation decreases in time for compression.

Finally, we obtain the optimal protocol for the stiffness in the compression process,

[TABLE]

Again, the initial jump in the stiffness goes in the “right” direction, it increases to $\kappa_{c}^{*}(t=0^{+})=1-m_{c}>1$ because $m_{c}$ as given by (48b) is negative.

Figure 2 is similar to Fig. 1 but for compression. We have chosen the parameter values $\kappa_{\max}=5$ and $\kappa_{\text{f}}=2$. The different curves correspond to different connection times $t_{\text{f}}$. Similarly to the decompression case, the optimal protocol $\kappa_{c}^{*}(t)\neq\kappa_{\text{u}}^{*}(t)$ except for the longest time, since for the remainder of them $\kappa_{\text{u}}^{*}(t)$ violates the inequality $\kappa\leq\kappa_{\max}$. Also, the matching time $t_{c}^{J}$ increases with $t_{\text{f}}$, $t_{c}^{J}\to 0$ as $t_{\text{f}}\to t_{c}^{\min}$ and $t_{c}^{J}\to t_{\text{f}}$ as $t_{\text{f}}\to t_{c}^{\text{u}}$. Similarly to the decompression case, (48) and (49) only make sense in case 2.

V Phase diagram and average work

V.1 Inaccessible and accessible states

Depending on the values of the parameters $\{t_{\text{f}},\kappa_{\text{f}},\kappa_{\max}\}$, we have three different “phases” when the stiffness is bounded, which correspond to each of the cases enumerated in the previous section. For each value of the maximum stiffness $\kappa_{\max}$, there are target points $(\kappa_{\text{f}},t_{\text{f}})$ that

are inaccessible; there is no control $\kappa(t)$ capable of linking the initial and final states, 2. 2.

can be reached by means of a matched solution; the optimal control moves partially over the boundary of the rectangle $S_{\text{b}}\equiv[0,t_{\text{f}}]\times[0,\kappa_{\max}]$, and 3. 3.

can be reached with the optimal control for the unbounded case; the standard deviation has the simple linear form $y_{1,\text{u}}(t)$ in (26).

In order to look into the different phases, it is worth going to the natural time scale for relaxation at the final stiffness $\kappa_{\text{f}}$, i.e we define

[TABLE]

In this time scale, the equilibration time is the same for all $\kappa_{\text{f}}$, $\tau_{\text{eq}}\simeq 3$, see (13). Therefore, the value of the connection time in the $\tau$ scale, $\tau_{\text{f}}=\kappa_{\text{f}}t_{\text{f}}$, directly gives the acceleration of the ESE process with respect to the STEP one. The times separating the different regions (inaccessible, matched, unbounded) are readily obtained from (34) and (35) in decompression

[TABLE]

and (43) and (44) in compression,

[TABLE]

In Fig. 3, we plot the different regions in the plane $(\tau_{\text{f}},\kappa_{\text{f}})$ for the specific case $\kappa_{\max}=50$. We have shaded regions in (i) grey, (ii) red and (iii) green. The dashed red lines separate regions (i) and (ii), i.e. they are given by $\tau_{\text{f}}=\tau_{p}^{\min}$, where $p=c$ or $d$ depending on the type of process, compression or decompression. The optimal protocol over these lines is an initial abrupt change from $\kappa(0)=\kappa_{\text{i}}=1$ to $\kappa=0$ in the decompression process (to $\kappa=\kappa_{\max}$ for compression) and another sudden jump from this value to the target stiffness $\kappa_{\text{f}}$ at the final time. The solid green lines separate regions (ii) and (iii), i.e. they are given by $\tau_{\text{f}}=\tau_{p}^{\text{u}}$, again $p=c$ or $d$ depending on the type of process. Over these lines, the unbounded solution becomes valid throughout the whole time interval, reaching the border ($\kappa=0$ for decompression, $\kappa=\kappa_{\max}$ for compression) at $t=t_{\text{\text{f}}}^{-}$.

The decompression case deserves further commenting. Both the minimum connection time $\tau_{d}^{\min}$ and the time above which the unbounded solution is valid $\tau_{d}^{\text{u}}$ are bounded, specifically

[TABLE]

As it is clearly seen in Fig. 3, this means that $\tau_{d}^{(2)}=1$ is a critical time in decompression: above it, the initial equilibrium state can always be connected to another equilibrium state corresponding to an arbitrary value of the stiffness with the protocol valid for the unbounded case. Moreover, $\tau_{d}^{(1)}=1/2$ is a second critical time in decompression: for $\tau_{d}^{(1)}<\tau_{\text{f}}<\tau_{d}^{(2)}$, all the final equilibrium states are accessible but the unbounded solution is only valid for weak enough decompression—meaning $\kappa_{\text{f}}$ smaller but not too far from unity, whereas for $\tau_{\text{f}}<\tau_{d}^{(1)}$ there appear inaccessible states. This is to be contrasted with the compression case. In this latter case, the three possible phases, inaccessible, matched and unbounded, are possible for all the connecting times $\tau_{\text{f}}$.

Note that the existence of an upper (resp. lower) bound on $\kappa$ does not affect the decompression (resp. compression) case. This stems from the monotonicity of the optimal protocols for the stiffness, as explicitly proven in Appendix B.

V.2 Properties of the mean work

At this point, it is worth looking into the optimal average work and elucidate how the problem changes upon constraining the stiffness $\kappa$. The optimal value for the irreversible work $W^{*}$ can be computed in regions (ii) and (iii), when the connection between the initial and final states is possible. In region (iii), the bound on the stiffness plays no role for calculating $W_{\text{irr}}^{*}$, $W_{\text{irr}}^{*}=W_{\text{irr},\text{u}}^{*}$ as given by (30). In region (ii), we have to use the matched solutions in Secs. IV.1 and IV.2 to derive the minimum work. We employ again $p=c$ or $d$ to label the kind of process. The integral in (12) is split into two parts: the first one from [math] to $t_{p}^{J}$, where $y_{1}(t)$ is linear in time with slope $m_{p}$, and the second one from $t_{p}^{J}$ to $t_{\text{f}}$, where $y_{1}(t)$ is given by the boundary solution $y_{1}(t)|_{\kappa=\kappa_{p}}$; $\kappa_{p}$ stands for the relevant boundary value of $\kappa$, $\kappa_{d}=0$ and $\kappa_{c}=\kappa_{\max}$. Then,

[TABLE]

Integrating over $y_{1}$ instead of $t$ in the second term of the rhs and making use of (9b) and the continuity of $y_{1}$ at the matching time, $y_{1,p}(t_{p}^{J-})=y_{1,p}(t_{p}^{J+})=y_{1,p}^{J}$, one gets

[TABLE]

Let us particularise (55) for decompression and compression. First, in the decompression case we have that

[TABLE]

in which $m_{d}$ and $t_{d}^{J}$ are given by (41) in terms of $(t_{\text{f}},\kappa_{\text{f}})$, and $y_{1,d}^{J}$ is the value of $y_{1}$ at the joining time $t_{d}^{J}$ as defined in (40a). Second, for compression we obtain

[TABLE]

where $m_{c}$ and $t_{c}^{J}$ are the solutions of the system of equations (48b), and $y_{1,c}^{J}$ is given by (48a).

In what follows, we plot with dashed lines the optimal work coming from the unbounded expression, as given by (30). Solid lines are used for the optimal work when the bound $0\leq\kappa\leq\kappa_{\max}$ is relevant, (56) for decompression and (57) for compression. In addition, we have shaded the different regions with the same colour code employed in the phase diagram. The solid lines are always above the dashed ones, because the minimum with no constraints is logically lower than the constrained one.

First, we investigate the optimal work as a function of the final time $\tau_{\text{f}}$, for different values of $\kappa_{\text{f}}$. Specifically, we consider a compression protocol with $\kappa_{\text{f}}=2$ and a decompression protocol with $\kappa_{\text{f}}=0.5$ in Fig. 4. The stiffness is bounded in the interval $0\leq\kappa\leq\kappa_{\max}=5$. The difference between the constrained and unconstrained optimal values of the work becomes more important as the connection time $\tau_{\text{f}}$ decreases, as discussed below.

Let us investigate the decompression case in more detail. We focus on the difference between the actual optimal work $W_{\text{irr},d}^{*}$ and its value for the unconstrained case $W_{\text{irr},\text{u}}^{*}$ for the minimum connection time $\tau_{\text{f}}\to\tau_{d}^{\min}$ (or $t_{\text{f}}\to t_{d}^{\min}$), which is given by (51). At this point, this difference reaches its maximum value. Therefore, the first term on the rhs of (56) for the optimal work $W_{\text{irr},d}^{*}$ does not contribute thereto, because $t_{d}^{J}\to 0$, and we have

[TABLE]

In Fig. 5 we plot the relative difference $(W_{\text{irr},d}^{*}-W_{\text{irr},\text{u}}^{*})/W_{\text{irr},\text{u}}^{*}$ as a function of $\kappa_{\text{f}}$. It remains small for $\kappa_{\text{f}}\gtrsim 0.3$, for instance for $\kappa_{\text{f}}=0.5$ it is below $1\%$. As $\kappa_{\text{f}}$ decreases it starts to grow; in fact, $W_{\text{irr},d}^{*}$ diverges as $\kappa_{\text{f}}\to 0$. For example, for $\kappa_{\text{f}}=0.1$ the relative difference is around $10\%$, for $\kappa_{\text{f}}=0.01$ it has increased to $40\%$ and for $\kappa_{\text{f}}=10^{-3}$ it exceeds $80\%$.

Second, we study the optimal work as a function of the compression ratio $\kappa_{\text{f}}$ for a fixed value of the connection time $\tau_{\text{f}}$. Similarly to the situation found when $\tau_{\text{f}}$ was varied for fixed $\kappa_{\text{f}}$, we have again inaccessible, matched and unbounded regions. Specially interesting is the decompression case, in principle the minimum value of the stiffness $\kappa_{d}^{\min}$ for having connected states and the value $\kappa_{d}^{\text{u}}$ above which the unbounded solution works should be obtained by using (51). Notwithstanding, the situation is a little more complex. Specifically we have that

[TABLE]

The piecewise definitions of $\kappa_{d}^{\min}$ and $\kappa_{d}^{\text{u}}$ are readily rationalised by looking at Fig. 3: obtaining $\kappa_{d}^{\min}$ (51) only makes sense as long as $\tau_{\text{f}}\leq\tau_{d}^{(1)}=1/2$, above it $\kappa_{d}^{\min}=0$ because all the states with $\tau_{\text{f}}>\tau_{d}^{(1)}$ are accessible. A similar reasoning applies to $\kappa_{d}^{\text{u}}$: for $\tau_{\text{f}}>1$, all the states can be connected with the unbounded solution. The most interesting region in the ESE context is $\tau_{\text{f}}<\tau_{d}^{(1)}$, which corresponds to the higher acceleration of the equilibration process, $\tau_{\text{f}}/\tau_{\text{eq}}<1/6$.

Figure 6 corresponds to the specific case $\kappa_{\max}=5$ and $\tau_{\text{f}}=1/3$. Therefore, the connection time is roughly one tenth of the equilibration time, $\tau_{\text{f}}/\tau_{\text{eq}}=1/9$, and we plot compression ($\kappa_{\text{f}}>1$) and decompression ($\kappa_{\text{f}}<1$) processes in the same graph. For these values of the parameters, the main effect of the bounds is the reduction of the effectively accessible region for $\kappa$, which is much smaller than the whole interval $[0,\kappa_{\max}]$. The matched solutions are needed in two layers close to the borders of the accessible region, but the differences between the bounded optimal work and the unbounded value are quite moderate.

We consider a value of the connecting time close to the critical value $\tau_{d}^{(1)}$ in Fig. 7, specifically $\tau_{\text{f}}=0.498$. The inaccessible region becomes very small, since $\kappa_{d}^{\min}=0.004$ but the bounded irreversible work $W_{\text{irr},d}^{*}$ is about $60\%$ higher than the unbounded irreversible value $W^{*}_{\text{irr},\text{u}}$ at $\kappa_{\text{f}}=\kappa_{d}^{\min}$. In fact, as $\tau_{\text{f}}\to\tau_{d}^{(1)}$ we have that $\kappa_{d}^{\min}\to 0$ and the corresponding $W_{\text{irr},d}^{*}$ diverges logarithmically whereas the bounded value remains finite, $W_{\text{irr},\text{u}}^{*}\to 2$, as expressed by (58). We further illustrate this fact by plotting both $W_{\text{irr},d}^{*}$ and $W_{\text{irr},\text{u}}^{*}$ at $\kappa_{\text{f}}=\kappa_{d}^{\min}$ as a function of $\tau_{\text{f}}$ in Fig. 8. It is observed that, consistently with the discussion above and the picture shown in Fig. 5, the difference between the two are largest for $\tau_{\text{f}}\to\tau_{d}^{(1)}$. In principle, it may seem strange that $W_{\text{irr},d}^{*}$ and $W^{*}_{\text{irr},\text{u}}$ tend to coincide in the limit as $\tau_{\text{f}}\to 0$. Looking once more at Fig. 3, it is seen that the inaccessible region grows as $\tau_{\text{f}}$ is decreased and fills the whole decompression region as $\tau_{\text{f}}\to 0$, i.e. $\kappa_{d}^{\min}\to 1$ in this limit. Therefore, a very large acceleration of the process is only possible in the linear response regime $1-\kappa_{\text{f}}\ll 1$, for which both works are infinitesimally small. In fact, their relative difference can also be shown to be very small.

VI Conclusions

In experiments with confined colloids, a natural constraint on the trap stiffness is that expressed by (5). This non-holonomic constraint makes it impossible to solve the minimisation problem of the work by employing the usual approach involving the Euler-Lagrange equations. Instead, it is necessary to address the problem by employing the tools of control theory, specifically Pontryagin’s maximum principle. Interestingly, a similar approach based on control theory has been recently applied to address the minimisation of entropy production in the trapped colloidal particle problem Muratore-Ginanneschi and Schwieger (2017), but with “bounded accelerations”. The relevance of these bounds, which were originally introduced to regularise the jumps of the stiffness at the initial and final times Aurell et al. (2012), for experiments is not obvious.

The bounds on the stiffness strongly modify the problem of minimising the work performed on the colloidal particle. The solution for unbounded stiffness, in which the standard deviation $y_{1}$ connects linearly the initial and final states, is no longer valid in general: the associated optimal stiffness $\kappa_{\text{u}}^{*}(t)$ violates the inequality (5) for short enough connecting times $t_{\text{f}}$. First, there appear minimum times for connecting the initial and final states, since it is impossible to compress (resp. decompress) the system with any control $\kappa(t)$ faster than with the one corresponding to $\kappa(t)=\kappa_{\max}$ (resp. $\kappa(t)=0$) for all times.

Second, and most importantly, for times longer than the minimum time but smaller than a certain time, there exists an optimal control $\kappa^{*}_{\text{b}}(t)$ but it is different from $\kappa^{*}_{\text{u}}(t)$. This is the significant time window for ESE protocols, since we need the connection to be possible but with the shortest possible time. The associated time evolution for the standard deviation comprises two branches. First, a linear branch, $y_{1,\text{lin}}(t)$ where $y$ denotes the position standard deviation, in the first part of the time interval while $\kappa(t)$ has not reached the bounds yet. Second, a branch corresponding to the solution for the appropriate boundary value of $\kappa$ ($\kappa_{\max}$ in compression, [math] in decompression) in the second part of the time interval. The two functions match smoothly, with $y_{1}(t)$ and $\dot{y}_{1}(t)$ being continuous, at the joining time.

Rather dramatic changes are observed in the decompression case, when the bound $\kappa\geq 0$ comes into play. This is not a mathematical bound but a physical one: with a harmonic trap, it is experimentally difficult to engineer a repulsive potential, and thus the stiffness has to remain non-negative. Most importantly, there appear two critical values for the connection time: for $\tau_{\text{f}}<\tau_{d}^{(2)}$ the bounded optimal work $W_{\text{irr},d}^{*}$ deviates from that of the unbounded problem $W_{\text{irr},\text{u}}^{*}$, and at $\tau_{\text{f}}=\tau_{d}^{(1)}$ we have that $W_{\text{irr},d}^{*}$ diverges.

In the last decade, stochastic heat engines have been designed by trapping a Brownian particle in a harmonic potential, the stiffness of which can be externally controlled Schmiedl and Seifert (2008); Blickle and Bechinger (2012); Martínez et al. (2016b); Ciliberto (2017), i.e. the physical system investigated here. The cycles considered in these Brownian heat engines typically comprise four branches, with two of them being isothermal compression and isothermal decompression processes. The work over these isothermal processes must be minimised to maximise the power delivered by the engine—the work performed by the system is minus the work performed on the system, which is the one considered throughout this paper.

The changes in the optimal work derived here for isothermal compression/decompression processes, which are entailed by the bounds in the stiffness, impinge on the optimal power of the Brownian heat engines. Specifically, the optimal power is lowered as compared with the value obtained for unbounded stiffness. In this respect, analysing in detail the impact of the bounds on the power of heat engines constitutes an interesting prospect for future research. Another relevant venue lies in optimising mixed quantities, such as a combination of the mean dissipated work and its standard deviation, which may exhibit phase transitions in protocol space Solon and Horowitz (2018). Also, in the realm of microfluidics Sajeesh and Sen (2014), it seems interesting to explore the extension of the ideas presented here to the design of optimal devices for separating and sorting particles in a desired time.

Appendix A “Surgery” method for the unbounded case

Here we deal with the optimisation of the work in the unbounded case from an alternative point of view. In absence of the non-holonomic constraint $0\leq\kappa\leq\kappa_{\max}$, one may hope to address the optimisation problem by employing the classical variational approach leading to the Euler-Lagrange equations. Below we show the difficulties that arise and how to cope with them by a physically appealing “surgery” procedure 333This “surgery” can be thought of as a particular case of the procedure explained in section 2.2 of Ref. Tolle (2012) for the Miele problem..

We start by writing the irreversible work as

[TABLE]

where the boundary conditions for $y_{1}$ are given by (11b). Therefore, this seems to be a “trivial” problem: the Euler-Lagrange equation for the optimal profile $y_{1}^{*}$ is simply $\ddot{y}^{*}_{1}=0$ and its solution is exactly (26). The issue arises now, because the optimal stiffness $\kappa(t)$ obtained from (9),

[TABLE]

does not verify the boundary conditions (11a). Note that these boundary conditions for $\kappa$ are equivalent to $\dot{y}_{1}(0)=\dot{y}_{1}(t_{\text{f}})=0$, i.e. they ensure that the system is properly equilibrated at both the initial and final states 444The fact that the linear solution (26) verifies the boundary conditions for $y_{1}$ but not those for $\dot{y}_{1}$ is not surprising mathematically: the Euler-Lagrange equation is a second order differential equation and minimises the irreversible work for given values of $y_{1}$ at the boundaries. Then, there seems to be no room for “extra” boundary conditions..

From a physical point of view, there should be an optimal procedure—in the sense that the irreversible work attains a minimum over it—to connect the initial and final equilibrium states in a finite time. Therefore, there should be a time evolution for $y_{1}$ that minimises the irreversible work and verifies both the boundary conditions for $y_{1}$ and $\dot{y}_{1}$, i.e. a solution of the overdetermined problem

[TABLE]

Below we show that this is indeed the case by explicitly building a solution. With this constructive procedure, what we basically reveal is that the extra boundary conditions for $\dot{y}_{1}$ do not change the solution in $(0,t_{\text{f}})$: it suffices to bring to bear the boundary conditions for $y_{1}$ and introduce suitable jumps in $\dot{y}_{1}$ at the boundaries. Note that we have omitted the asterisk in the solution of the variational problem, i.e. we have written $y_{1}$ instead of $y_{1}^{*}$ in order not to clutter (62).

To keep expressions simpler, first we introduce suitable rescalings for both $y_{1}$ and $t$,

[TABLE]

such that

[TABLE]

where the prime indicates derivative with respect to $s$, and the overdetermined problem in (62) is

[TABLE]

We build the family of functions in the half interval $s\in[1/2,1]$. We split the interval $[1/2,1]$ into two parts, $[1/2,1-\epsilon]$ and $[1-\epsilon,1]$, and write down the following family of piecewise defined functions

[TABLE]

with

[TABLE]

It is easily shown that that the functions $u_{\epsilon}(s)$ (i) satisfy the boundary conditions at the right endpoint $s=1$ in (65) for all $\epsilon$ and (ii) are continuous and have continuous derivative in $[1/2,1]$, including the connection point $s=1-\epsilon$. Nevertheless, in the limit as $\epsilon\to 0^{+}$ the “boundary layer” $[1-\epsilon,1]$ collapses with $u_{\epsilon}(s)$ remaining continuous at the endpoint $s=1$ but $u_{\epsilon}^{\prime}(s)$ becoming discontinuous. Specifically,

[TABLE]

because $\mu_{\epsilon}\to 1$ in the considered limit. Therefore, in the limit as $\epsilon\to 0$ we generate the discontinuity in the derivative of $u$—and thus of $\dot{y}_{1}$ and $\kappa$.

In the other half interval $s\in[0,1/2]$ the function is defined by a “mirroring” process (both left-right and up-down) with respect to the central point $s=1/2$, $u=1/2$, i.e.

[TABLE]

The boundary conditions at $s=0$ are automatically fulfilled as a consequence of the boundary conditions at $s=1$.

It is a matter of simple algebra to show that

[TABLE]

Therefore, the irreversible work for this family of functions approaches the minimum value for the standard problem—with only the values of $u$ fixed at the boundaries—as $\epsilon$ goes to zero. Since the minimum in the overdetermined problem, with extra conditions on the derivative, cannot be smaller than that for the standard problem, we conclude that the solution for the overdetermined problem is given by $\lim_{\epsilon\to 0^{+}}u_{\epsilon}(s)$. In other words, the solution for the standard problem with a sudden finite jump at the boundary.

Figure 9 shows the corresponding stiffness protocols $\kappa_{\epsilon}(t)$, as given by inserting the family $u_{\epsilon}(t)$ into (61), for several values of $\epsilon$. They are compared with the solution $\kappa_{\text{u}}^{*}(t)$ that we calculate in the main text by applying Pontryagin’s principle, which has finite jumps at the boundaries. It is neatly observed how the proposed surgery procedure recovers the solution $\kappa_{\text{u}}^{*}(t)$ in the limit as $\epsilon\to 0^{+}$, including the jump at the boundary—although $\kappa_{\epsilon}(t)$ is continuous and has continuous derivative everywhere for any finite $\epsilon$.

Appendix B Derivation of the solution for the bounded

case

Let us consider the solution for the bounded case in more detail. We focus on the decompression case because the calculations are simpler as a consequence of our choosing $\kappa_{\min}=0$. In the main text, we have built the optimal solution $\kappa^{*}(t)$ by assuming that (i) when $\kappa^{*}$ touches the boundary $\kappa=0$, then $\kappa^{*}$ remains over the boundary for longer times, and (ii) the upper bound $\kappa_{\max}$ plays no role in the decompression problem. In what follows, we show that this is indeed the case.

On the one hand, the solution of the system of equations (24) provides the optimal time evolution inside those time windows such that the corresponding stiffness $\tilde{\kappa}$ calculated from (23) remains non-negative, i.e.

[TABLE]

In those time windows, the optimal stiffness is $\kappa^{*}=\tilde{\kappa}$ and as a consequence $y_{1}$ is a linear function of time. Note that the left hand side of (71) above is simply $\tilde{\kappa}y_{1}^{2}$. On the other hand, inside the time windows for which

[TABLE]

we have that $\kappa^{*}=0$ and $y_{1}(t)$ satisfies the particularisation of (9) to $\kappa=0$.

In principle, there may appear a number of different time windows with several joining times $[0,t_{1}]$, $[t_{1},t_{2}]$, $\ldots$, $[t_{n},t_{\text{f}}]$, with the solution changing from linear to the $\kappa=0$ case (or vice versa) at each of the joining times $t_{k}$. Now we prove that there is only one joining time $t_{1}$ ($t_{d}^{J}$ in the notation of the main text) by establishing that once $\tilde{\kappa}=0$ at a certain time, the optimal control $\kappa^{*}$ remains over the boundary. To do this, we show that

[TABLE]

Since the condition $\frac{\psi_{1}}{2\psi_{0}}y_{1}+1\leq 0$ implies that $\tilde{\kappa}\leq 0$, it suffices to prove that $\frac{\psi_{1}}{2\psi_{0}}y_{1}+1$ is decreasing for $\kappa=0$. By using the evolution equations for that case, it is readily shown that

[TABLE]

Second, we explain why the upper bound plays no role in the decompression process. At stated in the main text, at $t=0^{+}$ the stiffness coming from the proposed solution is positive and lower than $\kappa_{\text{i}}=1$; (42) leads to $\kappa_{d}^{*}(t=0)=1-m_{d}$, and $0<m_{d}<1$. Moreover, in the “linear” time window $[0,t_{d}^{J}]$ the stiffness monotonically decreases, because

[TABLE]

Initially, this derivative is negative because $y_{1}(t=0)=1$. The term in parentheses increases linearly but at the joining time $m_{d}y_{1}(t_{d}^{J})=1$, see (40), so the derivative is still negative. Therefore, it does not change sign in the interval $[0,t_{d}^{J}]$, being always negative therein. Since once it touches the boundary, $\kappa_{d}^{*}$ remains constant, we have

[TABLE]

Thus, it is clear that the upper bound $\kappa_{\max}$ is irrelevant when finding the optimal stiffness protocol for decompression, $\kappa_{d}^{*}(t)<1$ for all times.

Along similar lines, it is shown that the solution given in the main text is the only one for compression: there is also only one connecting time and the lower bound $\kappa_{\min}=0$ is irrelevant in that case. The calculations are a little bit lengthier and are thus not given here.

Acknowledgements.

C.A.P and A.P. acknowledge the support of Universidad de Sevilla’s VI Plan Propio de Investigación through Grant PP2018/494. C.A.P. also acknowledges the support from the FPU Fellowship Programme of the Spanish Ministerio de Educación, Cultura y Deporte through Grant FPU14/00241 and from STARS2018 project through UNIPD. D.G.O. and E.T. acknowledge financial support from the Agence Nationale de la Recherche (research funding Grant No. ANR-18-CE30-0013-03). E.T. also acknowledges funding from the Investissement d’Avenir LabEx PALM program (Grant No. ANR-10-LABX-0039-PALM).