Free-energy estimates from nonequilibrium trajectories under varying-temperature protocols

Introduction.— The Jarzynski equality allows the calculation of free-energy differences by measuring the work done by a nonequilibrium protocol. For a system at fixed temperature $\beta^{-1}$ that is initially in equilibrium and driven out of it by the variation of a set of control parameters, the Jarzynski equality reads Jarzynski (1997, 2008)

Here $W$ is work, the angle brackets denote an average over many independent dynamical trajectories resulting from the protocol, and $\Delta F$ is the free-energy difference associated with the initial and final values of the control parameters. However, because it involves the average of an exponential, the number of trajectories required to estimate $\Delta F$ using (1) typically grows exponentially with the variance of work fluctuations Jarzynski (2006); Hummer and Szabo (2010); Rohwer et al. (2015) (see Section S1).

This problem motivates the search for nonequilibrium protocols that perform desired transformations while minimizing work or other path-extensive quantities Schmiedl and Seifert (2007); Vaikuntanathan and Jarzynski (2008); Davie et al. (2014); Blaber et al. (2021); Engel et al. (2022) ¹¹1Minimizing work does not necessarily minimize the fluctuations of work, but there is usually a strong correlation between these things. See e.g. panel (d) of Fig. 1.. However, while protocols of this nature can involve varying temperature – such as the protocol that reverses the magnetization of the Ising model with least dissipation Rotskoff and Crooks (2015); Gingrich et al. (2016) – Eq. (1) applies only at fixed temperature.

In this paper we consider a variant of (1) that applies to protocols whose temperature can vary with time. Such variants have been derived within the framework of Hamiltonian dynamics Williams et al. (2008); Chelli (2009) or Markovian dynamics specific to the Ising model Chatelain (2007). We follow Ref. Crooks (1998) and consider a general Markovian dynamics satisfying detailed balance. If we consider the protocol to involve a set of time-varying control parameters and a time-varying reciprocal temperature $\beta(t)$, with the latter starting and ending at a value $\beta$ ²²2$\beta$ with no time argument or subscript denotes the fixed reciprocal temperature at the start and end of the trajectory: we want to estimate the value $\beta\Delta F$ appearing in (1) by allowing the system at intermediate times to have a value $\beta(t)$ that may be different to the end-point value $\beta$., then (1) is replaced by

Here $\Omega\equiv\beta W+\beta Q-\Sigma$, where $Q$ is the heat exchanged with the bath and $-\Sigma$ is the path entropy produced by the trajectory (the instantaneous change of heat divided by the instantaneous temperature, summed over the trajectory). The quantity $\Omega-\beta\Delta F$ is the total entropy production ³³3Eq. (2) can be considered a special case – one where the temperatures at the trajectory endpoints are equal – of a varying-temperature version of the entropy-production fluctuation theorem Evans and Searles (2002); Seifert (2012). The angle brackets in (2) denote an average over nonequilibrium trajectories that start in equilibrium at temperature $\beta^{-1}$, that finish at the same temperature (not necessarily in equilibrium), and that otherwise involve an arbitrary change of temperature and other control parameters.

Under the same conditions, the fluctuation relation

holds. Here $P_{\rm F}(\Omega)$ denotes the probability distribution of $\Omega$ under the protocol, and $P_{\rm R}(-\Omega)$ the distribution of $-\Omega$ under the time-reversed protocol. For a fixed-temperature trajectory we have $\beta Q=\Sigma$ and so $\Omega=\beta W$, and (2) and (3) reduce to the Jarzynski equality (1) and the Crooks fluctuation relation Crooks (1999a)

respectively.

In what follows we sketch the derivation of (2) and (3). Details of the derivation are given in Sections S2–S4, which follow Ref. Crooks (1998) with minor notational changes ⁴⁴4We also show that the Jarzynski equality becomes the staged Zwanzig formula for free-energy perturbation if the trajectory remains in equilibrium, and becomes the formula for thermodynamic integration if, in addition, the control parameters change in infinitesimal increments; related limiting forms were derived in Ref. Jarzynski (1997) within the framework of Hamiltonian dynamics.. We use simulation models of a colloidal particle in an optical trap and a ferromagnet undergoing magnetization reversal to illustrate that the fluctuation relation (3) holds for varying-temperature protocols. Moreover, such protocols can give rise to smaller fluctuations of $\Omega$ and better convergence properties of (2) than do fixed-temperature protocols for $W$ and (1), allowing more accurate extraction of free-energy differences.

Sketch of derivation.— Consider a system at instantaneous temperature $\beta^{-1}(t)$. The system’s microscopic coordinates are given by a vector ${\bm{x}}$ and its energy function is $E({\bm{x}}|{\bm{\lambda}})$, where ${\bm{\lambda}}$ is a vector of control parameters. We define the protocol as the deterministic time evolution of ${\bm{\lambda}}(t)$ and $\beta(t)$. Starting in microstate ${\bm{x}}(0)={\bm{x}}_{0}$ with the parameter values ${\bm{\lambda}}(0)={\bm{\lambda}}_{0}$ and $\beta(0)=\beta_{0}$, a dynamical trajectory $\omega$ of the system consists of a series of alternating changes of the protocol, $\beta_{0},{\bm{\lambda}}_{0}\to\beta_{1},{\bm{\lambda}}_{1}\to\cdots\to\beta_{N},{\bm{\lambda}}_{N}$, and coordinates, ${\bm{x}}_{0}\to{\bm{x}}_{1}\to\cdots\to{\bm{x}}_{N}$. If $P[0\to N]$ is the probability of generating the trajectory $\omega$, and $P[0\leftarrow N]$ the probability of generating its reverse under the time-reversed protocol, then, for Markovian stochastic dynamics satisfying detailed balance, we have $P[0\to N]/P[0\leftarrow N]={\rm e}^{-\Sigma_{\omega}}$, where

is (minus) the path entropy produced by $\omega$. If we assume that $\omega$ and its reverse start in thermal equilibrium with respective control-parameter values $\beta_{0},{\bm{\lambda}}_{0}$ and $\beta_{N},{\bm{\lambda}}_{N}$, where $\beta_{0}=\beta_{N}\equiv\beta$, then the path-probability ratio can be written

Here $\Delta F$ is the Helmholtz free energy difference at temperature $\beta^{-1}$ corresponding to the change ${\bm{\lambda}}_{0}\to{\bm{\lambda}}_{N}$, and

and

are the work done and heat exchanged with the bath within $\omega$.

Summing (6) over all trajectories $\omega$ gives Eq. (2) (there we have dropped the path label subscript $\omega$ on $\Omega$). Multiplying both sizes of (6) by $\delta{\left(\Omega_{\omega}-\Omega\right)}$ and summing over trajectories $\omega$ gives Eq. (3).

Constant-temperature protocols.— Consider a model of an overdamped colloidal particle in an optical trap Schmiedl and Seifert (2007). The particle has position $x$ and the system has energy function $E(x|\lambda)=k(t)\left(x-\lambda(t)\right)^{2}/2$, where $\lambda$ specifies the trap center. The particle undergoes the Langevin dynamics

which satisfies detailed balance with respect to the system’s energy function Crooks (1999b); Dellago et al. (1998). The noise $\xi$ satisfies $\left\langle\xi(t)\right\rangle=0$ and $\left\langle\xi(t)\xi(t^{\prime})\right\rangle=2\beta(t)^{-1}\delta(t-t^{\prime})$. Initially we set $k(t)=1$ for all $t$.

We simulate (9) using the forward Euler discretization with timestep $\Delta t=10^{-3}$. Starting in equilibrium at temperature $\beta^{-1}=1$, with the trap center at $\lambda(0)=\lambda_{0}=0$, we consider the fixed-temperature protocol that moves the trap center to a final position $\lambda(t_{\rm f})=\lambda_{\rm f}=5$, in finite time $t_{\rm f}$, and that minimizes the work averaged over many realizations of the process. This protocol has the form $\lambda^{\star}(t)=\lambda_{\rm f}(t+1)/(t_{\rm f}+2)$, for $0<t<t_{\rm f}$, with jump discontinuities at the start $(t=0)$ and end $(t=t_{\rm f})$, and produces mean work $W^{\star}=\lambda_{\rm f}^{2}/(t_{\rm f}+2)$ Schmiedl and Seifert (2007). Note that $\Delta F=0$ for this protocol, because the energy function is translated but otherwise unchanged.

In Fig. 1 we verify that work distributions produced by the protocol $\lambda^{\star}$ obey the standard relations (1) and (4). In panel (a) we show the protocol $\lambda^{\star}$ for two trajectory lengths $t_{\rm f}$. In panel (b) we show for $t_{\rm f}=10$ the time-resolved distribution of particle positions $\rho(x)$ under this protocol, together with the energy function $E(x|\lambda)$ and the associated Boltzmann distribution $\rho_{0}(x|\lambda)$. Here and subsequently we calculate distributions and averages over $10^{6}$ independent trajectories. In panel (c) we verify that the work distribution produced by this protocol and its time reverse satisfy the work-fluctuation relation (4).

In panel (d) we compare protocols $\lambda^{\star}$ carried out for a range of trajectory lengths $t_{\rm f}$. Shown are the mean work $\left\langle W\right\rangle$; the variance $\sigma^{2}_{W}$ of the work distribution; the Jarzynski free-energy estimator $J_{W}=-\beta^{-1}\ln(N_{\rm traj}^{-1}\sum_{i=1}^{N_{\rm traj}}{\rm e}^{-\beta W_{i}})$, where $i$ labels trajectories and $N_{\rm traj}=10^{6}$; and the variance $\sigma^{2}_{{\rm B}_{W}}$ of the block average ${\rm B}_{W}=N_{\rm block}^{-1}\sum_{i=1}^{N_{\rm block}}{\rm e}^{-\beta W_{i}}$, where $N_{\rm block}=100$. The mean work satisfies $\left\langle W\right\rangle=W^{\star}$, as expected. The estimator $J_{W}$ returns $\Delta F=0$ for large values of $t_{\rm f}$, but for small values of $t_{\rm f}$ is imprecise; the fluctuation $\sigma^{2}_{W}$ is the source of this imprecision, and $\sigma^{2}_{{\rm B}_{W}}$ is one measure of its size.

Simple varying-temperature protocols.— In Fig. 2 we again consider the protocol $\lambda=\lambda^{\star}(t)$, for $t_{\rm f}=10$, but now set $\beta(t)=\beta^{\prime}$ for $0<t<t_{\rm f}$ (with $\beta(0)=\beta(t_{\rm f})=1$) and $k(t)=(1-t/t_{\rm f})k+(t/t_{\rm f})k^{\prime}$. The change of spring constant imposes a free-energy difference $\Delta F=\frac{1}{2}\ln(k^{\prime}/k)$: we choose $k^{\prime}=3k=3$, giving $\Delta F\approx 0.55$. We consider a range of $\beta^{\prime}$ either side of 1. In panel (a) we show that, as expected, the work fluctuation relation (4) is not obeyed for a varying-temperature protocol (here $\beta^{\prime}=0.64$), but the varying-temperature fluctuation relation (3) is. In panel (b) we we show that, as expected, the standard Jarzynski equality no longer applies: the green line is the free-energy estimator $J_{W}$, which for $\beta^{\prime}\neq 1$ does not equal $\Delta F$. By contrast, the estimator $J_{\Omega}=-\beta^{-1}\ln(N_{\rm traj}^{-1}\sum_{i=1}^{N_{\rm traj}}{\rm e}^{-\beta\Omega_{i}})$ provides a noisy estimate of $\Delta F$, indicating that (2) is obeyed.

In panels (c) and (d) we show, as a function of $\beta^{\prime}$, the variance $\sigma^{2}_{W}$ of $\Omega$ and the variance $\sigma^{2}_{{\rm B}_{\Omega}}$ of the block average ${\rm B}_{\Omega}=N_{\rm block}^{-1}\sum_{i=1}^{N_{\rm block}}{\rm e}^{-\beta\Omega_{i}}$, with $N_{\rm block}=100$. The minimum of the latter occurs around $\beta^{\prime}=0.65$, and is about a third of that at $\beta^{\prime}=1$ (where $\Omega=\beta W$ and the relations (2) and (3) reduce to (1) and (4)). While fluctuations tend to increase with temperature, the combination $\Omega\equiv\beta W+\beta Q-\Sigma$ and its fluctuations can be smaller than $W$ and its fluctuations. These effects compete, and for some range of $\beta^{\prime}\neq 1$ the fluctuations of $\Omega$ are smaller than those of $W$ at $\beta^{\prime}=1$, leading to better convergence of (2) than (1).

Varying-temperature protocols in the presence of a phase transition.— This difference in convergence properties is small, because the physics of the trap model does not change substantially with temperature. But when a system’s physics does change with temperature, the difference between the convergence properties of the fixed- and varying-temperature relations can be significant.

In Fig. 3 we consider the 2D ferromagnetic Ising model with fixed coupling $J=1$, on a $32\times 32$ lattice, simulated using Glauber dynamics for $t_{\rm f}=10^{3}$ Monte Carlo sweeps. Protocols start and end at parameter values $\beta(0)=\beta(t_{\rm f})=1$ and $h(0)=-h(t_{\rm f})=-1$, giving $\Delta F=0$, in which case $\Omega$ is equal to the total entropy production. We determine time-dependent protocols $\beta(t)$ and $h(t)$ for $0<t<t_{\rm f}$ by expressing these quantities using a deep neural network and training that network by genetic algorithm to minimize $\left\langle\Omega\right\rangle$ measured over $10^{4}$ independent trajectories Whitelam (2023). Protocols learned in this way effect magnetization reversal. We consider one set of learning simulations with the constraint that $\beta(t)$ remain constant, and another in which $\beta(t)$ is allowed to vary.

In panel (a) we show in parametric form the protocols learned in this way. The varying-temperature protocol avoids the critical point and the first-order phase transition line Rotskoff and Crooks (2015); Gingrich et al. (2016), while the fixed-temperature protocol is constrained to cross it. In panel (b) we show protocols as functions of time, together with time-ordered snapshots taken from typical trajectories. The fixed-temperature protocol effects nucleation and growth, accompanied by large values of dissipation: $\left\langle\Omega\right\rangle=\left\langle\beta W\right\rangle\approx 1550$. By contrast, the varying-temperature protocol is accompanied by much smaller dissipation, with $\left\langle\Omega\right\rangle\approx 11$.

As a result, distributions of work for the fixed-temperature protocol and its time-reverse are well-separated, while distributions of $\Omega$ for the varying-temperature protocol and its time reverse cross at a value that approximates $\Delta F=0$; see panel (c). Using $10^{5}$ trajectories, the free-energy estimator $J_{\Omega}$ yields a value 0.46, with block-average variance $\sigma^{2}_{{\rm B}_{\Omega}}\approx 13,700$. Thus (3) and (2) provide an accurate if imprecise measure of $\Delta F$. By contrast, if we are constrained to fixed temperature in order to apply the standard relations (1) and (4), measuring $\Delta F$ is not a realistic proposition. Fixed-temperature trajectories hundreds of times longer would be required to achieve convergence properties comparable to (2) and (3) using varying-temperature trajectories ⁵⁵5For the Ising model we could change $J$ at fixed $\beta$ to mimic temperature variation. Our model study is carried out to represent an experiment in which we cannot change the microscopic parameters of a material, but we can use temperature as a control parameter to take us across a phase boundary..

Conclusions.— We have shown, within the framework of Markovian stochastic dynamics satisfying detailed balance, that the relations (2) and (3) replace the Jarzynski equality (1) and Crooks work fluctuation relation (4), for trajectories influenced by a time-varying temperature that starts and ends at a value $\beta^{-1}$. We have used simulation models to show that free-energy differences can be calculated more accurately using the varying-temperature relations than the fixed-temperature ones, particularly when varying temperature gives us the freedom to avoid the large dissipation associated with a first-order phase transition. To measure the quantity $\Omega=\beta W+\beta Q-\Sigma$ that appears in (2) and (3) we must be able to measure the time-resolved heat flow, as well as total heat and work. For many experiments this will be technically demanding, but it is possible in principle.

Acknowledgments.— I thank Corneel Casert for comments on the paper. Code for the trap model can be found here tra . Code for doing neuroevolutionary learning of Ising model protocols can be found here dem , which accompanies Ref. Whitelam (2023) (for this paper we neglect the shear term, Eq. (7), used in that work). This work was performed at the Molecular Foundry at Lawrence Berkeley National Laboratory, supported by the Office of Basic Energy Sciences of the U.S. Department of Energy under Contract No. DE-AC02–05CH11231.