Exact distribution for work and stochastic efficiency of an isothermal machine

The function of an engine is to convert one form of energy into another form. For example, heat engines converts the energy from the heat reservoir to a useful work, windmills converts the wind energy to the electrical energy as a wind turbine or to pump the water as a windpump, a turbine connected to an electric generator utilizes the energy from the flowing water to convert into an electrical energy, a refrigerator pumps the heat from the cold environment to the hot environment, etc. The performance of an engine depends upon the amount of output power it delivers in the expense of the input power. For instance, a heat engine [1, 2, 3] extracts heat $Q_h$ from the hot reservoir at a temperature $T_h$ and dumps some amount of heat $Q_c$ in the cold reservoir at a temperature $T_c<T_h$ in a cyclic manner, and it generates useful work $W=Q_h-Q_c$. The efficiency $\eta$ of such engine is given by $\eta =W/Q_h$ and it is bounded above by the Carnot efficiency ${\eta }_c=1-T_c/T_h$, i.e, $\eta \le {\eta }_c$ where ${\eta }_c$ is maximum possible efficiency achieved by an engine operating in a quasi-static limit and in a reversible fashion. Hence, an engine operating at Carnot efficiency has power (output work per unit time) zero and is practically useless to do a work in a reasonable amount of time.

In the macroscopic thermodynamics, the system is composed of large number of degrees of freedom. Consequently, fluctuations are negligible, and the observables such as heat, work, entropy change, etc., attain a definite value. In contrast, for a system having small number of degrees of freedom, these observables become stochastic quantities [4, 5, 6]. The probability density function of a stochastic observable contains much wide information than their ensemble average value. In past two and half decades, there have been lot of experimental [7, 8, 9, 10, 11, 12, 13] and theoretical [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28] investigations done to understand the probability density function and their fluctuation relations [29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40] associated with them.

For a microscopic machine or engine converting the input power into the output power, efficiency which is the ratio of the output power and the input power is a stochastic quantity. Recently, a number of studies have been devoted to investigate the probability density function and large deviation function [41] of the stochastic efficiency of a microscopic engine. For example, Varley et al. [42] considered two models of microscopic engine namely Brownian work-to-work converter and photoelectric device to study the large deviation of the stochastic efficiency. In a similar work, Verley et al. [43] obtained a method to compute the large deviation function for stochastic efficiency from the characteristic function of input and output power. Subsequently, they used this method to study the efficiency statistics using a system consisted of two states coupled to two distinct temperature reservoirs. Gingrich et al. [44] studied the efficiency fluctuation and the large deviation for a time-asymmetric stochastic heat engine consisted of a two states where the temperature and the energy levels are varied cyclically in four consecutive steps. The probability density function and the large deviation function for the stochastic efficiency of effusion as a thermal engine were investigated by Proesmans et al. [45]. Polettini et al. [46] derived the full probability density function for the stochastic efficiency when the thermodynamics fluxes obeys the multivariate Gaussian distribution with cumulants proportional to the time of observation. Moreover, it was shown that the probability density function for efficiency has two maxima and one minimum. Proesmans et al. [47] studied an isothermal engine using a Brownian particle driven by two time-periodic external forces where one force serves as load and the other one plays the role of a drive. The statistics of the stochastic efficiency is obtained analytically and also verified with the experiment. In Ref. [48], authors generalized a model of isothermal work-to-work converter engine given in Ref. [42] in the underdamped limit using the stochastic external load and stochastic drive force instead of constant external forces. They obtained the large deviation function and large but finite time probability density function for stochastic efficiency of an isothermal engine. Some more studies in this area one can see in the Refs. [49, 50, 51, 52, 53]. Moreover, several experiments have also been performed to understand the efficiency of a microscopic engine [54, 55, 56].

In the paper, we consider a one dimensional isothermal engine composed of two Brownian particle (say A and B) interacting with each other with harmonic potential where particle B is confined in a harmonic trap. When time $t\le 0$, the minimum of the harmonic trap is stationary and is at the origin of $x$-axis. Hence, the system has equilibrium Boltzmann’s distribution at $t=0$. At $t>0$, the minimum of the harmonic confinement is dragged with a constant velocity, and a load is attached to the Brownian particle A. We study the work done on particle A ($W_A$) and particle B ($W_B$), and the stochastic efficiency $\eta$ of the isothermal engine which is the ratio of the work done against the load on particle A to the work done on particle B by dragging the minimum of the harmonic potential with constant velocity. The probability density function for $W_A$, $W_B$ and $\eta$ are obtained for all time.

The paper is organized as follows. In section 2, we give the model system and definitions of work done on particle A, particle B, and the stochastic efficiency. Section 3 contains the joint distribution for the work done on particle A and particle B. The transient fluctuation theorem [30] is studied for work done on particle A and B in section 4. In section 5, we discuss the transient fluctuation theorem for total work done on the system. The probability density function for the stochastic efficiency is computed and phase diagrams are shown which suggests the sign of efficiency where the peak of the density function occurs in section 6. We summarize our paper in section 7. Some of the results are given in Appendix.

Consider a model of a machine consists of two Brownian particles (say particle A and B) coupled by a harmonic spring of stiffness $k$. Suppose particle B is trapped in a harmonic confinement of stiffness $k_0$. The whole setup is immersed in a heat bath of a constant temperature $T$. The potential energy of the system is given by

where $z$ is the minimum of the harmonic trap which is time dependent, and $x_A$ and $x_B$ are the positions of particle A and B, respectively.

Suppose the minimum of the harmonic trap is moved with a constant finite velocity $v$ and is at $z=vt$ at time $t>0$, and a constant finite load $F$ is attached to particle A from time $t>0$ to $t=\tau$. The schematic diagram of the machine is shown in figure 1. Note that similar model for a Brownian particle confined in a harmonic potential ($k=0$) whose minimum is dragged with a given velocity is already studied both theoretically [19, 20, 21, 22] and experimentally [7, 8].

In this paper, our aim is to compute the efficiency of a machine which does work against the load $F$ attached to particle A by dragging the harmonic trap which confines the particle B from time $t>0$ to $t=\tau$. When time $t\le 0$, there was no load attached to particle A, and the minimum of the trap was at the origin, i.e., $z=0$. Thus, the system obeys the steady state distribution at time $t=0$:

where the row vector $U_0^T=\left[x_A\right(0),x_B(0\left)\right]$, and the correlation matrix

where $\delta =k/k_0$ is the dimensionless coupling parameter.

For $t>0$, the dynamics of the system given in (2) and (3), can be rewritten in a matrix form as

where column vectors $\zeta \left(t\right)=\left[{\xi }_A\right(t),{\xi }_B(t){]}^T$, $B\left(t\right)=(F/\gamma ,z/{\tau }_{\gamma }{)}^T$, and the matrix $A$ is

The solution of above equation at time $0<t\le \tau$ is given by

where the symmetric matrix $G\left(t\right)=e^{-\left(t/{\tau }_{\gamma }\right)A}$ is given in A.

Since the correlation between $W_A$ and $W_B$ is zero: $C_{AB}=0$, the joint distribution $P_{\tau }(W_A,W_B)$ can be written in factorize form: $P_{\tau }(W_A,W_B)=P_A\left(W_A\right)P_B\left(W_B\right)$, where

For convenience, we have dropped the subscript $\tau$ in $P_r\left(W_r\right)$.

In figure 2, we have shown the comparison of the analytical result of the probability density function given in (23), with the numerical simulation. The plots show that there is nice agreement between theory and numerical simulation.

For a system initially in the equilibrium and then driven away from the equilibrium using external driving, the stochastic quantity $\Omega$ is observed. Let $P\left(\Omega \right)$ be the probability density function of $\Omega$. When $\Omega$ satisfies the transient fluctuation theorem (TFT) [30], it obeys the following relation

The above relation states that the probability of getting positive values of $\Omega$ is exponentially favourable than that of negative values.

In our case, we analyze TFT for both $W_A$ and $W_B$. It is clear that TFT may not hold for all parameters (see C). However, we investigate TFT in the large and small time limit.

Using above equations, one can see that ${\tilde{C}}_{AA}=2{\tilde{\mu }}_A$ and ${\tilde{C}}_{BB}=2{\tilde{\mu }}_B$ in the limit $\delta \to 0$. Therefore, TFT for both work done ($W_A$ and $W_B$) is satisfied in the weak coupling limit and for small observation time ($u_{\tau }\ll {\delta }^{-1}$).

In figures 3 (a) and  3(c), the comparison between the analytical probability density function for work done on particle A and B given in (23) in which the cumulants $({\mu }_A,{\mu }_B,C_{AA},C_{BB})\to ({\tilde{\mu }}_A,{\tilde{\mu }}_B,{\tilde{C}}_{AA},{\tilde{C}}_{BB})$ are given in (25)–(28), with the numerical simulation is shown for coupling strength $\delta =0.001$ and observation time $\tau =10\ll {\delta }^{-1}$. TFT for work done on particle A and B is shown in figures 3(b) and  3(d) for respective parameters. In all plots, blue points are the obtained from numerical simulations and the red dashed lines are the analytical results. In figures 3(b) and  3(d), red dashed lines have slope unity. Therefore, in the weak coupling limit $(\delta \ll 1)$ and observation time $u_{\tau }\ll {\delta }^{-1}$, TFT is satisfied.

In figures 4(a) and  4(c), we have shown a comparison of analytical probability density function for work done on particle A and particle B, respectively, given in (23) in which cumulants $({\mu }_A,{\mu }_B,C_{AA},C_{BB})\to ({\overline{\mu }}_A,{\overline{\mu }}_B,{\overline{C}}_{AA},{\overline{C}}_{BB})$ are given as (29)–(32), with the numerical simulation for coupling strength $\delta =0.1$ and time $\tau =50\gg {\delta }^{-1}$. The variations of function $\ln \left[P_i\right(W_i\left)/P_i\right(-W_i\left)\right]$ are shown against $W_i$ in figures 4(b) and  4(d) where magenta solid line in figure 4(b) has slope unity. Thus, in the weak coupling limit $(\delta \to 0)$, ${\lambda }_2\sim \delta$, the violation of TFT can be seen for both work done ($W_A$ and $W_B$) for large observation time limit $u_{\tau }\gg {\delta }^{-1}$.

In above two cases, we have studied the TFT for $W_A$ and $W_B$ in the weak coupling limit ($\delta \to 0$) and showed that TFT would hold in the limit $u_{\tau }\ll {\delta }^{-1}$ whereas violation can be seen in the large observation time $u_{\tau }\gg {\delta }^{-1}$. This is because in the small time limit, the effect from the other particle will not affect the TFT of the work done on the observed particle. However, when the time of observation is large $(u_{\tau }\gg {\delta }^{-1})$, the effect from the other particle may appear as the term coupling ($\delta$) times relative separation $(x_A-x_B)$ becomes relevant [see (2) and (3)] which leads to the violation of TFT of the work done of the observed particles even in the weak coupling limit. Similar results one can see in Refs. [26, 28].

In the above section, we have analyzed TFT for work done on particle A and B. In this section, we study the TFT for the total work done on both particles.

Therefore, the probability density function for total work done $W$ is given by

Using the above relation, one can see that $P\left(W\right)/P(-W)=e^W$, i.e., TFT for total work done on both particles is satisfied. Therefore, it is clear that when degrees of freedom (DOFs) having same time of relaxation are coupled and driven out of equilibrium, total work done on all DOFs obeys the TFT for all parameters.

Figure 5(a) and  5(b) show the comparison of the analytical results for the probability density function for $W$ given in (37) and the function $\ln \left[P\right(W\left)/P\right(-W\left)\right]$ with the numerical simulation, and they have nice agreement. Figure 5(b) indicates that the total work done $W$ obeys the transient fluctuation theorem.

It can be seen that the function $I(\eta ,\tau )$ has two extrema, i.e., at $\overline{\eta }$ and ${\eta }^{\ast }$ where

Moreover, $I(\eta ,\tau ){|}_{\eta =\overline{\eta }}=0$. The function $I(\eta ,\tau )$ has a maximum and a minimum at $\eta =\overline{\eta }$ and $\eta ={\eta }^{\ast }$, respectively. Note that both $\overline{\eta }$ and ${\eta }^{\ast }$ do not depend upon the parameter $\alpha$ (relative strength acting on particle A, i.e., $\alpha =F/\sqrt{\gamma T}$). The efficiency $\overline{\eta }$ can have any sign depending upon the parameters $\tau$, $\delta$, and $\theta$. To understand the nature of $\overline{\eta }$, we have plotted phase diagram in $(\theta ,\delta )$ plane as shown in figure 6 for fixed time $\tau$. In figure 6, the light red shaded regions correspond to the areas where the efficiency $\overline{\eta }$ is positive. Note that $C_{AA}$ and $C_{BB}$ are positive [see (43)]. Therefore, $\overline{\eta }$ and ${\eta }^{\ast }$ have opposite sign. Hence, the unshaded regions in figure 6 represent the areas where ${\eta }^{\ast }$ is positive otherwise negative.

Figure 7 shows the comparison of analytical results of the probability density function $p_{\tau }\left(\eta \right)$ given in (38) with the numerical simulation results for both signs of $\theta$. While in figure 7(a) the comparison is shown for parameters $\theta =-1.0$, $\alpha =0.5$, and coupling parameter $\delta =0.3$ at $\tau =10$, $\tau =30$ and $\tau =50$, the parameters $\theta =0.5$, $\alpha =0.5$, and coupling parameter $\delta =0.1$ at $\tau =30$ are taken in figure 7(b). In figure 7(a), the red ($\tau =10$), orange ($\tau$=30) and magenta ($\tau$=50) dashed lines are the analytical result given by (38) whereas the blue dots are obtained from the numerical simulations at corresponding times $\tau$. Similarly, the red dashed lines correspond to the analytical result given by (38) and blue dots are obtained from numerical simulation in figure 7(b). Both of figures show that there is a nice agreement between theory and numerical simulation. From figure 7(a), it is clear that the peak of the density function $p_{\tau }\left(\eta \right)$ shifts from negative to positive side as time $\tau$ is increased from $\tau =10$ to $\tau =30$ as shown in figures 6(a)— 6(c).

We have also compared the large time ($u_{\tau }\gg {\lambda }_2^{-1}$) analytical results for the probability density function $p_{\tau }\left(\eta \right)$ given in (38) in which $({\mu }_A,{\mu }_B,C_{AA},C_{BB})\to ({\overline{\mu }}_A,{\overline{\mu }}_B,{\overline{C}}_{AA},{\overline{C}}_{BB})$ as given in (29)–(32), with the numerical simulations in figure 8(b). These comparisons indicate that there is a nice agreement between theory and numerical simulations.