Optimal protocols for minimal work processes in underdamped stochastic thermodynamics

I Introduction

The free energy difference $\Delta F$ between two equilibrium states is an important quantity in isothermal statistical mechanics. Strategies to extract $\Delta F$ from experiments or computer simulations are traditionally based on either thermodynamic integration or thermodynamic perturbation Zwanzig (1954) which use one infinitesimally slow transition or many infinitesimally fast transitions, respectively, between the two equilibrium states. A decade ago, Jarzynski proposed the remarkable relation

which interpolates between these extreme cases using nonequilibrium work values $W$ obtained from trajectories of finite time transitions between the equilibrium states at temperature $T$ (with Boltzmann’s constant $k_{B}=1$ throughout the paper) Jarzynski (1997a, b). This exact relation which holds for any time-dependent driving described by an external control parameter $\lambda(t)$ has been extended to various fluctuation theorems Crooks (1999, 2000); Hummer and Szabo (2001); Hatano and Sasa (2001); Seifert (2005). Although these (necessarily irreversible) finite time transitions occur in nonequilibrium, the equilibrium quantity $\Delta F$ can be inferred from a sufficient number of trajectories either from computer simulations Hendrix and Jarzynski (2001); Hummer (2001); Shirts et al. (2003); Park and Schulten (2004); Maragakis et al. (2006) or real experiments Liphardt et al. (2002); Collin et al. (2005). However, the convergence of the involved exponential average causes problems for far out of equilibirium transitions where the work $W$ is substantially larger than the free energy difference $\Delta F$ Gore et al. (2003). In this regime, the exponential average is dominated by low work values which are very rarely sampled Jarzynski (2006). As a remedy, several path sampling techniques biasing the dynamics for low work have been proposed Ytreberg and Zuckerman (2004); Atilgan and Sun (2004); Lechner et al. (2006). It is, however, still under debate Oberhofer et al. (2005); Ytreberg et al. (2006); Lechner and Dellago (2007) for which systems fast growth techniques are superior to refined “conventional” approaches such as umbrella sampling Torrie and Valleau (1977) or flat histogram methods Wang and Landau (2001). Though valuable for computer simulations, it is hard to imagine how to bias dynamics in real experiments, where, however, apparatus drift may prevent long measurements necessary for thermodynamic integration Collin et al. (2005); Maragakis et al. (2007) and thus render fast growth methods competitive.

Both for thermodynamic integration and “fast growth” methods employing Jarzynski’s equality (or some variant), efficiency gains can be achieved by optimizing the driving scheme $\lambda(t)$. For thermodynamic integration, where the work $W\geq\Delta F$ is taken as an estimator for $\Delta F$, it is obvious that a minimal work gives the best result de Koning (2005). In the case of fast growth methods, the statistics for free energy estimates quite generally also improves with smaller mean work Atilgan and Sun (2004); Jarzynski (2006). Fast growth simulations even allow to combine data from different driving schemes in a straightforward way Minh (2006).

The minimization of the work spent in a finite time process can, however, also be seen in the context of minimal energy dissipation. On a macroscopic scale, the optimization of the work (or power) exerted in a macroscopic cyclic process has been discussed for quite a while Curzon and Ahlborn (1975); Band et al. (1982); Andresen et al. (1984); Hoffmann et al. (2003). On a microscopic scale, fluctuations will also affect optimal cyclic processes Schmiedl and Seifert (2008) which may become relevant for constructing optimal nano machines.

For overdamped Langevin dynamics, the optimal protocol leading to a minimal mean work in a finite time $t$ has been calculated analytically for harmonic potentials Schmiedl and Seifert (2007). Surprisingly, the optimal protocol shows jumps at the beginning and end of the finite time transition. Since most molecular dynamics (MD) simulations of the dynamics of biomolecules are on time-scales where inertia plays an important role (see Sotomayor and Schulten (2007) for a review on steered MD), it is an interesting question how these results transfer to underdamped dynamics. In particular, it is important to know whether the jumps are a result of having neglected inertia.

In this paper, we investigate the role of inertia for two previously introduced paradigmatic processes. In Sect. II, we calculate optimal driving schemes for an underdamped Brownian particle dragged through a viscous fluid by harmonic optical tweezers. In Sect. III, we study an underdamped Brownian particle subject to an optical trap with time-dependent stiffness. In both cases, we compare our findings with the corresponding results in the overdamped limit Schmiedl and Seifert (2007). We find that the optimal protocol still involves jumps. Even more surprisingly, the optimal protocol includes delta peaks at the beginning and end of the process.

II Case Study I: The moving trap

We consider a Brownian particle of mass $m$ dragged through a viscous fluid with friction coefficient $\gamma$ by a harmonic potential

where $k$ is the (constant) trap stiffness. The focus of the optical trap $\lambda(t)$ is changed time-dependently from an initial position $\lambda_{i}=0$ to a final position $\lambda_{f}$. Including inertial effects, the Langevin equation reads

where thermal fluctuations are modeled by Gaussian white noise

The mean work spent in the process of total duration $t_{f}$ is given by

where the average $\langle\dots\rangle$ is over the intitial thermal distribution and over the noise history. In the present case the total (mean) work reduces to

where, for simplicity in the notation, we have defined the mean position of the particle as $u\equiv\langle x\rangle$. This quantity $u(t)$ depends on the whole history of $\lambda(t)$ and thus, the work $W$ is a non-local functional of the protocol $\lambda(t)$. However, in analogy to the overdamped limit Schmiedl and Seifert (2007), we can express the work as a local functional of the mean particle position $u$. By averaging the evolution equation (3) we have

The only term remaining in the integral, $\dot{u}^{2}$, is identical to the one in the overdamped limit, while the boundary terms are different. In complete analogy to the overdamped case, we now proceed in two steps. First, we calculate the optimal shape $u(t)$ minimizing only the integral given initial values $u(0^{+})=0$ and $\dot{u}(0^{+})=A$. Note that despite the initial equilibrium value $\dot{u}(0)=0$, we are free to choose $\dot{u}(0^{+})=A$ since the necessary “kink” in $u(t)$ at $t=0$ does not contribute to the integral. Similarly, at the end of the protocol (at $t=t_{f}$) there can be another jump in the velocity. In a second step, we adjust the constant $A$ to yield the minimal total work. First, from the Euler-Lagrange equation corresponding to the Lagrangian $\dot{u}^{2}$ (and subject to the initial conditions just mentioned), we find

for $0<t<t_{f}$. In contrast to the overdamped case, we cannot determine all the boundary terms at $t_{f}$ from the evolution equation. Thus, $C\equiv\dot{u}(t_{f})$ is another free parameter. With $\ddot{u}(t_{f})=[k(\lambda_{f}-At_{f})-\gamma C]/m$, we get the total work as a function of the yet unknown constants $A$ and $C$

The work is clearly minimal for $C^{*}=0$, where the asterisk will denote optimal from now on. The remaining terms then read

which, surprisingly, is exactly the same expression that was found in the overdamped limit. Minimizing this expression with respect to $A$ leads to

Inserting Eq. (9) into Eq. (7), we find the optimal protocol

for $0<t<t_{f}$ implying symmetrical jumps

at the beginning and at the end of the process.

Superficially, this optimal protocol looks like the expression in the overdamped case Schmiedl and Seifert (2007). There is, however, a subtle difference arising from the presence of inertia terms. The optimal protocol forces the mean velocity to instantly jump at the beginning of the process from its initial equilibrium value $\dot{u}(0)=0$ to $\dot{u}(0^{+})=A^{*}$. At the end of the protocol, the optimal strategy consists in setting back the mean velocity to zero $\dot{u}(t_{f})\equiv C^{*}=0$. Due to the second time derivative in the equation of motion such jumps in the velocity, which require delta functions in the acceleration, imply a delta-type singularity in the protocol. Specifically, in Eq. (7), the jumps in $\dot{u}$ imply a $\delta$-function for $\ddot{u}$ and hence a $\delta$ function in $\lambda(t)$. The optimal protocol [Eq. (14)] thus becomes

as shown in Fig. 1. In the overdamped limit, $m\to 0$, the delta peaks vanish.

II.2 Physical origin of singularities in the optimal protocol

The benefit of having jumps in the optimal protocol can be understood intuitively as follows. From the perspective of minimal dissipation, it is obvious that the particle should be dragged at a constant (mean) velocity from the beginning rather than being accelerated during a finite time. This initial jump in the velocity of the particle can only be achieved by a finite initial difference $\lambda(0)-u(0)$, corresponding to a jump in $\lambda$ at $t=0$. In the present underdamped regime, a velocity jump corresponds to a jump in the (mean) particle momentum which can only be achieved by a delta peak in the force, corresponding to a delta peak in the protocol.

The final jump is harder to grasp intuitively. In fact, it stems from focussing on the minimal work rather than on the minimal (mean) dissipation (or entropy production). If we had searched for the minimal entropy production (as defined in Seifert (2005)), we would have found an optimal protocol without a final jump. In the present minimization, at the final time $t_{f}$, the particle is still in non-equilibrium with respect to the final potential $V(x,\lambda(t_{f}))$. Relaxation to equilibrium leads to further dissipation after time $t_{f}$ which has, however, already been paid for by the total work since at constant $\lambda$ no work is exerted anymore. A smaller final particle position $u(t_{f})$ leads to a longer relaxation time which can decrease the total dissipation of the combined process (nonequilibrium transition and relaxation).

The final delta peak corresponds to setting the final velocity to zero. This decreases the kinetic energy of the particle and thus is beneficial for a small work. It also explains the surprising fact that, according to Eq. (13), we do not have to pay any extra cost for having inertia. During the initial singularity, the exerted work is stored in the (mean) kinetic energy of the particle. This contribution is fully recovered during the final singularity where the kinetic energy of the particle is set back to the equilibrium value.

II.3 Comparison to a linear protocol

Without prior knowledge, one might have expected a continuous linear protocol

Solving the second order differential equation of motion (7) using the linear protocol $\lambda^{\rm lin}(t)$, we find the ratio:

III Case Study II: The Stiffening trap

In the first case study, only the averaged quantity $u=\left\langle x\right\rangle$ appeared in the work and thus the same result could have been obtained from a deterministic damped dynamics. We next examine a second case study where fluctuations are important. We consider a Brownian particle of mass $m$ in an optical trap with time dependent stiffness $\lambda(t)$ which is driven from an initial value $\lambda(0)=\lambda_{i}$ to a final value $\lambda(t_{f})=\lambda_{f}$ in a finite time $t_{f}$. The time dependent potential

leads to the underdamped Langevin equations

where $p$ is the momentum of the particle and the noise $\eta(t)$ has the same properties introduced in the first case study. Again our main goal is to find the protocol for which the corresponding total (mean) work

is minimal. Note that the mean squared position

of the particle is non-trivially coupled to the mean squared momentum

Their time evolution is governed by the set of coupled differential equations

Unlike both the moving trap (with and without inertia) and the stiffening trap in the overdamped limit, the present case is much more involved since one cannot eliminate the protocol and write the work as a function of one variable only. We thus express the work as a time-local functional of $w(t)$ and $z(t)$. Solving Eqs. (27) and (29) for $\lambda$ yields

which, inserted in Eq. (23) and after partial integration, leads to

We proceed in two steps analogously to the moving trap. We first minimize the integral in Eq. (31) for given initial conditions and then optimize with respect to remaining free parameters. The integrand $\mathcal{L}$ depends on $w$ (and $\dot{w}$) but also on $z$. The variables $w$ and $z$ are not independent. Eliminating $\lambda$ and $y$ from the equations of motion (27), (28), and (29), we find the physical constraint

A detailed analysis of the solution of this optimization problem using Euler-Lagrange equations is given in Appendix A.

The results for both the rescaled optimal protocol $\lambda^{*}(t/t_{f})/\lambda_{i}$ and the optimal work $W^{*}/T$ depend on the dimensionless quantities

with optimized parameter $c$ and optimized final delta peak, and (iii) a protocol leading to

IV Concluding Discussion

In summary, we have calculated optimal protocols yielding the minimal mean work for underdamped Langevin dynamics in two different model potentials. Surprisingly, these optimal protocols involve jumps and delta peaks at the initial and final times $t_{i}=0$ and $t_{f}$. While we have shown that the singularities in the optimal protocol appear for harmonic potential, there is no reason to believe that this feature generically vanishes for anharmonic potentials. In fact, in the overdamped limit, a recent study has shown that initial and final jumps are also present in a simple anharmonic potential Then and Engel (2008). At first sight, such singularities seem to be unphysical since neither jumps nor delta peaks can be implemented in real experiments. Still our theoretical result is an important insight because it implies that there exists no optimal continuous protocol. Every such protocol could be improved by even steeper gradients mimicking the jumps and delta peaks at the beginning and end. If there was an experimental constraint on the allowed maximum rate of change in $\lambda$, $|\dot{\lambda}|<r$, the minimal work would still be achieved by a protocol which looks roughly like the optimal one, with the jumps and delta peaks replaced by their best approximation consistent with $|\dot{\lambda}|<r$ (e. g. steep straight lines instead of jumps). Thus, it should be possible to exploit our results for real experiments.

Our results may also be used in steered MD simulations. Even though we here have calculated optimal protocols for underdamped Langevin dynamics, there is no reason to believe that other thermostats frequently used in MD simulations would yield qualitatively different results for the optimal protocol. We have neglected memory effects by assuming white noise. While there are systems for which the underdamped Langevin equation is an appropriate physical description Douarche et al. (2006), it might still be interesting to see how our results are altered when considering memory effects.

The optimal protocol for a minimal mean work is not strictly equivalent to a protocol leading to an optimized free energy estimate. However, it has been found that the latter shares the same features (jumps at the boundaries) for overdamped Langevin dynamics Dellago (2007). Moreover, the optimal protocol leads to improved estimates of the free energy difference in both of our (underdamped) case studies. For case study I, the work distribution is Gaussian Mazonka and Jarzynski (1999); Tooru Taniguchi, E. G. D. Cohen (2008) for which it has been shown that the error in the estimate of the free energy difference decreases with decreasing mean work Gore et al. (2003). In our second case study, for which the work distribution is no longer Gaussian, the error in the estimate of the free energy difference is indeed lower for the optimal protocol compared to a linear protocol, see Fig. 4 and its caption for technical details. For thermodynamic integration, it is obvious that a minimum mean work leads to the best estimate of the free energy difference. We thus conjecture that appropriate singularities at the boundaries generically improve free energy calculations from either fast growth methods or thermodynamic integration.

For determining the optimal protocol for an unknown potential, we envisage an adaptive procedure in which trial protocols (including estimated singularities) are successively improved in an iterative fashion guided by the monitored work values. It might also prove beneficial to use the optimal moving trap protocol (case study I) rather than a linear protocol in simulations of (protein) pulling experiments.

Appendix A Solution of the optimization problem in case study II

In this appendix, we give a detailed analysis of the numerical solution of the optimization problem. In order to minimize the integral in Eq. (31), the constraint [Eq. (33)] is included in an effective Lagrangian $\mathcal{L}_{\rm{eff}}\equiv\mathcal{L}-\alpha(t)\mathcal{G}$ through a Lagrange multiplier $\alpha(t)$. Then the Euler-Lagrange equations whose solutions minimize the integral in Eq. (31) are obtained from

which, together with the constraint $\mathcal{G}=0$, define a system of three differential equations for $w$, $z$ and $\alpha$. By defining the useful new variable

we can write the initially cumbersome differential equations (37) after a tedious manipulation in the following reduced form

These equations have no analytical solution but they can easily be solved numerically for given initial conditions $w(0^{+})$, $\dot{w}(0^{+})$, $\mu(0^{+})$ and $\alpha(0^{+})$. It is important to note that some of these initial conditions are not fixed by the initial equilibrium conditions $w(0)=T/\lambda_{i}$, $\dot{w}(0)=0$, $\ddot{w}(0)=0$, $z(0)=mT$, but can be realized by additional discontinuities in the respective quantities at the boundaries. If such discontinuities do not change the value of the integral in Eq. (31), they do not affect the optimization of the integral via the Euler-Lagrange equations and hence the respective initial conditions should (in a first step) be treated as free parameters. Since the Langrangian does not depend on $\ddot{w}(t)$, discontinuities in $\dot{w}(t)$ and $\ddot{w}(t)$ can occur at the boundaries. However, a jump in the mean squared position $w(t)$ would affect the integral in Eq. (31) and thus $w(t)$ must be chosen to be continuous at the boundaries, enforcing $w(0^{+})=w(0)\equiv w_{0}=T/\lambda_{i}$. Likewise, discontinuities in $z(t)$ can occur at the boundaries. However, the initial values $z(0^{+})$ and $\dot{w}(0^{+})$ are related by the constraint $\mathcal{G}=0$. Integrating this constraint

We consider a (possible) discontinuity through the parameter $s_{1}$ in

With Eq. (44), the jump in the derivative of $w$ at the initial time as a function of $s_{1}$ becomes

In the case in which $\lambda_{i}<\lambda_{f}$, the correct sign is the negative one. Note that the last equation implies $s_{1}>1$, so that at the initial time and given the equilibrium initial distribution, it is not possible to have a decrease in the mean squared momentum. From Eqs. (44) and (38) we also find

Secondly, we define a new free parameter $s_{2}$ in

which, from the evolution equation (29), directly yields $\ddot{w}(0^{+})=\frac{2T}{m}s_{2}$. Then, writing Eq. (39) at $t=0^{+}$ and inserting the above values, the initial value of the Lagrange multiplier needed to solve the Euler-Lagrange equations is

Last, from the evolution equations (27)-(29) we find the relative value of the initial jump in the protocol as a function of $s_{1}$ and $s_{2}$:

At the end of the process, the value of $z(t_{f})$ is allowed to jump again. Recalling Eq. (43) applied now at the final time $t=t_{f}$ and isolating $z(t_{f})$, we obtain

Every quantity on the right hand side of the last equation except for $\dot{w}(t_{f})$ is fixed by the solution of the Euler-Lagrange equation. The minimum value for $z(t_{f})$, which leads to the minimal contribution to the work in Eq. (54), is reached for $\dot{w}(t_{f})\equiv s_{3}=0$.

For a comparison of the present case with its overdamped analogue Schmiedl and Seifert (2007), one can formally integrate the differential equations for $\mu$ and $\alpha$ and plug them into Eq. (39) to obtain the following integro-differential equation for $w$,

where $f(t)\equiv\frac{w(t)}{w_{0}}e^{2\gamma t/m}$ and

Combining Eqs. (27)-(29), the work $W$ [Eq. (23)] can be written as

To calculate the integral, we insert the solution of the Euler-Lagrange equations for $z$, which depends on $s_{1}$ and $s_{2}$. Then, we need to insert the boundary values for $w$ and $z$ at $t=0$ and $t=t_{f}$. In a last step, the work is optimized with respect to the free parameters $s_{1}$ and $s_{2}$.

References

Figure captions

Figure 1 : Scheme of the optimal mean position $u^{*}(t)$ and protocol $\lambda^{*}(t)$.