Correlation-enhanced Stability of Microscopic Cyclic Heat Engines

For cyclic heat engines operating in a finite cycle period, thermodynamic quantities have intercycle and intracycle correlations. By tuning the driving protocol appropriately, we can get the negative intercycle correlation to reduce the fluctuation of work through multiple cycles, which leads to the enhanced stability compared to the single-cycle operation. Taking the Otto engine with an overdamped Brownian particle as a working substance, we identify a scenario to get such enhanced stability by the intercycle correlation. Furthermore, we demonstrate that the enhancement can be readily realized in the current experiments for a wide range of protocols. By tuning the parameters within the experimentally achievable range, the uncertainty of work can be reduced to below $\sim 50 \%$.

For cyclic heat engines operating in a finite cycle period, thermodynamic quantities have intercycle and intracycle correlations. By tuning the driving protocol appropriately, we can get the negative intercycle correlation to reduce the fluctuation of work through multiple cycles, which leads to the enhanced stability compared to the single-cycle operation. Taking the Otto engine with an overdamped Brownian particle as a working substance, we identify a scenario to get such enhanced stability by the intercycle correlation. Furthermore, we demonstrate that the enhancement can be readily realized in the current experiments for a wide range of protocols. By tuning the parameters within the experimentally achievable range, the uncertainty of work can be reduced to below ∼ 50%.
Nevertheless, many studies of cyclic heat engines so far consider single-cycle operation and focus on the performance within a single cycle. In these studies, fluctuations of the thermodynamic quantities usually include only the intracycle correlation. In the quasistatic limit, since the thermal noise erases the correlation among thermodynamic quantities in different cycles [13,38], it is sufficient to describe the fluctuations focusing on a single cycle. However, to get the nonzero power output, we need to operate engines in a finite cycle period. In this case, the effect of the intercycle correlation becomes non-negligible. Therefore, for engine operations over multiple cycles, assessing the performance within a single cycle is insufficient. Instead, assessments of the engine performance should address the global process over multiple cycles to include intercycle correlations.
However, the role of the time correlation in fluctuations of thermodynamic quantities has not been thoroughly explored. Since engines are supposed to operate over multiple cycles consecutively with a finite cycle period in practical situations, there is a great demand for a scheme to prevent the degradation of performance in multiple cycles by the intercycle correlation effect. In this Letter, by clarifying the effect of time correlation of work in microscopic heat engines with a finite cycle period, we identify such a scheme to reduce the fluctuation of work output. Since the fluctuation of work output is comparable to or even bigger than the average of work output in current experiments of small heat engines [4,5], reducing the fluctuation of work output is a crucial issue. Taking an example of the Otto engine using a Brownian particle as a working substance, we demonstrate that the reduction of the fluctuation of work output can be realized in a robust manner in the current experiments, and this reduction can be more than 50%.

Setup.
We study a small cyclic heat engine whose working substance is in contact with a heat bath with the controllable temperature T (t) (we set the Boltzmann constant k B = 1 throughout the Letter). The working substance is described by the Hamiltonian H(Γ,t) with an external control parameter λ (t), where Γ is the microstate of the working substance in the phase space. The engine is driven by time-periodically modulating T and λ with period τ, i.e., T (t) = T (t + τ) and λ (t) = λ (t + τ). Under such a protocol, we assume the engine is already driven into a periodic state with the probability distribution function (PDF) satisfying p(Γ,t) = p(Γ,t + τ) after running many cycles [20]. Therefore, we can represent time t by the phase as θ = 2πt/τ, and assign the initial phase θ 0 for the staring point of the cycle.
The work W (n) θ 0 extracted through n cycles with the initial Relation between the single-cycle and multicycle uncertainties. To discuss the relationship between the uncertainties within a single cycle ∆ (1) θ 0 and infinite cycles ∆ ∞ , we consider an overdamped Brownian particle trapped in a onedimensional harmonic oscillator potential with the Hamiltonian Here, λ (t) is the stiffness of the potential which serves as a mechanical control parameter and x(t) is the position of the Brownian particle. This system is described by the Ornstein-Uhlenbeck process [65]. The correlation function is derived from the solution of the Itô stochastic differential equation for this process [66]. The resulting correlation function φ (t,t ′ ) for t < t ′ is given by where µ is the mobility. In addition, since φ (t,t) is periodic in time, we have The covariance function of power becomes [66]. From Eqs. (6) and (7), we get the following properties of the covariance function: Therefore, C(t,t ′ ) decays exponentially in time when |t − t ′ | ≫ τ, and the correlation time of work is given by τ corr = 2µ τ 0 dtλ (t)/τ. From the above properties of C(t,t ′ ), we can write the intercycle correlation C (2) θ 0 within the two successive cycles as C (2) In the same way, we can write C (n) θ 0 in terms of a, γ(θ 0 ), and Var W (1) θ 0 [66]. Then, the uncertainty for n cycles reads [66] n∆ (n) where s n ≡ (1 − a n )/(1 − a) ≥ 1. For infinite cycles, we get For finite n cycles, the uncertainty derived from Eqs. (9) and (10) reads If the intercycle correlation C (10) and (11) [67]. This means that the negative covariance of work between two successive cycles, C θ 0 < 0, indicates the reduction of uncertainty of work in multiple cycles. It is vice versa for the positive intercycle correlation. It is noted that the essential point of the above discussion is the exponential decay in time of the correlation functions. Even if the effect of inertia is non-negligible beyond the overdamped limit, the correlation functions can still be exponential in time with a smaller correlation time in the overdamped regime [66]. In addition, in the strongly underdamped regime, the correlation functions can be well approximated by an exponentially decaying function with a large correlation time τ corr ≃ γ −1 obtained by averaging over the rapid oscillation [66]. Therefore, for the both cases, the above results can still apply, but with a different value of τ corr .
It is possible to observe the enhanced stability due to the negative intercycle correlation when τ τ corr . To show this effect, below we consider a simple Brownian Otto engine, where the analytical result can be obtained. However, a similar result is also obtained for the Carnot cycle [66].
Results for the Brownian Otto cycle. Next, taking the Brownian Otto engine as an example, we demonstrate that the negative intercycle correlation can be realized in a wide range of parameters in the driving protocol. We still consider an overdamped Brownian particle in a harmonic oscillator potential described by the Ornstein-Uhlenbeck process. In this model, since the PDF p(x,t) of any periodic state is Gaussian, we can define the effective temperature T eff of the Brownian particle given by T eff = λ x 2 [18]. The Brownian Otto engine consists of two isochoric and two isentropic strokes as shown in Fig. 1 [68]. During the hot (cold) isochoric strokes, the temperature T of the bath and the parameter λ are fixed at T h and λ h (T c  com . Insets of (b) show cycle diagrams on the λ -T eff plane for each typical point. Here, we set T h = 1, λ h = 1, and λ c = 0.5 for all the three cycle diagrams. and λ c ), respectively, for the duration τ h (τ c ) with λ c < λ h . During the isentropic strokes, T and λ are quenched simultaneously in a way such that the Shannon entropy S ≡ − ln p is unchanged [18]. We assume that the isentropic strokes are instantaneous, so that the cycle period is given by τ = τ h + τ c .
For each mth cycle, we assign an odd integer i = 2m − 1 for the isentropic expansion stroke and an even integer i = 2m for the isentropic compression stroke (see the strokes labeled "1" and "2" in Fig. 1 for m = 1). Since work is done only in the isentropic strokes, the fluctuation of work can take two values depending on whether θ 0 is in the hot or cold isochoric strokes. Therefore, the analysis can be divided into two cases according to the initial phase: the cycle starts before the isentropic expansion or compression. Then, we get the variance of work for the two cases, Var[W (1) exp ] = ∑ i, j=1,2 C i j and Var[W (1) com ] = ∑ i, j=2,3 C i j , respectively. Here, the subscript θ 0 in W (1) θ 0 is replaced by "exp" and "com" for clarity, and In this example, the correlation function φ i j ≡ φ (t i ,t j ) is analytically solvable [66].
From the analytical solution of φ i j , one can find that the uncertainties ∆ com , and ∆ ∞ depend on three parameters [66]: Here, e h and e c are measures of the incompleteness of the equilibration in the hot and cold isochoric strokes, respectively, and φ r describes the spread of the width of the PDF of the Brownian particle during the hot isochoric strokes. Since we are interested in the heat engine, the mean value of work should be positive, W (1) = (λ h − λ c )(φ 11 − φ 22 )/2 > 0, and thus φ r > 1. In addition to the condition φ r > 1, the region of φ r is upper bounded as φ r < 1/e c because the parameters e h , e c , and φ r are constrained by [66] ( where R ≡ T c λ h /(T h λ c ) describes the reversibility with Since 0 < e c < 1, 0 < e h < 1 < φ r , and 0 < R < 1, we get φ r < 1/e c from Eq. (12).
com , and ∆ ∞ is smaller than the others. Regions II and III are of our interest, where the uncertainty ∆ ∞ is smaller than those for a single cycle irrespective of the starting point of the cycle. Figure 2(a) tells that, if the equilibration in the cold isochoric strokes is sufficient with e c < 1/φ 2 r , we can get the reduction of the uncertainty for multiple cycles. It is noted that, to obtain this reduction, only the degree of equilibration in the cold isochoric strokes matters, but not that in the hot isochoric strokes.
We can provide a physical understanding of Fig. 2(a). An example of the protocol λ (t) of the Brownian Otto engine starting before the isentropic expansion (stroke 1) is shown in  θ 0 satisfy C i,i+1 = −e c C ii for odd i and C i,i+1 = −e h C ii for even i. The correlation decays with n as C i, j+2n = a n C i j , where a = e c e h < 1. Therefore, C (2) exp is proportional to C 11 a −C 22 e h ∝ φ 2 r − 1/e c . As we have discussed, the necessary and sufficient condition for the reduction of uncertainty is C (2) θ 0 < 0, which gives e c < 1/φ 2 r corresponding to regions II and III. In the same way, for cycles starting before the isentropic compression, the intercycle correlation C given by C r /e h , but it is always smaller than zero. Therefore, we have ∆ ∞ < ∆ (1) com for arbitrary e h . Summarizing the results for the above two cases, we get ∆ ∞ < ∆ com which depends on the intracycle correlation. Finally, we discuss the role of the temperature of the bath and the experimental feasibility to get the reduction of the fluctuation by the intercycle correlation. First, we consider the mean value of the power P. As obtained in Ref.
[71], P depends on six parameters: T h , T c , λ h , λ c , τ h , and τ c [66]. At any point in the region of 0 < e c < 1 and 0 < e h < 1, the power can be set to any positive value for a given φ r by tuning the remaining free parameters, such as T h , T c , and λ h . Figures 4(b)-4(e) show the power for different values of T h /T c . It can be seen that the point in the e c -e h plane giving the maximum power can be located in region I or II by tuning T h /T c . It is noted that we have R > T c /T h for the Otto engine with P > 0 (η > 0) from Eq. (13). Second, we discuss the role of the temperature ratio T h /T c in the correlation-enhanced stability. Figure 4(a) shows a contour plot of R as a function of e c and e h for a fixed value of φ r . As can be seen from Fig. 4(a), if R is larger than that at e c = 1/φ 2 r and e h = 0, it is guaranteed that we are in either region II or III. From Eq. (12), we find that this condition is R > 1/(φ r + 1), or A sufficient condition to satisfy this inequality is T h /T c < 2, which is easy to realize in experiments. In experiments of microscopic heat engines with Brownian particles [4-6, 8, 9], one of the heat bath temperatures (commonly T c ) is usually set to be the room temperature: T c ∼ 300 K. In such a case, if T h is 300 K < T h < 600 K which is indeed the case in typical experiments [4,5], it is guaranteed that the fluctuation of work in the Brownian Otto cycle is always reduced for multiple cycles irrespective of the other parameters. To demonstrate the large reduction of ∆ ∞ by the intercycle correlation, we plot ∆ ∞ /∆ (1) exp and ∆ ∞ /∆ (1) com as functions of T h /T c in Fig. 5 for parameter values accessible in current experiments. Since the work output is zero at T h /T c = λ h /λ c and increases with T h /T c , the region of T h /T c shown in Fig. 5 gives positive work output. It is noted that, compared to the above-mentioned sufficient condition, T h /T c < 2, for ∆ ∞ < ∆ com , we can obtain this reduction of ∆ ∞ in a much wider temperature region of T h /T c 5.4. Furthermore, the reduction of ∆ ∞ over the single-cycle uncertainties can be very large by appropriately tuning the parameters and protocol. At T h /T c ≃ 2.3 where the red and blue lines cross, we have the same reduction rate for an arbitrary starting point. In this case, the uncertainty ∆ ∞ can be reduced to less than 60% of the single-cycle uncertainties. If we set the starting point of the cycle before the isentropic compression stroke (i.e., the case of the red line), ∆ ∞ can be reduced to even below 50% of the single-cycle uncertainty.
Conclusion. Our work has clarified the consequences of time correlation of work over different cycles in cyclic heat engines. If the cycle period is finite, focusing on one cycle is insufficient to discuss fluctuations of the performance of the microscopic heat engines. In particular, taking advantage of the intercycle correlation, the stability of the work output for the multicycle operation can be improved over the singlecycle one. Since such an improvement can be realized in a wide range of protocols, one can further optimize the other performance of the engine (such as efficiency, power, and uncertainty within each cycle). Furthermore, we have demonstrated that our findings can be readily realized in the current experiments. By tuning the parameters within the experimentally achievable range, the uncertainty of work output for infinite cycles can be reduced to less than 50% of the uncertainty for each single cycle. Since the fluctuation of work output can be even larger than the average of the work output in the current experiments [4,5], our result should provide an important step toward the realization of microscopic heat engines for practical use. The effect of time correlation in other kinds of heat engines, such as autonomous heat engines and selfoscillating heat engines [24], is an interesting future problem.  θ 0 and C ∞ θ 0 are the same, so that the conditions C (2) θ 0 < 0 and C ∞ θ 0 < 0 for ∆ ∞ < ∆ (1) θ 0 are equivalent.
[68] In experiments of the Brownian heat engines, adiabatic strokes are often replaced by isentropic strokes [5,8,72] since the working substance is always in contact with the environment (water), so that it is impossible to thermally isolate from the environment. In these isentropic strokes, the parameter λ and the temperature are controlled to keep the mean value of the entropy of the  −2 and the radius of the Brownian particle r = 0.5 µm given by Ref.
[69], the mobility µ is obtained by µ = 1/(6πηr). Motion of a Brownian particle in a one-dimensional (1D) harmonic oscillator potential is a Gaussian process x(t) described by the following Itô stochastic differential equation [S1]: where dW is the Wiener noise. The solution of Eq. (S1) is Thus the mean value x(t) is given by Because of the periodicity of the phase-space distribution function of the Brownian particle with the cycle period τ, we have The correlation function φ (t,t ′ ) is given by From the periodicity, x 2 (τ) = x 2 (0) , and Eq. (S2), we have Therefore, the correlation function φ (t,t ′ ) reads

The effect of inertia
As can be seen from Eq. (S7), the correlation function decays exponentially in the overdamped regime. Our results, Eqs. (10) and (11) of the main paper, are obtained from such exponentially decaying correlation function. Here, we show that, even if the inertia is non-negligible, the correlation functions can still be exponential in time in the strongly overdamped regime with γ ≫ ω, where ω = λ /m is the frequency of the harmonic oscillator trapping potential and mγ = µ −1 is the frictional coefficient. In addition, in the strongly underdamped regime with γ ≪ ω, the correlation function can be well approximated by an exponentially decaying function in time obtained by averaging over the rapid oscillation. Therefore, for the both cases, the similar argument in the main paper can still apply with the effect of inertia, but τ corr becomes smaller (larger) in the strongly overdamped (underdamped) regime compared to that in the overdamped limit. As a result, it is harder (easier) to observe the correlation-enhanced stability in the strongly overdamped (underdamped) regime because of the necessary condition: τ τ corr .
We assume that λ changes in a timescale λ /λ much larger than 2π/ω, and the correlation function oscillates with period 2π/ω much smaller than τ in the strongly underdamped regime with γ ≪ ω [S3]. Averaging over the coarse-grained timescale, which is sufficiently larger than 2π/ω but sufficiently smaller than τ, the correlation function of work is given by where "--" means the average over the coarse-grained timescale. In the timescale much smaller than λ /λ , where the change of λ ≡ mω 2 is negligible, we have the following expression from Eq. (S10): Averaging over each oscillation period, we obtain Therefore, also in the strongly underdamped regime of γ ≪ ω, the correlation function of work is well approximated by an exponentially decaying function in time with the large correlation time τ corr = γ −1 .
in Fig. S1. In accordance with the periodicity of C(t 1 , t 2 ), C(t 1 + τ,t 2 + τ) = C(t 1 ,t 2 ), equivalent subdomains are denoted by D i with the same integer i. Now we introduce and we have and For n cycles, the variance of work is given by Here, we have used the following summation formula of the series: with s n ≡ (1 − a n )/(1 − a). From Eq. (S23), the uncertainty for n cycles reads n∆ (n) From Eq. (S23), we can readily identify C (n) C. C: Covariance function of power for an overdamped Brownian particle trapped in a one-dimensional harmonic oscillator potential For a Gaussian process, x(t), all the higher-order correlation functions can be decomposed into products of two-point correlation functions using Wick's theorem [S7]. For example, we have Therefore, for an overdamped Brownian particle trapped in a one-dimensional harmonic oscillator potential with the Hamiltonian, H(x,t) = 1 2 λ (t) x(t) 2 , the covariance function of power is given by com , and ∆ ∞ depend on three parameters: e h , e c , and φ r Here, we consider the Brownian Otto engine with the protocol defined in the main paper, and show that the uncertainties ∆ com , and ∆ ∞ depend on e h , e c , and φ r . First, let us consider the case starting before an isentropic expansion stroke. The mean value W (1) of work is given by , where λ h and λ c are the stiffness of the potential during the hot and cold isochoric stroke, respectively. The variance Var W exp of work is given by 2 11 with e c ≡ exp(−2µλ c τ c ) and τ c being the duration of the cold isochoric stroke, the uncertainty ∆ Similarly, in the case of starting before an isentropic compression stroke, the variance Var W com of work is given by Finally, we discuss ∆ ∞ . Since ∆ ∞ is given by θ 0 , we focus on γ(θ 0 ), which is defined as Eq. (S21). For clarity, we write γ(θ 0 ) in the case of starting before an isentropic expansion (compression) stroke as γ exp (γ com ). Here, we consider γ exp as an example. Figure S2 shows the covariance matrix C i j ≡ C(t i , t j ), which appears in the integrand of the expression of γ(θ 0 ) [Eq. (S21)]. For γ exp , the elements of C i j enclosed by the dashed lines contribute, so that γ exp can be written as Similarly, γ com can be written as Thus the uncertainty ∆ ∞ for infinite cycles is given by Therefore, from Eqs. (S31), (S33), and (S36), the uncertainties ∆ com , and ∆ ∞ depend on the three parameters: e h , e c , and φ r . In addition, for n cycles, we obtain E. E: Constraint on e h , e c , and φ r : derivation of Eq. (12) We consider the protocol starting from the adiabatic compression with λ (t) = λ h for 0 < t < τ h , λ (t) = λ c for τ h < t < τ, T (t) = T h for 0 < t < τ h , and T (t) = T c for τ h < t < τ, where λ (t) = λ (t + τ) and T (t) = T (t + τ). From Eq. (S6), the correlation function φ 22 ≡ φ (τ, τ) = φ (0, 0) is given by Here, we have used the following results: In the same way, the correlation function φ 11 ≡ φ (τ h , τ h ) is given by The ratio φ r ≡ φ 11 /φ 22 reads which gives a constraint on e h , e c , and φ r for the Otto engine: F. F: The role of e h , e c , and φ r on the power In the main paper, we have discussed about the correlation-enhanced stability in the Brownian Otto engine, and we have found that the uncertainty of work depends only on the three parameters: e h , e c , and φ r . Here, we show that the mean value of power for given e h , e c , and φ r can reach arbitrary value for 0 < e c < 1 and 0 < e h < 1. In addition, we give a detailed calculation of the power plotted in Figs. 4 Power of the Brownian Otto engine is given by [S5] which depends on six parameters: T h , T c , λ h , λ c , τ h , and τ c . At any point in Fig. 2(a), only three parameters are given: e h = exp(−2µλ h τ h ), e c = exp(−2µλ c τ c ), and φ r . Here, φ r can be substituted by R = T c λ h /(T h λ c ) from Eq. (12). Therefore, we can still choose arbitrary positive values for the remaining three of the six parameters. For example, when λ h τ h , λ c τ c , and R are given, we can still choose arbitrary positive values for T h , T c , and λ h with T c /T h < R. Here, the constraint on T c /T h is from the fact that η = 1 − λ c /λ h = 1 − (T c /T h )/R > 0 for the Otto engine. In this case, we can rewrite Eq. (S43) as follows: We discuss the effect of intercycle correlation in the Carnot cycle with the model by Schmiedl and Seifert [S6]. We calculate the uncertainty of work output at maximum power for the Carnot cycle, and show that the negative intercycle correlations can be obtained in some regions of the parameter space for the Carnot cycle as well.
In this model, the working substance is an overdamped Brownian particle trapped in a harmonic oscillator potential. According to Ref. [S6], the cycle consists of two adiabatic jumps and two isothermal strokes. The durations of the adiabatic strokes are negligible and those of the hot and cold isothermal strokes are τ h and τ c , respectively. During the hot (cold) isothermal stroke, 0 < t < τ h (τ h < t < τ = τ h + τ c ), the water temperature T h (T c ) is constant. The heat engine protocol λ (t) is optimized to yield the maximum work output with given boundary values of the variance of the particle position x: σ min ≡ σ (0) = σ (τ) and σ max = σ (τ h ), where σ ≡ x 2 − x 2 . Since the power is nonzero at each moment,Ẇ = 0, the correlation of work for the Carnot cycle is more complicated than that for the Otto cycle considered in the main paper. From our numerical result, we find that the uncertainty of work at maximum power depends only on the ratio σ r ≡ σ min /σ max and η C ≡ 1 − T c /T h . Figure S3 shows that there are some regions (regions II and III) in the parameter space where we have the correlation-enhanced stability.  [S6]. ∆ ∞ is compared to ∆ (1) θ 0 with the starting point t 0 = 0 + , τ − h , τ + h , or τ − corresponding to the node in the Carnot cycle. In regions II and III, ∆ ∞ is smaller than ∆ (1) θ 0 for any of the four choices of t 0 .