Energetics and performance of a microscopic heat engine based on exact calculations of work and heat distributions
We investigate a microscopic motor based on an externally controlled two-level system. One cycle of the motor operation consists of two strokes. Within each stroke, the two-level system is in contact with a given thermal bath and its energy levels are driven with a constant rate. The time evolution of the occupation probabilities of the two states are controlled by one rate equation and represent the system’s response with respect to the external driving. We give the exact solution of the rate equation for the limit cycle and discuss the emerging thermodynamics: the work done on the environment, the heat exchanged with the baths, the entropy production, the motor’s efficiency, and the power output. Furthermore we introduce an augmented stochastic process which reflects, at a given time, both the occupation probabilities for the two states and the time spent in the individual states during the previous evolution. The exact calculation of the evolution operator for the augmented process allows us to discuss in detail the probability density for the performed work during the limit cycle. In the strongly irreversible regime, the density exhibits important qualitative differences with respect to the more common Gaussian shape in the regime of weak irreversibility.
pacs:02.50.Ey, 05.40.-a, 05.70.Ln
Non-equilibrium phenomena in the presence of time-varying external fields are of vital interest in many areas of current research [1, 2, 3]. Examples are aging and rejuvenation effects in the rheology of soft-matter systems and in the dynamics of spin glasses, relaxation and transport processes in biological systems such as molecular motors, ion diffusion through membranes, or stretching of DNA molecules, driven diffusion systems with time-dependent bias, and nano-engines. With minimization of the system size thermal fluctuations become increasingly relevant. In these systems it is useful to introduce microscopic heat and work quantities as random variables whose averages yield the common thermodynamic quantities. Averages over functions of these microscopic heat and work quantities yield generalized fluctuation theorems [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. In this context mesoscopic engines operating between different heat baths under non-equilibrium conditions have received increasing attention. The variety of models can be roughly classified according to the dynamical laws involved. In the case of the classical stochastic heat engines, the state space can either be discrete or continuous (c.f., for example, [15, 16, 17] and the references therein). Examples of the quantum heat engines are studied, e.g., in [18, 19].
The traditional consideration of efficiency of heat engines operating between two baths at temperatures and leads to the Carnot upper bound . The bound is only achieved under reversible conditions where the state changes require infinite time and hence the power output is zero. Real heat engines generate a finite power output , i.e., they perform work during a cycle of a finite duration . Thus an appropriate way to characterize the engines is to compare their efficiencies at maximum power. On the macroscopic level this quantity is roughly bounded by the Curzon-Ahlborn value . Alternative expressions for quantifying efficiency have been discussed  which are based on mean quantities, e.g., on the mean work done during the operational cycle. On the mesoscopic level, the work is inherently a fluctuating quantity and one should be able to calculate not only its mean value but also its fluctuation properties.
In this paper we study a simple model of mesoscopic heat engines operating between two different heat baths under non-equilibrium conditions. The working medium consist of a two-level system. The cycle of operation includes just two isothermal branches, or strokes. Within each stroke, the system is driven by changing the energies of the two states and we assume a constant driving rate, i.e. a linear time-dependence of the energies. The response of the working medium is governed by a master equation with time-dependent transition rates. The specific form of the rates guarantees that, provided the two energies were fixed, the system would relax towards the Gibbs equilibrium state. Of course, during the motor operation, the Gibbs equilibrium is never achieved because the energies are cyclically modulated. At a given instant, the system’s dynamics just reflects the instantaneous position of the energy levels. After a transient regime, the engine dynamics approaches limit cycle with the periodicity of the driving force. We will focus on the properties of this limit cycle. In particular, we calculate the distribution of the work during the limit cycle.
Our two-isotherm setting imposes one important feature which is worth emphasizing. As stated above, at the end of each branch we remove the present bath and we allow the thermal interaction with another reservoir. This exchange of reservoirs necessarily implies a finite difference between the new reservoir temperature and the actual system (effective) temperature. Even if the driving period tends to infinity, we shall observe a positive entropy production originating from the relaxation processes initiated by the abrupt change of the contact temperature. Differently speaking, our engine operates in an inherently irreversible way and there exists no reversible limit.
The paper is organized as follows. In section 2 we solve the dynamical equation for the externally driven working medium. For the sake of clarity we first give the solution just for an unrestricted linear driving protocol using a generic driving rate and a generic reservoir (section 2.1). Thereupon, in section 2.2, we particularize the generic solution to individual branches and, using the Chapman-Kolmogorov condition, we derive the solution for the limit cycle. In section 3 we employ the recently derived  analytical result for the work probability density under linear driving. Again, we first give the result for the generic linear driving and then we combine two such particular solutions into the final work distribution valid for the limit cycle. The results from section 2 and section 3 enable a detailed calculation of the energy and entropy flows during the limit cycle in section 4 and allow for a discussion of the engine performance in section 5.
2 Description of the engine and its limit cycle
Consider a two-level system with time-dependent energies , , in contact with a single thermal reservoir at temperature . In general, the heat reservoir temperature may also be time-dependent. The time evolution of the occupation probabilities , , is governed by the master equation  with time-dependent transition rates specified by the reservoir temperature and by the external parameters. To be specific the dynamics of the system is described by the time inhomogeneous Markov process assuming the value , , if the system resides at time in the th state. Explicitly, the master equation reads
where is the unity matrix and the transition matrix with elements , . These elements are the conditional probabilities
If we denote by the initial state at time with the occupation probabilities , the occupation probabilities at the observation time are described by the column vector .
Due to the conservation of the total probability the system (1) can be reduced to just one non-homogeneous linear differential equation of the first order. Therefore the master equation (1) is exactly solvable for arbitrary functions , . The rates are typically a combination of an attempt frequency to exchange the state multiplied by an acceptance probability. We shall adopt the Glauber form
where sets the elementary time scale, and . The rates in equation (1) satisfy the (time local) detailed balance condition.
The resulting propagator satisfies the Chapman-Kolmogorov condition
for any intermediate time . Its validity can be easily checked by direct matrix multiplication. The condition simply states that the initial state for the evolution in the time interval can be taken as the final state reached in the interval . This is true even if the parameters of the process in the second interval differ from those in the first one. Of course, if this is the case, we should use an appropriate notation which distinguishes the two corresponding propagators. This procedure will be actually implemented in the paper. Keeping in mind this possibility, we shall first analyze the propagator for a generic linear driving protocol.
2.1 Generic case – linear driving protocol
Let us consider the linear driving protocol , and , where denotes the energy of the first level at the initial time , and is the driving velocity (energy change per time). The rates (3) can then be written in the form
where is the temperature-reduced driving velocity, and incorporates the initial values of the energies.
Under this linear driving protocol one can evaluate the definite integral in (4) explicitly and rewrite the propagator as
Here we have introduced the dimensionless ratio of the attempt frequency characterizing the time scale of the system’s dynamics and the time scale of the external driving, respectively. Naturally, this ratio will describe the degree of irreversibility of the process. Depending on the value , the explicit form of the function reads
where denotes the Gauss hypergeometric function .
2.2 Piecewise linear periodic driving
We now introduce the setup for the operational cycle of the engine under periodic driving. Within a given period, two branches with linear time-dependence of the state energies are considered with different velocities. Starting from the value , the energy linearly increases in the first branch until it attains the value , at time and in the second branch, the energy linearly decreases towards its original value in a time (see figure 1). We always assume , i.e.
This pattern will be periodically repeated, the period being .
As the second ingredient, we need to specify the temperature schedule. The two-level system will be alternately exposed to a hot and a cold reservoir, which means that the function in equation (3) will be a piecewise constant periodic function. During the first branch, it assumes the value , during the second branch it attains the value .
This completes the description of the model. Any quantity describing the engine’s performance can only depend on the parameters , , , , and . In the following we will focus on the characterization of the limit cycle, which the engine will approach at long times after a transient period.
Here the matrixes evolve the state vector within the respective branches and have the form
Notice, both propagators and are given by the generic propagator (4). In order to get , we replace in equation (5) the initial position of the first energy by , the driving velocity by , and we set . Analogously, the propagator follows from the generic propagator, if we replace by , by , and by .
The system state probabilities at the ends of the periods form a Markov chain and we are interested in its fixed point behavior. If we take the stationary state as the initial condition, the system revisits this special state at the end of the limit cycle. Therefore it suffices to solve the eigenvalue problem . Solving the algebraic equation, the fixed point probabilities at the beginning (or end) of the limit cycle are
These probabilities, and hence also the specific form of the limit cycle, depend solely on the model parameters.
We now put aside the transitory regime and we focus entirely on the limit cycle. Generally speaking, the parametric plot of the occupation difference (the response) versus the energy of the first level (the driving) exhibits two possible forms which are exemplified in figure 1. First, we have a one-loop form which is oriented either clockwise or anticlockwise. For clockwise orientation, the work done by the engine on the environment during the limit cycle is negative, while for counter-clockwise orientation it is positive. Secondly, we can obtain a two-loops shape exhibiting again either positive or negative work on the environment. Slowing down the driving, the branches gradually approach the corresponding equilibrium isotherms . We postpone the further discussion to the section 4.
3 Probability densities for work and heat
Heuristically, the underlying time-inhomogeneous Markov process can be conceived as an ensemble of individual realizations (sample paths). A realization is specified by a succession of transitions between the two states. If we know the number of the transitions during a path and the times at which they occur, we can calculate the probability that this specific path will be generated. A given paths yields a unique value of the microscopic work done on the system. For example, if the system is known to remain during the time interval , , in the th state, the work done on the system during this time interval is simply . Accordingly the probability of the paths gives the probability of the work. Viewed in this way, the work itself is a stochastic process and we denote it as . We are interested in its probability density , where denotes an average over all possible paths.
As a technical tool we introduce the augmented process which simultaneously reflects both the work variable and the state variable of the Markov process. The augmented process is again a time non-homogeneous Markov process. Actually, if we know at a fixed time both the present state variable and the work variable , then the subsequent probabilistic evolution of the state and work is completely determined. The work done during a time period , where , simply adds to the present work and it only depends on the succession of the states after the time . And this succession by itself, as we know from section 2, cannot depend on the dynamics before time .
The one-time properties of the augmented process will be described by the functions
where . We represent them as the matrix elements of a single two-by-two matrix ,
Using this matrix notation, the Chapman-Kolmogorov condition for the augmented process assumes the form
Here the matrix multiplication on the right hand side amounts for the summation over the intermediate states at the time , and the integration runs over all possible intermediate values of the work variable . The equation is valid for any intermediate time . Similarly to the preceding section 2, the Chapman-Kolmogorov condition can be used to connect two different propagators describing the time evolution of the augmented process within two branches of the driving cycle. Before we address this point, we focus on the generic situation.
We need an equation which controls the time dependence of the propagator and which plays the same role as the rate equation (1) in the case of the simple two-state process. It reads
where the initial condition is . This is a hyperbolic system of four coupled partial differential equations with time-dependent coefficients. It can be derived in several ways. For example, as explained in reference , one considers at the time the family of all realizations, which display at that time the work in the infinitesimal interval and, simultaneously, which occupy a given state. During the infinitesimal time interval , the number of such paths can change due to two reasons. First, while residing in the given state, some paths enter (leave) the set, because the energy levels move and an additional work has been done. Secondly, some paths can enter (leave) the described family because they jump out of (into) the specified state. These two contributions correspond to the two terms on the right hand side of equation (LABEL:Gdynequation). Another derivation  is based on an explicit probabilistic construction of all possible paths and their respective probabilities.
Similar reasoning holds for the random variable describing the heat accepted by the system from the environment, and for the internal energy . The variable is described by the propagator with the matrix elements
It turns out that there exists a simple connection between the heat propagator and the work propagator . Since for each path, heat and work are connected by the first law of thermodynamics, we have for any path which has started at time in the state and which has been found at time in the state . Accordingly,
where . This relation can be written in the form of the symmetry relation
3.1 Generic case–linear driving protocol
For the linear driving protocol the first term in the curly brackets in equation (LABEL:Gdynequation) is time-independent. As for the second term, we use again the Glauber rates (7). Thereby the evolution equation (LABEL:Gdynequation) assumes the form
with the parameters and introduced in connection with equation (7).
We shall now employ the method described in reference  by taking the double Laplace transformation with respect to the variables and . As shown in reference , a special difference equation results, which can be solved exactly. Moreover, it is possible to carry out the final double inverse Laplace transformation. However, in , only the case has been studied. In the present context we need the solution for a general initial energy difference . It turns out that such generalization represents a nontrivial task. The constant cannot be simply scaled off because it enters only the second term on the right hand side of equation (LABEL:Gdynequation). In order to overcome this difficulty, we had to modify the procedure from reference . However, in view of the specific topic of the present paper, we refrain from giving the technical details and we proceed with the description of the final result.
For the presentation of the result it is convenient to introduce the reduced work variable and the reduced time variable . Moreover, it is helpful to use the abbreviations
For , the result is
Here is the Dirac delta-function, and is the Heaviside unit step function. The solution for follows from interchanging the indices and in equations (27)-(30). If , then , and our results coincide with the formulae (49)–(52) in reference .
3.2 Piecewise linear periodic driving
The generic result (27)-(30) immediately yields the work and heat propagators for the individual branches in the protocol according to equation (11). We simply carry out the replacements described in the text following equation (15). We denote the corresponding matrices as and . Then the Chapman-Kolmogorov condition (20) yields the propagator
As demonstrated above, the heat propagator for the limit cycle is connected with the work propagator through simple shifts of the independent variable . Specifically, we get with , , and
In the last step we take into account the initial condition (17) at the beginning of the limit cycle and we sum over the final states of the process . Then the probability density for the work done on the system during the limit cycle reads
Similarly, the probability density for the head accepted within the limit cycle is
These two functions represent the main results of the present Section. They are illustrated in figures 2-4. We discuss their main features in section 5.
4 Engine performance
As shown in section 2, the occupation probabilities during the limit cycle are with given by equation (12). These probabilities are ensemble averaged quantities and cannot describe fluctuations of the engine’s performance. But they render the energetics in terms of mean values as we discuss now.
During the limit cycle, the internal energy changes as
Here is the mean heat received from the reservoirs during the period between the beginning of the limit cycle and the time . Analogously is the mean work done on the system from the beginning of the limit cycle till the time . If , the positive work is done by the system on the environment. Therefore the oriented areas enclosed by the limit cycle in figure 1 and in figure 4 represent the work done by the engine on the environment per cycle. These areas approach maximal absolute values in the quasi-static limit. The internal energy, being a state function, fulfills . Therefore, if the work is positive, the same total amount of heat has been transferred from the two reservoirs during the cycle. The case cannot occur if both reservoirs would have the same temperature. That the perpetuum mobile is actually forbidden can be traced back to the detailed balance condition in (1).
We denote the system entropy at time as , and the reservoir entropy at time as . They are given by
Upon completing the cycle, the system entropy re-assumes its value at the beginning of the cycle. On the other hand, the reservoir entropy is controlled by the heat exchange. Owing to the inherent irreversibility of the cycle we observe always a positive entropy production per cycle, . The total entropy increases for any . The rate of the increase is the larger the stronger is the representative point in the diagram deviates from the corresponding equilibrium isotherm (a strong deviation, e.g., can be seen in the diagram in figure 4c). Due to the instantaneous exchange of the baths at times and in the model considered here, a strong increase of always occurs after these time instants. A representative example of the overall behavior of the thermodynamic quantities (mean work and heat, and entropies) during the limit cycle is shown in figure 5.
An important characteristics of the engine is its power output and its efficiency . They are defined as
where is the total heat absorbed by the system per cycle. The performance of the engine characterized by the output work, efficiency, output power, and entropies from equations (37) and (38) are shown in figure 6 and figure 7.
In figure 6 the performance is displayed as a function of the cycle duration for . With increasing , the output work and the efficiency increase whereas the output power and the entropy production first increase up to a maximum and thereafter they decrease when approaching the quasi-static limit (). Notice that the maximum efficiency and output power occur at different values of . In figure 6a) we show also the standard deviation of the output work, which was calculated from the work probability density . Finally, let us note that the values and used in figure 6 give the Carnot efficiency . This should be compared with the efficiency of the engine for a long period , that is, with the value . As discussed above, the Carnot efficiency cannot be reached here even for , due to the immediate temperature changes at times and .
In figure 7 we have fixed and plotted the behavior as function of the time asymmetry (or time splitting) parameter . As can be seen from the upper three panels in figure 7, there exist also a maximal efficiency and a maximal output power with respect to a variation of the time asymmetry parameter (as long as the engine performs work, i.e., ). Again, the optimal parameter , where these maxima occur, is different for the efficiency and output power. In a reversed situation, considered in the lower three panels in figure 7, where the work is performed on the engine (), minima of the efficiency and output power occur.
The overall properties of the engine critically depend on the two dimensionless parameters . We call them reversibility parameters. For a given branch, say the first one, the parameter represents the ratio of two characteristic time scales. The first one, , is given by the attempt rate of the internal transitions . The second scale is proportional to the reciprocal driving velocity. Contrary to the first scale, the second one is fully under the external control. Moreover, the reversibility parameter is proportional to the absolute temperature of the heat bath.
Let us first consider the work probability density (34) within the first stroke in the case , c.f. figure. 2a). In essence, is given by a linear combination of the functions equations (27)-(30). It vanishes outside the common support which broadens linearly in time. Besides the continuous part located within the support, the diagonal elements display a singular part represented by delta functions at the borders of the support. The delta functions correspond to the paths with no transitions between the states. Specifically, the weight of the delta function located at represents the probability that the system starts in the first state and remains there up to time . The weight corresponding to the first level decreases with increasing time and vanishes for . On the contrary, the weight of the delta function at approaches the nonzero limit for , which is the probability that the a path starts in the second state and never leaves it.
Within the second stroke, the density results from the integral of the propagators for the individual strokes, c.f. equation (31). Due to the integration, the singular parts of the cycle propagator are now situated inside the support, at the values and . The two delta functions approach each other and, upon completing the cycle, they coincide at the point . The nonsingular component of the density is no more continuous . The jumps are located at the positions of the delta functions and their magnitudes correspond to the weights of the delta functions (for a discussion of the origin of these jumps, see ).
If both reversibility parameters are small, the isothermal processes during both branches strongly differs from the equilibrium ones. The signature of this case is a flat continuous component of the density and a well pronounced singular part. The strongly irreversible dynamics occurs if one or more of the following conditions hold. First, if is small, the transitions are rare and the occupation probabilities of the individual energy levels are effectively frozen during long periods of time. Therefore they lag behind the Boltzmann distribution which would correspond to the instantaneous positions of the energy levels. More precisely, the population of the ascending (descending) energy level is larger (smaller) than it would be during the corresponding reversible process. As a result, the mean work done on the system is necessarily larger than the equilibrium work. Secondly, a similar situation occurs for large driving velocities . Due to the rapid motion of the energy levels, the occupation probabilities again lag behind the equilibrium ones. Thirdly, the strong irreversibility occurs also in the low temperature limit. In the limit , the continuous part vanishes and .
In the opposite case of large reversibility parameters , both branches in the plane are located close to the reversible isotherms. The singular part of the density is suppressed and the continuous part exhibits a well pronounced peak. From general considerations , the density must approach a Gaussian shape. Our results allow a detailed study of this approach. Let us denote as the free energy of a two level system with energies at temperature , i.e.,