Convergence of thermodynamic quantities and work fluctuation theorems in presence of random protocols
Rahul Marathe, Sourabh Lahiri

TL;DR
This paper investigates how random variations in protocols affect the convergence of thermodynamic quantities and work fluctuation theorems in non-equilibrium systems, using analytical and numerical methods on Brownian particles.
Contribution
It introduces the analysis of random protocols in fluctuation theorems, extending the understanding of convergence and symmetry functions beyond fixed protocols.
Findings
Random protocols often do not affect convergence compared to fixed protocols.
Analytical and numerical results support the robustness of fluctuation theorems under protocol randomness.
Results are applicable to experimental verification in microscopic systems.
Abstract
Recently many results namely the Fluctuation theorems (FT), have been discovered for systems arbitrarily away from equilibrium. Many of these relations have been experimentally tested. The system under consideration is usually driven out of equilibrium by an external time-dependent parameter which follows a particular {\it protocol}. One needs to perform several iterations of the same experiment in order to find statistically relevant results. Since the systems are microscopic, fluctuations dominate. Studying the convergence of relevant thermodynamics quantities with number of realizations is also important as it gives a rough estimate of number of iterations one needs to perform. In each iteration the protocol follows a predetermined {\it identical/fixed} form. However, the protocol itself may be prone to fluctuations. In this work we are interested in looking at a simple…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Convergence of thermodynamic quantities and work fluctuation theorems in presence of random protocols
Rahul Marathe
Department of Physics, Indian Institute of Technology, Delhi, Hauz Khas 110016, New Delhi, India.
Sourabh Lahiri
Department of Physics, Birla Institute of Technology Mesra, Ranchi, Jharkhand 835215, India.
Abstract
Recently many results namely the Fluctuation theorems (FT), have been discovered for systems arbitrarily away from equilibrium. Many of these relations have been experimentally tested. The system under consideration is usually driven out of equilibrium by an external time-dependent parameter which follows a particular protocol. One needs to perform several iterations of the same experiment in order to find statistically relevant results. Since the systems are microscopic, fluctuations dominate. Studying the convergence of relevant thermodynamics quantities with number of realizations is also important as it gives a rough estimate of number of iterations one needs to perform. In each iteration the protocol follows a predetermined identical/fixed form. However, the protocol itself may be prone to fluctuations. In this work we are interested in looking at a simple non-equilibrium system namely a Brownian particle trapped in a harmonic potential. The center of the trap is then dragged according to a protocol. We however lift the condition of fixed protocol. In our case the protocol in each realization is different. We consider one of the parameters of the protocol as a random variable, chosen from some known distribution. We study the systems analytically as well as numerically. We specifically study the convergence of the average work and free energy difference with number of realizations. Interestingly, in several cases, randomness in the protocol does not seem to affect the convergence when compared to fixed protocol results. We study symmetry functions. A Brownian particle in a double well potential is also studied numerically. We believe that our results can be experimentally verified.
I Introduction
The field of non-equilibrium statistical mechanics, especially when the system is far from equilibrium, has attracted a lot of works in the last two decades or so. It is of great practical importance, given the fact that in nature, one almost always deals with irreversible or non-equilibrium processes. The framework of equilibrium thermodynamics is not so useful to study such processes, and we need separate recipes for dealing with them. One of the major developments in this area has been a group of relations collectively called “Fluctuation Theorems”. They constitute relations that have the rare distinction of being exact equalities, no matter how far the system has been driven away from equilibrium. Let the system be initially at equilibrium with a thermal bath of inverse temperature at a given initial value of an externally controlled parameter or protocol . At time , the time dependence of this parameter is switched on, as a result work is done on the system. Such a process is called forward process, in contrast to the reverse process, in which the initial value of parameter is at which the system is at thermal equilibrium, and thereafter the reverse protocol is applied for the same time of observation . Two well-known equalities hold for such processes:
[TABLE]
Here, is the stochastic work done along individual trajectories, and the average (denoted by angular brackets) is the ensemble average. The first equality was proved by C. Jarzynski jar97_prl and the second by G. E. Crooks cro98_jsp ; cro99_pre .
A large volume of works followed thereafter. In the meantime, Sekimoto established the field of Stochastic Thermodynamics on a firm ground, thereby providing clear definitions for thermodynamic quantities along individual phase space trajectories (i.e. in individual experiments) sek98 . Fluctuation theorems for total entropy production were derived by Seifert sei05_prl ; sei08_epjb . Fluctuation theorems for systems undergoing transitions from one non-equilibrium steady state to another were studied hat01_prl ; spe07 . FTs have been derived for systems following Hamiltonian as well as stochastic dynamics sei12_rpp . They have been proved for both closed and open quantum systems han11_rmp .
Our focus in this article would however be on the convergence of thermodynamic quantities like work, free-energy and Jarzynski and the Crooks work theorems, Eq. (1) and Eq. (2). Till now, all the works have studied FTs with a fixed protocol. In other words, although we are repeating the experiment a large number of times, the functional form is the same for all these experiments. The stochasticity of work arises entirely because the system traces out a different trajectory in phase space every time because of thermal noise, and not because of variations in protocol.
We, in this work, would like to lift the condition of deterministic protocol and check whether the theorems are robust even if the protocol is different in different experiments of the ensemble. The parameters used in the protocols will be sampled from some distribution. Our objective would be to compare the effect of three such distributions: the uniform, Gaussian and exponential, and to numerically observe which option provides better convergence to analytical values. In fact, it is important to state that we use the same number of realizations of the experiment as in the case of a fixed protocol to obtain the value of .
Our protocols will mainly consist of dragging a colloidal particle through a medium, using different types (time dependence) of protocols as well as using different random distributions of the parameters involved for the same functional form of the protocol. Several experimental and theoretical works have been performed using dragged colloidal particles. Wang et. al. performed an experiment wang02_prl with a dragged colloidal particle, where he demonstrated the validity of the Crooks work fluctuation theorem. Subsequently, the authors of zon03_pre explicitly showed the validity of the Crooks fluctuation theorem and the violation of the steady state fluctuation theorem in such systems. The model has also been studied in detail in tre04_pnas ; jar99_arxiv ; spe05_epjb ; lah09_pre .
We have also used sinusoidal protocols for two examples where the work distribution becomes non-Gaussian.
II The model
We will consider a particle in a medium at temperature , that is placed in a harmonic trap of stiffness constant and follows the overdamped Langevin equation of motion:
[TABLE]
Here, is the usual Gaussian white noise, having the properties and , is the friction coefficient of the medium, and overhead dot implies total time derivative. The center of the harmonic trap is dragged according to a protocol . For the analysis carried below we fix the values of parameters to be , , . To arrive at these dimensionless parameters, the time has been scaled by the relaxation time , the position by the square root of initial thermal width , and energies are measured in units of , i.e. the following replacements have been made:
[TABLE]
The values of these parameters in any system of units can therefore be obtained by multiplying the dimensionless values by the appropriate factors listed above. The typical dimension of position at mesoscopic scales is nanometers, of time is seconds and of energy is piconewton-nanometers.
The system is initially allowed to equilibrate at (the initial value of the dragging protocol) and then the external time-dependent perturbation is switched on. Note that for equilibrium initial distributions, the ergodicity of the system prevails throughout the time of observation, so that in simulations one can use time averaging if that is computationally less expensive as compared to ensemble averaging tha18_pre , a fact that can often come handy in simulations. The final value of the perturbation (dragging protocol), , will be chosen from a random distribution, thus making the protocols random in nature.
We consider two forms of the protocols:
Ramp protocol: varies linearly from at time to the final value at time . The form is given by
[TABLE]
The parameter is a random number, sampled from a given distribution, while without loss of generality, we chose . 2. 2.
Optimal protocol: The other form that we consider is the optimal protocol described in sei07_prl . This protocol is important since it minimizes the work done between the given initial and the final protocol values.
[TABLE]
As is obvious, the parameter does not reach its final value at time , nor does its initial value correspond to . This means that the optimal protocol requires a jump in the value of external parameter at the initial and the final times. Once again, , while is a random number chosen from a distribution.
We have chosen the most probable value of same as the final value reached in the fixed protocol case (with which we have provided comparisons). Throughout this article, we choose this value to be .
Also, since for any form of the dragging protocol, we do not need to keep the initial and final values of the parameter fixed for all protocols. Here, is a constant for a particular experiment, but varies randomly from one experiment to another. The random values, as mentioned earlier, will be sampled from uniform, Gaussian and exponential distribution and the resulting outcomes would be compared. The initial value of the parameter is always taken to be zero.
In the case of optimal protocol, where the initial and final values of the parameter are given, there must be jumps at initial and final times as shown in sei07_prl . This is not required in the ramp protocol, where the initial and final values of the protocol are connected simply by a straight line (obviously, this is not an optimal protocol, since we are not allowing for the jumps to occur). However, in either case, as long as protocol is a simple dragging of the center of the trap, the free energy remains zero and the functional form of the protocol is immaterial. Fig. (1) shows an example of such protocols.
To study how the convergence of the Jarzynski equality is affected when the protocol is randomized, we write the Jarzynski equality as
[TABLE]
Here, the index labels the individual realizations. In practice, the relation “converges” if is large enough, so that we can obtain an accurate value (up to the allowed tolerance level) of , that agrees with the “correct” value that emerges in the limit (which can be calculated analytically in some simple cases). We would use different random distributions to sample the parameters of our protocol and calculate the minimum value of at which the relation converges. The distribution corresponding to which is smaller, is practically more useful, since we can obtain results in quicker time.
We will further find the convergence of the mean work obtained from simulation with the theoretically predicted value of mean work. In the next section, we calculate the theoretical expressions for the mean works with which we wish to compare the numerically obtained values.
II.1 Ramp protocol
For the ramp form of protocol, we first note that the work done will be given by sek98
[TABLE]
where the angular brackets , denote the average over thermal noise.
Let us choose the initial value of the parameter to be fixed: , and the form of the protocol to be . The parameter is sampled from a random distribution, which we choose to be either uniform, Gaussian or exponential:
[TABLE]
For all distributions we choose as mentioned earlier. For simplicity, we have chosen the width of the Gaussian distribution , to be unity. Similarly, we have chosen the domain of the uniform distribution to be unity, so that the sampled random numbers are in . The value of depends on the specific example at hand. In our case, , so that the mean is 5, which is also the value of for the fixed protocol. We also note that for the exponential distribution, mean is equal to the standard deviation. The mean position is given by (here represents averaging over the thermal noise, whereas the overbar represents average over .)
[TABLE]
The mean work is given by (see Eq.7)
[TABLE]
As mentioned earlier the protocol does not smoothly extend to the given initial and final values of the protocol (which are fixed), so there are sudden jumps at the ends in order to abide by this constraint. This introduces an added contribution to the work. However, in the above case (dragging the center of trap), the value of equilibrium change in free energy identically vanishes, irrespective of the initial and final values of protocol. So we do not adhere to a fixed value of in our analysis.
II.2 Optimal protocol
If we use the protocol used in sei07_prl , then it is necessary to fix the final value of the parameter, a principle on which the derivation of the form of optimal protocol is based. In this case, the work done is given by
[TABLE]
where is the contribution coming due to the jumps at the initial and final points of the protocol. Below we calculate these contributions. At the jumps, the work done on the system is given by the change in internal energy (since no heat is dissipated in a sudden quench). At the initial time , we then have
[TABLE]
Using the fact that , averaging over the noise, and using , the right hand side becomes simply . Now, integrating over the protocol distribution, we get
[TABLE]
Similarly, at the final time ,
[TABLE]
where we have used (see sei07_prl ). In sei07_prl , it was also shown that the position of the particle, when averaged over thermal noise, is given by . Therefore, after averaging over thermal noise, the right hand side becomes . Thus
[TABLE]
Finally, adding both the contributions we get
[TABLE]
Setting in accordance with our choice, the total contribution of jumps, when averaged over noise as well as is given by
[TABLE]
thus
[TABLE]
We also note that for both the cases discussed above the work distributions are Gaussian for fixed protocols, and they satisfy Eqs. (1) and (2). Then it is easy to show that the variance . However in case of protocols chosen from a distribution , work distributions are non-Gaussian and can be calculated using the Bayes’ theorem:
[TABLE]
This in general can not be integrated to get a closed form. But still satisfies Eqs. (1) and (2).
III Numerical simulations
III.1 Uncorrelated protocol noise
We first study the case of dragging protocol, where the equation of motion is given by Eq. (3). We have used the velocity Verlet algorithm by Ermak ermak ; allen_tildeslay to integrate our equations of motion with time step , and have generated trajectories to construct a full ensemble.
Using the Jarzynski equality, , one can calculate for each case. Ideally, an infinitely large number of trajectories need to be considered for taking the average appearing on the left hand side, but practically for a large enough number of trajectories we can obtain a convergence that is consistent with our required tolerance. Thus, if for a given protocol, if we find the value so obtained as a function of the ensemble size , then an increase in would lead to a more accurate result.
If the analytical value of is known (which, for the dragging protocols, are zero), then we can find the required number at which the numerically obtained value falls within the tolerance range. A protocol that requires very high is much less efficient (in terms of convergence) when compared to one that provides convergence much faster. We set and throughout. Similarly, one can compute the number of trajectories required for the convergence of mean work, whose expressions are given by equations Eq. (12) and Eq. (20).
In figure 2, we show this behavior when the protocol is of type (1), but the parameter is either fixed, or sampled from a uniform and a Gaussian distribution. When it is random, we set the parameters of the distribution such that . In the case of uniformly distributed , we adjust the range of the distribution to , so that the mean value is . For the normally distributed , we simply set the peak of the distribution at 5. The theoretical mean values of work are given by (fixed protocol), 1.191 (uniform distribution) and 1.235 (Gaussian distribution).
In figure 3(a), where the convergence of mean work has been tested for protocol (2) i.e. the optimal protocol, the scenario for the convergence to theoretical value remains the same as in the case of ramp protocol, i.e. the values of are comparable in all the three cases and are roughly . The theoretical values of are 1.136 (fixed protocol), 1.140 (uniform distribution) and 1.181 (Gaussian distribution). However, as we observe from figure 3(b) for , we find that the convergences for uniform and Gaussian distributions become appreciable only beyond . We will observe later (see figure 6) that plotting the symmetry function is a better way to find efficiently.
Figure 4(a) and (b) show the convergence of mean work and free energy for protocol (1), respectively, when is exponentially distributed. Again, the convergence of mean work is appreciable roughly beyond . The convergence of free energy change is poor, even for , an obvious consequence of the large variance () as compared to uniform and Gaussian distributions. The solid line in fig. 4(a) corresponds to the analytical value . Figure 5 (a) and (b) displays the similar plots for protocol (2), for exponentially distributed . This time, the convergence in mean work appears to be quicker (figure 5(a)), although that in the free energy change (figure 5(b)) is again quite inefficient. The solid line in figure 5(a) corresponds to the analytical value .
Finally, we turn our attention to the validity of Crooks Fluctuation Theorem for work (see Eq. (2)). This can be done by taking the logarithm of Eq. (2):
[TABLE]
The quantity on the left hand side is called the symmetry function, and according to the above relation, if CFT is valid, must yield a straight line of slope unity as a function of .
In figure 6, we have plotted the symmetry functions for uniform distribution and Gaussian distributions corresponding to each type of protocol: ramp (protocol (1), figure 6(a)) and optimal (protocol (2), figure 6(b)). In (a), we again find excellent agreement for uniform protocol. Even for the Gaussian protocol, we find very good agreement. In (b), we again find very good agreement for the protocol (2) having uniformly distributed . The Gaussian, however, shows deviations beyond , which can be accounted for by the poor convergence of protocol (2) for large values of work (sampling of rare trajectories becomes difficult). Note that the convergence is not a problem for uniform protocol, since the values of are bounded.
IV Numerical results when work distributions are non-Gaussian even for deterministic protocol
IV.1 Numerical results for time-dependent trap strength
In all the cases discussed above the work distributions are Gaussian. Here we give numerical results for a potential of the form , where , is a half sinusoidal cycle. The equation of motion then becomes
[TABLE]
The parameter is chosen from a given distribution as earlier. Unlike earlier, in this case to make sure that for the full duration of the cycle. We choose . For uniform distribution we choose . For Gaussian distribution we choose both mean and variance to be . Note that for Gaussian distribution we also need to take to make sure for the full duration of the cycle. The work done is given by
[TABLE]
which is clearly non-Gaussian, being quadratic in , even for fixed .
Figure 7 shows the convergences for free energy change (which is again zero, since identically). It is clear that even with such a non-Gaussian distribution of work, the convergences are nearly same as in the case of the corresponding deterministic protocol.
We perform a more stringent verification of our results in the next section, where we replace the harmonic trap by a bistable potential.
IV.2 Numerical results for bistable potential
In order to check that the efficiency of random protocol is not a consequence of the simplistic form of quadratic potential, we now numerically simulate a bistable or double-well potential perturbed by a sinusoidal drive, of the form . The Langevin equation now is
[TABLE]
The expression for work becomes
[TABLE]
Here, again the work distribution is clearly non-Gaussian in nature since the distribution of itself is non-Gaussian, even for fixed .
The perturbation is carried out for one full time-period of the drive: , so that in the process. The results are shown in figure 8. We clearly observe that the blue curve (where is distributed normally with variance of unity, with the mean value fixed at , being the value of the parameter for a deterministic protocol) takes an extremely long time to converge, and as such is highly inefficient. Nevertheless, the uniform distribution does almost as good a job as the fixed one.
V Heuristic explanation for the observed behavior
If we consider the averaging over protocols and trajectories as a double integral, we might naively expect that the convergence would require about realizations of the experiment, where is the number of trajectories required for convergence with a fixed protocol and is the size of the ensemble of protocols. Thus, if realizations are needed for convergence under a fixed protocol, realizations should be required for random protocols. However, we find that this is never the case, and convergence appears much earlier. In our opinion, the following reason accounts for this apparently unusual result. In general, all values of sampled from a particular distribution need not require equal orders of . For instance, in the Gaussian distribution, the values of that are closer to would represent realizations where the protocol is close to being a fixed one, so that the convergence would occur earlier. On the other hand, those appearing at the tails would require an that is higher by an order of magnitude or more. However, since the weightage of these “rare values” is much less as compared to the ones that have values of comparatively closer to , the value of must be reduced considerably. Note that this happens clearly for the uniform and the Gaussian case, where the variances are equal to 0.083 and 1 respectively, which are small compared to . For the exponential distribution, however, the most probable value of is appreciably different from its mean. Thus, a major contribution comes from values of that are far from , and due to this reason it is difficult to obtain convergence within the size of the ensemble that we have been able to simulate.
Our analysis corroborates this heuristic argument. We note that the approach could be useful, for example, when there is an inherent error in the control of the external parameter, so that to some extent the protocol gets randomized. The nature of work distributions have been shown in figure 9, where we find that the Gaussian randomness in protocol gives rise to an exponential distribution, the uniform one leads to faster decay and the exponential one leads to slower decay. This again shows that an exponential randomness in requires a larger ensemble for convergence. In the Gaussian case, however, the rate of convergence will become slow if the variance is large as compared to the mean, which happens when the rate parameter in (note that this can be different on either sides of the peak) becomes much smaller than unity sab11_epl .
VI Relation to incorrect work measurement
We can also view the problem in the light of lah16_pre , where instead of randomizing the protocol, there were random errors in the measurement of work itself. In lah16_pre it was shown that even if the work done in a process is measured inaccurately (might be due to low resolution of the apparatus measuring the particle’s position), in several situations that occur frequently in experiments one can deduce the correct value of free energy change, using the symmetry functions obtained from the Crooks’ work fluctuation theorem. We will show that the work done in our case can be expressed as the work done using a fixed protocol, and an error term (see eq. (29) below). However, the error term turns out to be non-Gaussian in nature, which produces a more complex situation as compared to lah16_pre , whose analytical techniques would be useful only if one of the terms in the error becomes negligible as compared to the other terms, as detailed below.
Let us consider the ramp protocol (), where values are normally distributed. The work done is given by
[TABLE]
Now, we can define the deviation of from its mean as , so that , so that
[TABLE]
where gives the last three terms on the RHS. The above relation provides the general expression for work, irrespective of the distribution of protocol. Here,
[TABLE]
Note that is not the work done by a fixed protocol, because of the dependence of on the protocol. Nevertheless, the solution for the Langevin equation Eq. (3) is
[TABLE]
Here, is the value of the position of the particle at , which follows a Boltzmann distribution. Therefore, can now be related to the work done for fixed protocol:
[TABLE]
Collecting all results, we have
[TABLE]
where .
However, if is Gaussian (of zero mean in the present case), then we find that the last term in is non-Gaussian. As a result, a direct connection with lah16_pre is not on offer. However, if the variance is small enough, this term can be ignored, and we should recover the results of lah16_pre .
VII Conclusions
The work fluctuation theorems use deterministic protocols that are exactly the same for all experimental realizations constituting the ensemble. In this work, we study how the results are affected when the protocol is not deterministic but has a random part, so that their time dependence is different for different experimental realizations. We find that if the initial and final values of the externally manipulated parameter are held fixed, then irrespective of the type of randomness incorporated, the Jarzynski Equality is satisfied when the ensemble is large enough. We find that even experiments (number of realization needed for convergence when a deterministic or fixed protocol is applied) suffice to provide a good convergence of our results, which is surprising.
To check the convergence of Jarzynski Equality, we have studied a dragged harmonic trap, given by Eq. (3), where the equilibrium free energy is known to remain constant (). As a further check, we have analytically computed the mean works done in all the cases, and have observed the convergence of obtained from numerics with the theoretical values, as a function of the ensemble size.
We have numerically tested our results in a couple of examples where the work distribution follows a non-Gaussian distribution. These cases include a harmonic trap of varying stiffness, and a bistable potential. In both cases, we find that the convergence rates with number of realizations are comparable.
We also believe that this study can be extended to steady state fluctuation theorems. This work is currently underway.
VIII Acknowledgment
SL thanks DST-SERB, India (sanction order no. ECR/2017/002607) for financial support. RM thanks C. Jarzynski for initial discussions on this topic. Authors also thank D. Lacoste for careful reading of the manuscript.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1(1) C. Jarzynski, Phys. Rev. Lett. 78 , 2690 (1997).
- 2(2) G. E. Crooks, J. Stat. Phys. 90 , 1481 (1998).
- 3(3) G. E. Crooks, Phys. Rev. E 60 , 2721 (1999).
- 4(4) K. Sekimoto, Prog. Theor. Phys. Supp. 130 , 17 (1998).
- 5(5) Udo Seifert, Phys. Rev. Lett. 95 , 040602 (2005).
- 6(6) U. Seifert, Eur. Phys. J. B 64 , 423 (2008).
- 7(7) T. Hatano and S. Sasa, Phys. Rev. Lett. 86 , 3463 (2001).
- 8(8) T. Speck, V. Blickle, C. Bechinger, and U. Seifert, Europhys. Lett. 79 , 30002 (2007).
