Thermodynamics of a continuous quantum heat engine: Interplay between population and coherence

We present a detailed thermodynamic analysis of a three-level quantum heat engine coupled continuously to hot and cold reservoirs. The system is driven by an oscillating external field and is described by the Markovian quantum master equation. We use the general form of the dissipator which is consistent with thermodynamics. We calculate the heat, power, and efficiency of the system for the heat-engine operating regime and also examine the thermodynamic uncertainty relation. The efficiency of the system is strongly dependent on the structure of the dissipator, and the correlations between different levels can be an obstacle for ideal operation. In quantum systems, the heat flux is decomposed into the population and coherent parts. The coherent part is specific to quantum systems, and in contrast to the population part, it cannot be expressed by a simple series expansion in the linear-response regime. We discuss how the interplay between the population and coherent parts affects the performance of the heat engine.


I. INTRODUCTION
One of the main objectives in the field of quantum thermodynamics is to find the quantum signatures of heat engines. The thermodynamical description is applied even when the system has a few degrees of freedom, which allows us to study the quantum effects on heat and work from a microscopic point of view.
Since the seminal work by Scovil and Schulz-DuBois [1], the microscopic heat engines with a few discrete energy levels have been studied in many works [2][3][4][5][6][7][8]. An upper bound on the efficiency of the heat engine is given by the Carnot efficiency, as expected from the general arguments [9,10]. Experimental examples of quantum heat engines have been realized in a wide range of quantum systems such as trapped ions and ensembles of nitrogen-vacancy centers in diamond [11][12][13].
It is an important problem to answer the question of whether the quantum coherence enhances the power output of the heat engine. Several works have shown that the answer is positive [7,[14][15][16], which represents a promising feature for the design of microscopic devices.
Generally speaking, quantum effects arise when the density operator has off-diagonal components in the basis of the Hamiltonian operator. The quantum coherence is characterized by the off-diagonal parts, and the nonequilibrium entropy production can be decomposed into the classical population part and the quantum coherent part [17,18]. The Markovian dynamics is described by the Gorini-Kossakowski-Lindblad-Sudarshan (GKLS) equation [19][20][21]. By decomposing the dissipator into two parts, we can introduce the corresponding heat for each part [22].
In the present work, we revisit the three-level quantum heat engine developed by Geva and Kosloff [3] in order to examine the roles of the population and coherent parts of the heat flux. We also study how the result depends on the details of the dissipator part of the GKLS equation. The choice of the dissipator is important for finding thermodynamically consistent results, and we treat a possible general form of the dissipator. This paper is organized as follows. In Sec. II, we describe the settings of our model and define the heat, work, and efficiency according to the standard scenario. Section III summarizes the stationary solution of the model. We obtain the explicit forms of the heat, power, and efficiency. In Sec. IV, we introduce the concept of population heat currents and coherent heat currents to discuss how the quantum nature affects the performance. We also discuss the thermodynamic uncertainty relation (TUR) in Sec. V. The conclusion is summarized in Sec. VI.

A. GKLS equation
We define a continuous quantum heat engine by using a three-level maser as shown in Fig. 1. The system consists of three states |0 , |1 , and |2 , and each level has energy ω 0 , ω 1 , and ω 2 , respectively. We set ω 0 < ω 1 < ω 2 . Two of the states, |0 and |1 , are coupled to the cold reservoir with temperature T c = 1/β c , and |0 and |2 are coupled to the hot reservoir with T h = 1/β h , where T h > T c . We study the performance of the heat engine when the system is continuously coupled to the two heat reservoirs. The system is operated by applying a time-oscillating field with frequency ω. The field induces transitions between |1 and |2 , generating a heat flow, as we discuss below. The system Hamiltonian is given in the state basis aŝ where λ represents the field intensity. This Hamiltonian is diagonalized asĤ(t) = 2 n=0 ǫ n |ǫ n (t) ǫ n (t)|, and the energy eigenvalues are independent of t, as shown in Appendix A. We assume λ 2 < ω 10 ω 20 , where ω 10 = ω 1 − ω 0 and ω 20 = ω 2 − ω 0 , so that the order of the energy levels is unchanged, ǫ 0 < ǫ 1 < ǫ 2 , in the presence of the field.
The present setting is basically the same as discussed in related works [3][4][5][6]8]. A difference arises when we carefully treat the dissipation effect. The state of the system is described by the density operatorρ(t). We assume the Markovian dynamics, and the time evolution is described by the GKLS equation, The dissipation affects the quantum dynamics in three ways. Two of them are described by the dissipatorD α . It changes the eigenbasis of the density operator and the population of each basis element [22]. The remaining one shifts only the energy levels of the HamiltonianĤ(t) and is called the Lamb shift. Here we simply neglect the Lamb shift term or takeĤ(t) to be a renormalized Hamiltonian. The dissipatorD α is chosen so that the time evolution is a completely positive and trace-preserving map. The explicit form of the dissipator is dependent on coarsegraining procedures. In order to obtain a thermodynamically consistent description for the present problem, we use the form whereL ǫ α (t) represents projected jump operatorŝ The dissipator coupling γ α (ǫ) is generally nonnegative and is a system-dependent function of ǫ and β α . To construct a thermodynamically consistent theory, we assume the detailed balance condition Then, the Gibbs distribution becomes an instantaneous stationary solution [21,22]. We can derive the form of the dissipator in Eq. (3) from a microscopic model by using several approximations such as the Markov approximation and the rotatingwave approximation. In principle, the derivation implies that the present model is justified only within a certain range of parameters. However, the present setting without any additional constraints is completely consistent with the laws of thermodynamics, and we can discuss the performance of the quantum heat engine. The use of the jump operators projected onto the instantaneous eigenstates of the Hamiltonian is an important ingredient to find a thermodynamically consistent theory. It is contrasted to the "local" master-equation approach in which the jump operators are not projected to the eigenstate basis. The local approach is shown to be inconsistent with thermodynamics [23]. Although some ideas to overcome this shortcoming have been discussed [24,25], here we use the "global" approach defined by Eq. (3).
In principle, the present model has four types of the dissipator coupling, γ c (ǫ 10 ), γ c (ǫ 20 ), γ h (ǫ 10 ), and γ h (ǫ 20 ), where ǫ 10 = ǫ 1 − ǫ 0 and ǫ 20 = ǫ 2 − ǫ 0 . In previous works, only γ c (ǫ 10 ) and γ h (ǫ 20 ) were kept nonzero [3,4]. As we mentioned above, γ α (ǫ) is a system-dependent function and is obtained from the correlation function of a bath operator. A typical form is represented by using a Lorentzian function. In the following, we do not assume any functional form of γ α (ǫ) and consider three possible cases: The resonant-coupling case corresponds to the preceding works.

B. Heat, work, and efficiency
According to the standard scenario, we define the heat flux from the reservoirs to the system aṡ The last expression allows us to writeQ(t) = αQ α (t), whereQ α (t) represents the heat flux from each reservoir. When the system is operated periodically, the energy goes back to the original value after one period such that the work done by the system in the one period is given by where T 0 = 2π/ω. When the system acts as a heat engine, the relations Q h > 0, Q c < 0, and W > 0 hold, as we see in Fig. 2, and the efficiency is defined as The second law of thermodynamics is derived from the nonnegativity of the entropy production [10]. For the present model we obtain and as a result, the efficiency is bounded from above by the Carnot efficiency The nonnegativity of the entropy production holds even in the present model [22], which confirms the previous result on the upper bound on the efficiency [1][2][3][4][5][6].

A. Stationary solution
We use the stationary solution of the GKLS equation in the present setting to evaluate the performance of the heat engine. We note that the stationary solution means that the system settles down to a stable periodic behavior after transient evolutions during the first several periods. The eigenstate decomposition of the jump operator ǫ n (t)|L α |ǫ m (t) is shown to be time independent, which gives a simple stationary result. We describe the details in Appendix B. Here we summarize the result.
The stationary solution of the GKLS equation was studied in the preceding works. As we stressed above, the crucial difference is that the dissipator is represented by four types of dissipator couplings. Correspondingly, the explicit form of the dissipator is parametrized by four types of coupling functions, where θ is defined by the relation θ represents a rotation angle for the diagonalization of the Hamiltonian, as shown in Appendix A. In the resonant-coupling limit, relations g − 1 = e −βcǫ10 g 1 and g − 2 = e −β h ǫ20 g 2 hold, and the dissipator is characterized by g 1 and g 2 . In the general case, no trivial relations hold between g 1 , g 2 , g − 1 , and g − 2 , although we still have the detailed balance condition in Eq. (7).
We also introduce dimensionless nonnegative parameters to characterize thermodynamic quantities in the following. For the resonant-coupling case these quantities collapse to zero q 1 = q 2 = 0. At the stationary limit, the heat flux from each reservoir is given bẏ where ǫ mn = ǫ m − ǫ n and All the quantities introduced above are independent of t. ρ 0 and G are always positive, irrespective of the choice of the parameters. We show below that P represents the power of the heat engine and ρ 0 represents the groundstate component of the density operator, 0|ρ(t)|0 . The term ρ 0 P 0 represents a direct flow from the hot reservoir to the cold reservoir, and P 0 goes to zero when the temperature difference β c − β h disappears. Further details are discussed below.

B. Power and efficiency
Anticipating Q h > 0, Q c < 0, and W > 0, we can obtain the explicit form of the work done by the system. The power of the heat engine, which is defined by the work divided by the cycle period, is given by at the stationary limit. P is given in Eq. (25). The corresponding efficiency is obtained as where η SSD is reminiscent of the Scovil-Schulz-DuBois efficiency [1]. This expression is simplified when we consider the resonant-coupling limit γ c (ǫ 20 ) = γ h (ǫ 10 ) = 0. In this case, we find q 1 = 0, q 2 = 0, and P 0 = 0, and the efficiency coincides with the Scovil-Schulz-DuBois efficiency This result is consistent with the preceding works [2][3][4][5][6][7][8].
The efficiency in the general case is dependent on various parameters such as the temperatures and the frequency and is bounded from above by Eq. (33).

C. Heat engine conditions
The above results in the present section are exact and hold irrespective of the choice of parameters. When we require that the system works as a heat engine, the relations Q h > 0, Q c < 0, and P > 0 must hold.
As we see from Eq. (25), P > 0 holds when This condition is satisfied only when This is a necessary condition in general and is the necessary and sufficient condition in the resonant-coupling case. In the general case, Eq. (34) is rewritten as where The coupling-constant parameters in the dissipator are taken so that Eq. (36) is satisfied. We note that η SSD in Eq. (32) becomes positive in that case. The condition in Eq. (35) determines possible values of parameters in the Hamiltonian for a given η C as Equation (35) shows that the power becomes negative when the temperature difference is too small. In that case, the system does not work as a heat engine, and we can observe a heat flow from the low-temperature reservoir to the high-temperature one as Q h < 0 and Q c > 0. This behavior implies that the present system cannot be understood from the standard linear-response theory in which the heat flow arises due to the temperature difference. We discuss the origin of this quantum nature in the next section.

D. Plot of the results
At small frequency values, the power P is proportional to ω 2 . The efficiency η is also proportional to ω 2 provided P 0 > 0, as we see from Eq. (31). We see that the frequency-independent result in Eq. (33) is specific to the resonant coupling and is unusual. P as a function of ω has a Fano resonant form and is plotted in Fig 3. It is maximized at Interestingly, the efficiency η is also maximized at this frequency. The magnitude of the optimal frequency is basically determined by the energy gap ǫ 21 . The dissipation effect enhances the resonant frequency slightly. We plot the power P and the efficiency η for the optimal frequency ω in Eq. (41) in Fig. 4 (resonant coupling), We observe that the power is maximized at small ω 20 /ω 10 (> 1) and at a moderate value of λ. Although the possible parameter range for a heat engine is dependent on the dissipator couplings as well as the temperatures, the contour map is basically insensitive to the parameters.
In contrast to the power, the efficiency exhibits a stronger dependence on the dissipator coupling. As shown in Fig. 4(b) for the resonant coupling, the efficiency attains its maximum at the boundary where the power goes to zero. This behavior is totally reversed in comparison to the other dissipator coupling cases in Figs. 5(b) and 6(b), where the efficiency is minimized at the boundary. We also find that the efficiency and the power are reduced when we move away from the resonant coupling.
We plot the distribution (η, P ) in Figs. 4(c), 5(c), and 6(c). We find that the envelope curve of the distributions has a tendency towards being symmetric around the efficiency at maximum power as we approach the uniform coupling. We note that the efficiency at maximum power cannot be understood from Curzon-Ahlborn efficiency even in the linear-response regime [26,27]. As we mentioned above, the present system does not exhibit a heat-engine behavior at the linear response regime.
Summarizing the present result, we find that the performance of the present model as a heat engine worsens when we move away from the resonant coupling. In the next section, we discuss the origin of this behavior.

IV. DECOMPOSITION OF HEAT
Compared to the resonant coupling, we see that the decreasing of the efficiency in extended couplings is understood from the presence of ρ 0 P 0 in Eq. (31). As we also see in Eqs. (23) and (24), it represents a direct flow from the hot reservoir to the cold reservoir, which clearly reduces the performance of the heat engine. P 0 goes to zero when the temperature difference is zero and remains positive even at the zero-field limit λ → 0. This implies that the direct flow has a classical interpretation.
The quantum nature of the system can be understood by realizing that the basis of the density operator is different from that of the Hamiltonian. The density operator is diagonalized aŝ The eigenvalues of the density operator obey the masterequation-like relation and the basis of the density operator obeys the unitary time evolution whereξ(t) represents a generator of the time evolution [22]. The explicit forms of the eigenvalues and the eigenstates of the density operator are given in Appendix C. Accordingly, the heat flux in Eq. (11) is decomposed into the diagonal and nondiagonal parts asQ andQ nd α (t) is defined by the residual contribution oḟ Q α (t). The diagonal part is related to the population dynamics in Eq. (43), and the nondiagonal part is related to the coherent dynamics. In the present setting, the eigenvalue p n is independent of t. As a result, we obtainQ This relation shows that the diagonal part of the heat flux just goes through the system asQ d h = −Q d c and does not contribute to the work. This direct flow reduces the efficiency of the heat engine.
To quantify the diagonal and nondiagonal contributions, we decompose the efficiency as where η d denotes the contribution from the diagonal part of the heat flux, and η nd denotes the contribution from the nondiagonal part. Equation (46) shows that η d = 0 in the present case. The efficiency is given by and the explicit form of η nd is where . From Eq. (50), we can identify the diagonal and nondiagonal heat flows, respectively, as We plot 1/η nd as a function of λ in Fig. 7. From the value of 1/η nd , we can understand the pattern of the heat flow. When η nd > 1, both Q nd c and Q nd h are positive, and we can obtain the ideal heat flow as a heat engine. The case η nd > 1 is represented by pattern (ii) in Fig. 7. η nd has a purely quantum-mechanical origin and can exceed the Carnot efficiency. However, we cannot neglect the diagonal contribution, which reduces the total efficiency. As a result, even in the quantum system, the efficiency in Eq. (50) is bounded from above by the Carnot efficiency.
In the resonant-coupling case, η nd is always larger than unity, which leads to a high-performance result in Fig. 4. The situation is significantly changed when we move away from the resonant-coupling case. As we increase λ, the nondiagonal heat flow changes, and we observe the reduction of the efficiency as a result. Dependent on the parameters in the equation, we can observe pattern (i) and pattern (iii) in Fig. 7 within the heat-engine domain.
In Fig. 8, we plot the patterns of the heat flow in the uniform-coupling case. We can find patterns (i) and (iii) around the boundary where the efficiency becomes small. This behavior is consistent with that in Fig. 6.
The decomposition of the diagonal part and the nondiagonal part has been discussed in some works [17,18]. It is a difficult problem to observe each one as an independent contribution. However, the nondiagonal part goes to zero at λω → 0. The diagonal part is insensitive to the parameter and can be extracted in the weak-field regime.

V. THERMODYNAMIC UNCERTAINTY RELATION
As a final subject to study, we examine the TUR [28,29]. The standard form of the TUR is represented as where σ is the entropy production rate averaged over one cycle and var P is the variance of the power. The variance is bounded from below, and the bound is determined by the entropy production. Although this relation was shown in a broad range of classical systems, the relation does not necessarily hold, and the violation can be found especially in quantum systems. The quantum TUR is modified by a different bound [30]. For the present three-level system, the violation of the standard TUR was shown in [31,32]. Here we examine how this result is affected by the modification of the dissipator. The entropy production rate at each t is given bẏ σ(t) = −Tr ∂ tρ (t) lnρ(t) − α β αQα (t). The first term comes from the von Neumann entropy of the system and goes to zero when we take the average over one period. The average of the entropy production rate The variance of the power var P is calculated in Appendix D. It is decomposed as varP = (varP ) 1 + (varP ) 2 + (varP ) 3 + (varP ) 4 , and each part is respectively given as (varP ) 4 = 1 We note that (var P ) 1 and (var P ) 3 are nonnegative. We also see that (var P ) 2 is negative, and (var P ) 4 is zero at We plot var P in Fig. 9. The variance is basically an increasing function of λ/ω 10 and is insensitive to ω 20 /ω 10 . This result holds irrespective of the choice of the dissipator coupling. The corresponding behavior of the TUR is shown in Fig. 10. As we see Fig. 10, the TUR bound is strictly satisfied with the present choice of parameters. The bound is tight in the case of the resonant coupling and is loosened as we move away from the resonant coupling.
All of the results in Figs. 9 and 10 satisfy the standard TUR in Eq. (54). The optimal frequency in Eq. (41) is used there, and the result is changed by considering the frequency dependence in Fig. 11. We observe a violation These numerical results can be understood from the analytical expression. Basically, the bound comes from Eq. (56). In the resonant-coupling case, we find We numerically find that the other parts are very small, at least with the present choice of parameters. The small violation of the TUR in Fig. 11 is understood from a negative contribution of (var P ) 2 . We note that the contributions (var P ) 1 and (var P ) 2 are similar to the result in Ref. [31], where the variance takes the form A and B are positive, and the second term leads to a violation of the TUR in a certain range of parameters. We note that the local approach is used in Ref. [31] and the form of the dissipator is different from the present model. and (var P ) 4 , but we numerically find that these contributions are small and do not play any significant role.
In the case of the other dissipator coupling, we can understand the loose bound from the expression of the entropy production rate. The last term in Eq. (55) includes P 0 , which makes the uncertainty product σ var(P )/P 2 very large.

VI. CONCLUSIONS
We have presented a detailed thermodynamic analysis of a continuous quantum heat engine based on the global form of the quantum master equation. We found that the performance of the heat engine is strongly dependent on the form of the dissipator. The quantum coherence does not necessarily enhance the efficiency of the heat engine.
The quantum coherent heat flow cannot be understood from the laws of thermodynamics. It produces a nontrivial heat flow even in the linear-response regime. Although the coherent flow has the potential ability to enhance the efficiency of the heat engine, the efficiency is affected by the population heat flow, which prevents the system from violating the second law of thermodynamics. In the present model, the population heat flow does not contribute to the power of the heat engine but is required to keep the system within the heat-engine operation regime.
The decomposition of the efficiency in Eq. (47) allowed us to discuss the interplay between the population and coherent parts. We compared η d and η nd to find which part is important for the system to work as a heat engine. In the present setting, we found η d = 0. Microscopically, this property is due to the time independence of the eigenvalues of the Hamiltonian. When we consider a more general form of the Hamiltonian, the eigenvalue is generally dependent on t, and the population heat flow contributes to the power generation. It will be an interesting problem to study the performance of the heat engine and the interplay between the population and coherent parts in that case.
We obtain the stationary solution of the GKLS equation in Eq. (2). Multiplying ǫ m (t)| from the left and |ǫ n (t) from the right in Eq. (2), we obtain the expression Due to the property |ǫ 0 (t) = |0 and our choice of the jump operators in Eqs. (5) and (6), 0|ρ(t)|ǫ 1 (t) and 0|ρ(t)|ǫ 2 (t) , their conjugates do not couple to the other components and decay exponentially as a function of t. We write the off-diagonal component ǫ 1 (t)|ρ(t)|ǫ 2 (t) = e iωt (∆ 1 (t) + i∆ 2 (t)) , where ∆ 1 and ∆ 2 are real, to obtain ∂ t 0|ρ|0 = g 1 ǫ 1 |ρ|ǫ 1 + g 2 ǫ 2 |ρ|ǫ 2 − (g − 1 + g − 2 ) 0|ρ|0 , (B3) ∂ t ǫ 1 |ρ|ǫ 1 = −g 1 ǫ 1 |ρ|ǫ 1 + g − 1 0|ρ|0 − ω∆ 2 sin θ, (B4) ∂ t ǫ 2 |ρ|ǫ 2 = −g 2 ǫ 1 |ρ|ǫ 1 + g − 2 0|ρ|0 + ω∆ 2 sin θ, (B5) Due to the normalization of the density operator, the first three equations are not independent from each other. We can write setting χ α = χ and by using the stationary solution, we obtain This result is consistent with that in Eq. (30). Next, we consider the trace of the GKLS equation at second order given by In order to solve the equation, we need to find the explicit form ofρ χ 1 , not of Trρ χ 1 . The GKLS equation at first order in χ is given by where L is given in Eq. (B8). As we have shown above, ρ χ 1 is a linear function of t. This means that the first component of the vector on the right-hand side of Eq. (D8) is proportional to t and the other components are independent of t. By using Eqs. (B7) and (D4), we can write the solution as where This solution is inserted into Eq. (D7). Trρ χ 2 (t) is a quadratic function of t. The t 2 term represents the disconnected part of the fluctuation and is canceled out by subtracting the square of the average as in Eq. (D3). We obtain the expression of the variance given in the main text.