Suppression of work fluctuations by optimal control: An approach based on Jarzynski’s equality
Abstract
Understanding and manipulating work fluctuations in microscale and nanoscale systems are of both fundamental and practical interest. For example, aspects of work fluctuations will be an important factor in designing nanoscale heat engines. In this work, an optimal control approach directly exploiting Jarzynski’s equality is proposed to effectively suppress the fluctuations in the work statistics, for systems (initially at thermal equilibrium) subject to a work protocol but isolated from a bath during the protocol. The control strategy is to minimize the deviations of individual values of from their ensemble average given by , where is the work, is the inverse temperature, and is the free energy difference between two equilibrium states. It is further shown that even when the system Hamiltonian is not fully known, it is still possible to suppress work fluctuations through a feedback loop, by refining the control target function on the fly through Jarzynski’s equality itself. Numerical experiments are based on linear and nonlinear parametric oscillators. Optimal control results for linear parametric oscillators are also benchmarked with early results based on accelerated adiabatic processes.
I Introduction
Aspects of thermodynamics at the nanoscale have attracted great attention in recent years. In a small system with few number of particles, thermal fluctuations and quantum fluctuations in work, heat, and other quantities can be comparable to their ensemble mean values. Understanding and manipulating these fluctuations hence become an interesting and important topic. From a fundamental point of view, the connections between nonequilibrium statistical fluctuations with equilibrium properties, such as those established through fluctuation theorems (e.g., Jarzynski equality Jarzynski.97.PRL (); Jarzynski.97.PRE () and crooks theorem crooks ()), have laid a solid and fruitful foundation for nanoscale thermodynamics Campisi.11.RMP (). It is interesting to note that the existence of such type of fluctuation theorems does not rule out the possibility to control statistical fluctuations. That is, under appropriate conditions, it is still possible to manipulate nonequilibrium processes and alter the statistical fluctuations in thermodynamic quantities, without violating any fundamental law in thermodynamics or modifying the fluctuation theorems themselves. In the case of Jarzynski’s equality, i.e., Jarzynski.97.PRL (); Jarzynski.97.PRE (), where is work done on a system initially at thermal equilibrium with inverse temperature and is the free energy difference between the final and initial equilibrium states of the same , a change in fluctuations of can have an impact on how fast the statistics of may converge towards . This is the first motivation of this study.
Advances in our knowledge of thermodynamics of small systems have stimulated studies of efficient energy devices (classical or quantum) at microscale and nanoscale, sometimes consisting of few particles or few degrees of freedom. Of particular interest here is that several designs of microscale or nanoscale heat engines have been proposed theoretically Abah.12.PRL (); Bergenfeldt.14.PRL (); Zhang.14.PRL (). For this type of applications where heattowork conversion efficiency and the power of work output are apparently crucial, the stability or reliability of work output also becomes an important performance indicator Deng.13.PRE (); dario (). This is because the cycletocycle fluctuations in the work output are an unavoidable characteristic of nanoscale heat engines. Given two heat engines with the same mean work output per cycle, it seems desirable to prefer the one with less fluctuations and hence more uniform output. This understanding further motivates us to ask the following question: how to systematically suppress work fluctuations in a given protocol (such as in one step of a heat engine cycle that does not involve any heat exchange with a bath)?
It has been realized that work fluctuations in some controlled processes can be relatively smaller than those in bare uncontrolled processes. In particular, work fluctuations in an adiabatic process (here “adiabatic” means very slow, as compared with the system’s own natural time scale) are relatively small Deng.13.PRE (). However, an adiabatic process in the standard sense is simply too slow, and as such the power of work output based on adiabatic processes vanishes. For this reason it is natural to consider accelerated adiabatic processes to suppress work fluctuations and increase the power of work output at the same time Deng.13.PRE (); campoarxiv (). Certainly, accelerated adiabatic processes just represent a special type of controlled dynamics, and extensions based on alternative constructions of a control field should be possible. This paper is precisely one extension of Ref. Deng.13.PRE () in order to suppress work fluctuations in a wider variety of systems. The extension is based on an optimal control framework Shi.90.JCP (); Peirce.88.PRA (); Shi.91.CPC (); Sundermann.99.JCP (); Ohtsuki.04.JCP () that is far more versatile and far more flexible. We shall also use known results based on accelerated adiabatic processes to benchmark our optimal control results.
To see the necessity and advantages of an optimalcontrolbased extension, let us first make brief comments on early studies of accelerated adiabatic processes Demirplak.08.JCP (); Demirplak.05.JPCB (); Demirplak.03.JPCA (); Berry.09.JPAMT (); Masuda.10.PRSLSA (); Chen.10.PRL (); Chen.10.PRLa (); Choi.11.PRA (); Schaff.11.NJP (); Bason.12.NP (); Zhang.13.PRL (); Jarzynski.13.PRA (); Campo.11.PRA (). For the sake of discussions later on, we divide accelerated adiabatic processes known to date into two types (though these two types can even be regarded as being equivalent upon a transformation Deffner.14.PRX () and there is no clear distinction in the literature). In the first type, an additional control Hamiltonian is introduced to drive a system (within a short time scale), such that the evolution of the system, either classical or quantum, still follows the adiabatic evolution of the original bare system. Here we call this type of control as fastforward adiabatic driving (FFAD) Masuda.10.PRSLSA (). The extra control Hamiltonian in realizing FFAD may be found for very simple systems but in general, its analytical form is not available and and it is also challenging to precisely implement it experimentally. The second type of accelerated adiabatic processes, which we call “shortcuts to adiabaticity” (STA) mugareview (), relies on special time dependence of system parameters. In this sense, the required form of the control Hamiltonian is already contained in the original bare Hamiltonian. Within a finite time, the system under STA control can reach the same final states as those reached in conventional adiabatic processes. However, during an STA process the system may make nonadiabatic transitions Tu.14.PRE (); campoarxiv (). The realization of STA has been established in several simple systems such as twolevel systems or linear parametric oscillators. As Ref. Deffner.14.PRX () showed, STA in known examples can be understood as a result of a canonical transformation of FFAD (or vice versa). Though Campo and Jarzynski extended STA to more complicated systems such as manybody systems Campo.13.PRL (), a general discussion by them also indicates that this method can only apply to systems with the socalled “scaleinvariant driving” Deffner.14.PRX (). For these reasons both FFAD and STA require full knowledge of the system Hamiltonian, and this requirement presents a limitation when we attempt to use accelerated adiabatic processes to suppress work fluctuations in general situations.
In this work we restrict ourselves to classical systems (but our optimal control approach can be easily extended to quantum systems). As such only work fluctuations across a thermal ensemble are to be suppressed. Our classical results should be also useful in guiding possible manipulations of work fluctuations in quantum systems, especially in cases with a relatively high temperature. Specifically, we develop a method to design control fields to suppress work fluctuations based on the wellknown optimal control theory (OCT). The peculiarities of our OCT here lie in two aspects. First, we need to handle a thermal ensemble in order to consider work fluctuations. Thus our control is ensemblebased optimal control. Second, we need to design a useful control target function in order to reach our goal. As shown below, our control target function directly exploits Jarzynski’s equality, and the constructed target function is biased against deviations of individual values of from their ensemble average, i.e., . We show that our optimal control method built in this manner can be applied to a wide variety of systems, including highly nonlinear systems. When it is applied to previously known simple systems, it has essentially the same quantitative performance as that based on accelerated adiabatic processes, but now the control fields have many possible forms. Furthermore, our OCT method can be extended to implement a feedback mechanism, so as to suppress work fluctuations in systems with unknown system parameters. This is done by guessing the control target function iteratively through Jarzynski’s equality again.
This paper is arranged as follow. In Sec. II, we apply OCT to the topic of suppression of classical work fluctuations, with necessary details. In Sec. III, we test our approach using a simple system and compare our results with those obtained based on accelerated adiabatic processes. In Sec. IV, we consider a nonlinear oscillator and show that our OCT approach can still operate well. In Sec. V, we consider a nonlinear oscillator with some unknown system parameters. We show that, in the spirit of feedback optimal control, it is possible to refine the control target function on the fly, and then in the end it is still possible to find a control field to suppress the work fluctuations. The last section of this paper gives a brief summary. An appendix regarding numerical implementation of our OCT is also given.
Ii Theory of Optimal Control of work fluctuations: use of Jarzynski’s equality to construct a control target function
Let us consider a general timedependent classical system with Hamiltonian , where represents phase space coordinates (we assume onedimensional problems), and is a timedependent parameter due to a work protocol starting from to . For example, this can represent the frequency of a parametric oscillator, the length of a pendulum, or the width of a box, etc. Such a classical system is initially prepared at thermal equilibrium, with the phase space probability density given by
(1) 
where is the partition function associated with . To suppress work fluctuations in this protocol defined by , an extra control field is applied. Then the total Hamiltonian of the system becomes
(2) 
where is the control Hamiltonian and is assumed to be the timedependent amplitude of a control field. The evolution of obeys the following Hamilton’s equation of motion: and .
Next we consider a certain physical quantity evaluated at , which is written as . Consistent with this notation and for convenience, we also define and from now on. The minimization of the ensembleaveraged value of will be important in our OCT. The thermal ensemble average of is given by
(3) 
where represents the entire phase space. We call defined above a target function since it will be a quantity we hope to minimize. For a control problem, typically a cost function is also needed to reflect a costrelated constraint. Here we define the cost function using the amplitude of the control field, i.e.,
(4) 
where is a weightage factor that may depend on time. The larger is, the more constraints posted on the control field due to cost considerations. The overall target function can then be defined as . That is, the problem becomes to minimize under the general dynamical constraints reflected by Hamilton’s equations of motion. To that end we introduce two Lagrange multiplier as functions of , denoted by and . We then minimize instead, with
(5) 
Let be the variation in due to , an arbitrary variation in , . Then the variation in due to is found to be
(6) 
To minimize we let . Since the variation is arbitrary, one has the following relations:
(7) 
The list of relations in Eq. (7) can be numerically solved by an iteration procedure and some details are presented in Appendix.
After outlining the general steps in OCT, we turn to our main objective, which is to minimize the statistical work fluctuations if a work protocol is applied to a thermal ensemble. For the abovedefined protocol of during which there is no interaction with a heat bath, the inclusive work as a function of is given by Jarzynski.97.PRL (); Jarzynski.07.CRP (); Campisi.11.RMP ():
(8) 
where is simply the phase space coordinates at the end of the protocol. The ensemble average of is
(9) 
In order to suppress work fluctuations across the thermal ensemble, one may now choose an appropriate target function to minimize via OCT. A naive and simple choice of the target function can be just the variance squared of . That is, one may choose the ensemble average as the target function. However, it turns out that, different control fields can yield different values of , and as a result one does not know beforehand. Thus such type of intuitive choice of the control target function is not adopted here. At the point, Jarzynski’s equality comes as a rescue. Note first that
(10) 
where refers to the free energy difference associated with the total Hamiltonian . So long as the initial and final values of are zero (this is easily realized in the next section) , then will be the same as , the free energy difference associated with and . This being the case, no matter what the tobefound control field is, the ensembleaveraged value of is fixed, which is given by . We are thus motivated to design a control target function as follows:
(11) 
The abovedefined form of the control target function is intriguing, because minimization of this function is then to directly suppress the derivations of possible individual values of from their ensembleaverage value predicted by Jarzynski’s equality. In this regard, our OCT framework directly exploits Jarzynski’s equality. It is also expected that the found optimal control field can remove those rare trajectories with rare values of , which slow down the convergence of simulation results towards Jarzynski’s equality.
Iii Optimal Control of Work Fluctuations in Linear Parametric Oscillators
As a benchmark step, in this section we consider parametric (linear) oscillator systems. In particular, for such systems, the control Hamiltonian to realize accelerated adiabatic processes can be found easily, for both scenarios of fastforward adiabatic driving (FFAD) and shortcuts to adiabaticity (STA). We can hence compare the performance of our optimal control fields with those found in FFAD and STA.
The Hamiltonian of a parametric oscillator with a timedependent frequency, all in dimensionless units, is described by
(12) 
For FFAD, the extra control field is Deng.13.PRE (); Jarzynski.13.PRA ()
(13) 
For STA, the control field is found to be Campo.13.PRL ()
(14) 
As an example, we choose a frequency protocol Campo.13.PRL ()
(15)  
such that the above two control fields are indeed zero at or . Besides, note that
(16) 
So during the entire protocol if we choose . In our calculations we choose and in dimensionless units. Then, it is easy to show that
(17) 
That is, along an arbitrary trajectory, the work value is always positive. We further set to be as small as (as compared with ), such that the process will be highly nonadiabatic were there no control fields.
To benchmark our OCT approach outlined above with FFAD and STA, we introduce two kinds of control fields, i.e., and (with to be found numerically), in parallel with the above two control fields and to realize FFAD and STA. The weightage factor in Eq. (4) is set to be small, since here we are not much concerned with the cost of the control field. Nevertheless, the timedependence of can be designed to further alter the profile of the field amplitude. Here we propose to use , where is a small constant. Such timedependence of makes diverge at and , and as a consequence the numerically found will become zero automatically at these two boundary times due to the cost function in the OCT.
With the form of the control field and the weightage factor in OCT both specified, numerical iterations based on Eq. (7) can then yield explicit solutions of . The results are shown in Fig. 1 for both type and type optimal control. Their respective time dependence is also compared with the corresponding results to realize FFAD and STA. Note that in all these cases the control field amplitude is zero in the beginning or at the end. It is also seen that found from our optimal control is different from those in FFAD or STA, though the difference for the shown computational example is not yet dramatic.
Next we attempt to quantitatively characterize the performance in suppressing the work fluctuations, for our OCT approach along with FFAD and STA. In particular, we randomly sample initial phase space points (through a standard importance sampling prcedure) according to the initial thermal probability distribution and then evolve them under the total Hamiltonian . Individual values of are denoted , and the fluctuations in and in are characterized by
(18) 
where the total number of simulation trajectories is chosen to be . The variance in both and in itself under different control schemes are all presented in Table 1. A few interesting observations are in order. First, within expected statistical error due to a finite , the bare system, FFAD, STA, OCT with type field, and OCT with type field all yield the same . Second, though the found time dependence from OCT is different from FFAD or STA, the variances in and in obtained from OCT with type field are all the same as those obtained in FFAD and STA. This confirms that our OCT framework is doing an excellent job in suppressing work fluctuations, to the same degree as accelerated adiabatic processes can reach. This also hints that in the parametric oscillator example here, accelerated adiabatic processes already reach an optimized suppression of work fluctuations. Third, even the mean work from type OCT agrees with those obtained from FFAD and STA. Fourth, the two variances shown in the last row of Table 1, which is for OCT with type control, are however slightly above those in other cases. This is because for type control, as also shown in Fig. 1, the required amplitude of the control field is very high, so a small weightage factor chosen for the cost function can still cause a minor difference. This is also manifested in the last row of Table 1, which is again relatively higher than those obtained in OCT with type field, in FFAD, or in STA.
.
Process  

bare system  
FFAD  
STA  
optimal control of type  
optimal control of of type 
Having successfully benchmarked our OCT approach, we now shed light on the flexibility of our optimal control approach. First of all, we can introduce different time dependence to the weightage factor when accounting for the cost of the control field so as to obtain different , the time dependence of a control field. For example, instead of , we have considered
(19) 
for our type optimal control. The obtained is presented in Fig. 2 in comparison with our previous result. Interestingly, although the timedependence of the control field varies significantly with changes in the cost function weightage factor , all the three cases shown in Fig. 2 yield the same variance of , i.e., . This is a clear demonstration that there are many possible solutions in suppressing work fluctuations to a certain level.
So far we have chosen , which is more or less to simulate an instantaneous limit for the bare system. That is, the frequency of the parameter oscillator is changed rapidly as compared with the system’s own time scale. As shown above, in such a parameter regime our optimal control can suppress work fluctuations to the same degree as that achieved in accelerated adiabatic processes. To further check the usefulness of OCT, we now consider a slower protocol in which and . In this case, during the protocol the system’s own bare Hamiltonian will be important in the time evolution. Our optimal control fields are found to perform also very well in this regime. In particular, without a control field decreases to . This is because nonadiabatic effects and hence work fluctuations are weaker for a slower protocol. Interestingly, with FFAD or optimal control, we still find . One example of found from our type optimal control is presented in Fig. 3 in comparison with the field for FFAD. As seen from Fig. 3, for a slow protocol here the required control fields to suppress fluctuations are much weaker than that presented in Fig. 1
Iv Optimal control of work fluctuations in nonlinear oscillators
We are now ready to apply our OCT approach to systems of nonlinear oscillators under certain protocols. Explicit solutions to realize FFAD and STA cannot be found for general nonlinear systems (excluding scaleinvariant systems), but our OCT framework equally applies. In particular, let us consider the following system Hamiltonian with a timedependent nonlinear term:
(20) 
where is still the timevarying system parameter to model the same protocol defined in Eq. (15). In this system, the nonlinear term vanishes at . This choice is just for computational convenience when we sample the initial states from a thermal ensemble ensemble (which is still Gaussian). Below only type optimal control are presented (note that, as observed earlier, the required field strength for type control is much larger). In addition, for the sake of comparison, we also examine the parallel performance of the previous FFAD field obtained in the absence of the nonlinear term. To stress that the previous FFAD will no longer strictly give rise to accelerated adiabatic processes here due to the nonlinear term in the Hamiltonian, we call such a control approximate FFAD.
As the first computational example, we set and , , , and then numerically compute the fluctuations of . For approximate FFAD, we have ; whereas in our optimal control we have . Since these results are almost identical with our previous result in the absence of a nonlinear term, there must be a reason. The reason is simple. When the temperature is not high enough, the system is largely confined to a small neighborhood of the lowest energy state and hence the nonlinear term only plays a negligible role. This is consistent with a recent experiment Rohringer.13.a1 (), which showed that approximate FFAD can perform well in the presence of some degree of nonlinearity. For the same value of and , we have also considered a slow protocol with and . In this case both optimal control and approximate FFAD still yield .
To enhance the nonlinear effect, we next consider a case with a much higher temperature, i.e., (still with , , and ). The timedependence of the found optimal control field is presented in Fig. 5, which is seen to be very different from that of approximate FFAD. To understand the details better, we also count the number of trajectories that may give a negative work output. As shown previously, for the protocol here the work would be always positive if the nonlinear term were not introduced. So the presence of negative work values in the bare system does indicate the presence of nonlinearity. Detailed computational results are shown in Table 2. It is seen that the performance of optimal control is much better than that achieved by approximate FFAD. While the optimal control field has essentially suppressed almost all negative work values (but one out of one million), approximate FFAD increases probabilities of negative work. The reason why our optimal control approach is so effective in removing negative work values is quite simple. For cases with a positive , the value of with a negative would be too drastically larger than its “target” value , so it will be rejected by our optimization algorithm as extreme or rare values of . Note also that the variance in work under the optimal control field is less than one quarter of that for the bare system, and less than half of that obtained with approximate FFAD. The mean work under the optimal control field is also significantly smaller than that in the bare system or in the FFAD case. Because is an interesting quantity called the dissipated work, a decrease in indicates less dissipated work. We also show in Fig. 5 the work probability distribution . There it is seen that the optimal control field most effectively suppresses very large or very small work values, in addition to an almost complete removal of negative work values.
Process  Probability of negative work  

bare system  
FFAD control  
optimal control of type 
V Optimal Control of Work Fluctuations in Systems with Unknown Parameters
The flexibility in optimal control and early studies of feedback control motivated us to ask whether OCT can be used to suppress work fluctuations in those systems that have unknown system parameters. At the first thought, since with unknown system parameters, it is impossible to calculate and thus our optimal control approach seems not applicable. However, we show in this section that it is possible to refine the control target function on the go, by iteratively guessing the target function through Jarzynski’s equality. This possibility should be good news for the application of Jarzynski’s equality itself, as we can now indeed predict from nonequilibrium work values, with suppressed work fluctuations and hence better performance in reaching the correct based on a finite number of trajectories.
Let us first examine how our OCT framework depends on the knowledge of the bare system Hamiltonian. The first two relations in Eq. (7) require information of the initial thermal distribution and the target function. In other words, the full knowledge of the Hamiltonian in the beginning and at the end of the protocol is needed there. The third and fourth relations in Eq. (7) describe continuous time evolution under the total Hamiltonian including the control field. However, if the protocol is fast, then the main component of the total Hamiltonian can be the control field, so the bare system Hamiltonian may not play an important role there.
To have some quantitative ideas, we consider again the anharmonic Hamiltonian in Eq. (20), i.e., , with and . For this designed case, the quartic term vanishes at and , and hence the target function in our OCT is still fully known. Next we set , but in our construction and implementation of OCT we do not use this piece of knowledge. Instead, we use some wrong values of when computationally search for the OCT field. The knowledge of is used for checking only, after we have obtained the OCT field and start to look into the actual work fluctuations with the control field thus obtained.
As shown in Fig. 6, with used in our optimal control algorithm varying from to , the fluctuations in is not changing significantly. That is, for a wide range of incorrectly assumed values of , the associated optimal control field can still effectively suppress work fluctuations, with ranging from to . Additional numerical investigations of the evolving trajectories further indicate that the total time duration is too short for to play a role, thus confirming our qualitative insights above.
This computational example is enlightening, but we assume there that the (bare) system Hamiltonian is fully known at two boundary times and . Since in such situations can be exactly calculated from the bare system Hamiltonian, we are still one step away from predicting using the work values in a nonequilirium protocol. So our final question is the following, if some parameter in the (bare) system Hamiltonian is indeed unknown to us, how to construct the optimal control target function, suppress the work fluctuations, and eventually predict ?
Borrowing the idea from feedback optimal control theory Judson.92.PRL (), we aim to propose a useful computational feedback procedure as illustrated in Fig. 7. It consists of the following steps:

A certain small number of initial states are first sampled according to system’s thermal distribution, and then evolved in accord with a guessed control field.

Based on the previous step, a rough estimate of is obtained to yield , which is then used to yield/update the control target function.

To optimize the control field based on Eq. (7), one may neglect the effect of the bare system Hamiltonian or preassume some wrong parameter values because the evolution is mainly dictated by the control field during a rapid work protocol.

The control field is then updated by the output from an optimal control algorithm and then all the previous steps are repeated until some convergence threshold is met.
Process  

bare system  
Feedback OCT  
Feedback OCT  
Feedback OCT 
To demonstrate our strategy we consider the following system
(21) 
Note that in this case, the quartic term is not vanishing at , with assumed to be “unknown” when we seek the optimal control field. That is, we will not use this value to construct our control field. Instead, we use some preassumed wrong values of when solving Eq. (7). We still set the quartic term to be zero at for the convenience in initial state sampling. In the first iteration of our numerical experiment, we use a null control field to start with. We update the target function and the optimal control field based on trajectories only. Remarkably, the optimal control field can already converge well after only five iterations of the above four steps. In particular, we compare in Table 3 the results from trajectories, in terms of , , , and . In obtaining the numerical results we have used three different preassumed values of (all are much different from the real value). As seen from the third column of Table 3, for all three cases, the fluctuations in are effectively well suppressed (as compared with the bare case without a control field). The ensembleaverage in all the three cases are also close to the true theoretical value 0.3272 (obtained by numerically computing the partition function of with ). Interestingly, the value of for the bare system case (second column, second row of Table 3) is still slightly away from this theoretical value. This suggests that in the bare system case a highquality convergence towards Jarzynski’s equality has not been achieved with trajectories. Thus, the presence of a control field suppressing work fluctuations is seen to have, albeit slightly, enhanced the convergence of the simulation towards Jarzynski’s equality. To see this more clearly, we show in Fig. 8 how the numerically obtained average value of gradually converges towards the theoretical value , as the number of trajectories increases to . The horizontal line in Fig. 8 represents the true theoretical value. It is seen that all the three cases of feedback control have converged to the true theoretical value with about trajectories, but the bare system case still has a nonnegligible shift from the theoretical result.
Returning to Table 3, from the last two columns it is also seen that the mean work (hence also the dissipated work) as well as the variance of is suppressed by the optimal control field by about one order of magnitude. This significant control over the work output and its fluctuations is achieved even though part of the system parameters is unknown to us! Interestingly, this does not mean that the variance of (see second column of Table 3) will be also significantly reduced by the control field. Qualitatively, note that a higher temperature induces a wider initial probability distribution and hence larger thermal fluctuations in work values, but on the other hand, when calculating the larger work values are still scaled down by the inverse temperature . As such, from the explicit computational example here it is learned that one should not underestimate the implications of a seemingly “small” suppression in the variance of .
Vi Concluding Remarks
In this paper we have proposed an optimal control approach aiming at the suppression of work fluctuations associated with a protocol applied to a thermal ensemble. Our approach is original in the sense that we directly exploits Jarzynski’s equality in defining our control target function. Indeed, the control target function is constructed to bias against the deviation of individual values of from their ensemble average , where is the free energy difference. This approach is shown to be very effective. In the case of a parametric oscillator, the performance of our optimal control approach is simply as good as previous methods based on accelerated adiabatic processes. More importantly, our optimal control approach can equally apply to rather arbitrary nonlinear systems. We also note that a recent study Lu.14.PRA () considered a somewhat related optimal control approach for nonlinear systems, but assuming weak anharmonic potential to validate a perturbative treatment. Our approach can however be applied to systems with strong nonlinear effects, as we have shown through nonlinear oscillator systems with a large quartic term. One might also think that since is needed in the construction of a target function, we may have to know the full information of the system Hamiltonian. Through simple numerical experiments, we showed that this intuition is not true, because a feedback mechanism can help us to refine the control target function on the fly. As such, we claim that the suppression of work fluctuations via a control field is possible even when we do not have full knowledge of the system.
We believe that this study is just a motivator and a starting point along the general issue of controlling and manipulating nonequilibrium thermodynamic properties. Many follow up studies are on the way. For example, it is of immediate interest to extend this study to quantum systems in order to account for quantum effects and quantum fluctuations. So far we have assumed sufficiently rapid work protocols during which systembath interaction is neglected. In the future it is necessary to extend our approach here to open systems. To that direction we note an early study considering classical overdamped motion under a special type of optimization Schmiedl.07.PRL ().
Vii Ackowledgements
We are grateful to Peter Hanggi, Dario Poletti, and Manas Mukherjee for interesting discussions.
Viii Appendix: The Iteration Procedure of OCT
In this section we discuss numerical iteration details in seeking optimal control fields from Eq. (7). The procedure is rather standard in OCT Shi.90.JCP (). First, we divide the relations in Eq. (7) together with Hamilton’s equations of motion into four parts.
(22) 
(23) 
(24) 
(25) 
Next, in numerical calculations, we can consider a sufficiently large but finite phase space area, by considering the initial thermal distribution in Eq. (1). For each initial phase space point sampled, one integrates Hamilton’s equation from to , using a guessed field time dependence . Then the final values of the Lagrange multipliers and can be computed from Eq. (22. With this, Eq. (23) can be integrated backwards. Finally, with , , and all solved numerically, the field time dependence for the next iteration can be obtained by referring to Eq. (24). This iterative procedure is repeated until certain convergence criteria are satisfied. In our implementation of this procedure we simply set the starting to be zero. At the end of each iteration, we update the field time dependence for the next iteration as follows:
(26) 
Basically, the iteration equation means that a difference between final state and target state is calculated after each iteration, and then a modification of the control field is made.
References
 (1) C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
 (2) C. Jarzynski, Phys. Rev. E 56, 5018 (1997).
 (3) G. E. Crooks, Phys. Rev. E60, 2721 (1999).
 (4) M. Campisi, P. Hänggi, and P. Talkner, Rev. Mod. Phys. 83, 771 (2011).
 (5) O. Abah, J. Roßnagel, G. Jacob, S. Deffner, F. SchmidtKaler, K. Singer, and E. Lutz, Phys. Rev. Lett. 109, 203006 (2012).
 (6) C. Bergenfeldt, P. Samuelsson, B. Sothmann, C. Flindt, and M. B¨¹ttiker, Phys. Rev. Lett. 112, 076803 (2014).
 (7) K. Zhang, F. Bariani, and P. Meystre, Phys. Rev. Lett. 112, 150602 (2014).
 (8) J. Deng, Q. Wang, Z. Liu, P. Hänggi, and J. B. Gong, Phys. Rev. E 88, 062122 (2013).
 (9) Y. Zheng and D. Poletti, Phys. Rev. E 90, 012145 (2014).
 (10) A. del Campo, J. Goold, and M. Paternostro, arXiv:1305.3223.
 (11) S. Shi and H. Rabitz, J. Chem. Phys. 92, 364 (1990).
 (12) A. P. Peirce, M. A. Dahleh, and H. Rabitz, Phys. Rev. A 37, 4950 (1988).
 (13) S. Shi and H. Rabitz, Comput. Phys. Commun. 63, 71 (1991).
 (14) K. Sundermann and R. de VivieRiedle, J. Chem. Phys. 110, 1896 (1999).
 (15) Y. Ohtsuki, G. Turinici, and H. Rabitz, J. Chem. Phys. 120, 5509 (2004).
 (16) M. Demirplak and S. A. Rice, J. Chem. Phys. 129, 154111 (2008).
 (17) M. Demirplak and S. A. Rice, J. Phys. Chem. B 109, 6838 (2005).
 (18) M. Demirplak and S. A. Rice, J. Phys. Chem. A 107, 9937 (2003).
 (19) M. V. Berry, J. Phys. A: Math. Theor. 42, 365303 (2009).
 (20) S. Masuda and K. Nakamura, Proc. R. Soc. London Ser. A 466, 1135 (2010).
 (21) X. Chen, I. Lizuain, A. Ruschhaupt, D. GuéryOdelin, and J. G. Muga, Phys. Rev. Lett. 105, 123003 (2010).
 (22) X. Chen, A. Ruschhaupt, S. Schmidt, A. del Campo, D. GuéryOdelin, and J. G. Muga, Phys. Rev. Lett. 104, 063002 (2010).
 (23) S. Choi, R. Onofrio, and B. Sundaram, Phys. Rev. A 84, 051601 (2011).
 (24) J. F. Schaff, P. Capuzzi, G. Labeyrie, and P. Vignolo, New J. Phys. 13, 113017 (2011).
 (25) M. G. Bason, M. Viteau, N. Malossi, P. Huillery, E. Arimondo, D. Ciampini, R. Fazio, V. Giovannetti, R. Mannella, and O. Morsch, Nat. Phys. 8, 147 (2012).
 (26) J. Zhang, J. H. Shim, I. Niemeyer, T. Taniguchi, T. Teraji, H. Abe, S. Onoda, T. Yamamoto, T. Ohshima, J. Isoya, and D. Suter, Phys. Rev. Lett. 110, 240501 (2013).
 (27) C. Jarzynski, Phys. Rev. A 88, 040101 (2013).
 (28) A. del Campo, Phys. Rev. A 84, 031606 (2011).
 (29) S. Deffner, C. Jarzynski, and A. del Campo, Phys. Rev. X 4, 021013 (2014).
 (30) E. Torrontegui, S. Ibanez, S. Martí nezGaraot, M. Modugno, A. del Campo, D. GuéryOdelin, A. Ruschhaupt, X. Chen, and J. G. Muga, Adv. At. Mol. Opt. Phys. 62, 117169 (2013).
 (31) Z. C. Tu, Phys. Rev. E 89, 052148 (2014).
 (32) A. del Campo, Phys. Rev. Lett. 111, 100502 (2013).
 (33) C. Jarzynski, C. R. Phys. 8, 495 (2007).
 (34) W. Rohringer, D. Fischer, F. Steiner, I. E. Mazets, J. Schmiedmayer, and M. Trupke, arXiv 1312.5948 (2013).
 (35) R. S. Judson and H. Rabitz, Phys. Rev. Lett. 68, 1500 (1992).
 (36) X. J. Lu, X. Chen, J. Alonso, and J. G. Muga, Phys. Rev. A 89, 023627 (2014).
 (37) T. Schmiedl and U. Seifert, Phys. Rev. Lett. 98, 108301 (2007).