Quantum Otto engine with quantum correlations

We consider a quantum engine cycle that uses a single mode of the quantized radiation within an optical cavity as its working substance and operates between hot and cold heat reservoirs. During the cold isochoric process, the thermal bath is composed of an infinite collection of noninteracting bosonic modes, but during the hot isochoric stroke, the optical cavity interacts resonantly with a beam composed of thermally entangled pairs of two-level atoms [see Fig. 2(a)] that play the role of the hot nonthermal reservoir. The interaction of the cavity with these atoms is realized by sending only one of a pair of atoms to pass through the cavity, which means that only one of the atoms in the pair interacts with the radiation field Lutz09 . The engine model under consideration as a quantum version of the Otto cycle is drawn schematically in Fig. 2(b). This model consists of two unitary strokes denoted by $A\rightarrow B$ and $C\rightarrow D$, where the system is isolated from the two heat reservoirs, and two isochoric branches denoted by $B\rightarrow C$ and $D\rightarrow A$, along which the optical cavity with the constant Hamiltonian is weakly coupled to the hot reservoir or the cold reservoir, respectively. We now describe the four consecutive steps of the proposed Otto engine cycle as follows.

(i) Unitary compression $A\rightarrow B$: The system is isolated from the heat reservoir and undergoes a unitary compression in time $\tau_{ch}$. The system’s Hamiltonian is changed from $H^{ca}(\omega_{c})=\omega_{c}(\hat{a}^{{\dagger}}\hat{a}+1/2)$ to $H^{ca}(\omega_{h})=\omega_{h}(\hat{a}^{{\dagger}}\hat{a}+1/2)$ with a driving function $\omega_{ch}(t)=\omega_{c}\omega_{h}\tau_{ch}/[(\omega_{c}-\omega_{h})t+\omega_{h}\tau_{ch}]$, where $\hat{a}$ and $\hat{a}^{{\dagger}}$ are the annihilation and creation operators of the oscillator, respectively. Because no heat is exchanged, the change in the system’s internal energy is equal to the work done on the system. We use a two-projective-measurement scheme Lutz20 to calculate the probability distribution of the quantum work as follows:

$$p({w}_{ch})=\mathop{\sum}\limits_{n,m}\delta[{w_{ch}}-(\varepsilon_{m}^{h}-\varepsilon_{n}^{c})]p_{n\rightarrow m}^{{ch}}p_{n}^{A},$$

(7)

where $\varepsilon_{n}^{c}$ and $\varepsilon_{m}^{h}$ are the measured energies at the beginning and the end of this stroke, respectively. Here, $p_{n}^{A}$ is the initial occupation probability and $p_{n\rightarrow m}^{{ch}}=|\langle n|U_{ch}|m\rangle|^{2}$ with the unitary operator $U_{ch}$ is the transition probability from state $|n\rangle$ to $|m\rangle$.

(ii) Isochoric heating $B\rightarrow C$: For this stroke, we consider the case in which only one atom of the correlated pair with the constant inverse temperature $\beta_{h}$ is sent to pass through the cavity. The Hamiltonian of the single radiation mode in the cavity is constant, with a fixed frequency $\omega_{h}$. Therefore, the stochastic heat injection is equivalent to an increase in the system eigenenergy, with no work being produced along the stroke. The probability distribution for the heat absorbed by the system during this stroke, given that the stochastic compression work is $w_{ch}=\varepsilon_{m}^{h}-\varepsilon_{n}^{c}$, can then be expressed as Lutz20 ; jiao21

$$p({q}_{h}|{w}_{ch})=\mathop{\sum}\limits_{k,l}\delta[q_{h}-(\varepsilon_{l}^{h}-\varepsilon_{k}^{h})]p_{k\rightarrow l}^{{}_{h}}p_{k}^{B},$$

(8)

where $\varepsilon_{k}^{h}$ and $\varepsilon_{l}^{h}$ are the measured energies at the beginning and the end of the heating stroke, respectively. $p_{k}^{B}=\delta_{km}$ indicates the probability that the system is in eigenstate $|m\rangle$ after the second projective measurement during the unitary compression, and $p_{k\rightarrow l}^{{h}}$ is the probability that the system will collapse into the state $|l\rangle$ after the second projective measurement for this hot isochore. In the case where the system reaches the unique steady state at the end of the isochore Lutz20 ; jiao21 , the probability satisfies the following generalized canonical form: $p^{h}_{k\rightarrow l}=e^{-\beta_{h}^{\mathrm{eff}}H^{ca}(\omega_{h})}/\mathrm{Tr}[e^{-\beta_{h}^{\mathrm{eff}}H^{ca}(\omega_{h})}]$, where $\beta_{h}^{\mathrm{eff}}$ is called the effective temperature and will be defined in Eq. (11).

The atom that passes through the optical cavity interacts with the cavity via a resonant Jaynes-Cummings coupling, which is given by $H^{int}=-\gamma(\hat{a}\sigma^{+}+\hat{a}^{{\dagger}}\sigma^{-})$ Nar80 with a coupling constant $\gamma$. When an ensemble composed of many atoms is sent to the optical cavity, the cavity state dynamics can be described by Lutz09 ; Sch01

	$\displaystyle\frac{d\rho_{t}}{dt}$	$\displaystyle=$	$\displaystyle{{i}[H^{ca}(\omega_{h}),\rho_{t}]}+r_{1}^{h}(\gamma\tau)^{2}(\hat{a}^{\dagger}\rho_{t}\hat{a}-\frac{1}{2}\hat{a}\hat{a}^{\dagger}\rho_{t}-\frac{1}{2}\rho_{t}\hat{a}\hat{a}^{\dagger})$		(9)
		$\displaystyle+$	$\displaystyle r_{2}^{h}(\gamma\tau)^{2}(\hat{a}\rho_{t}\hat{a}^{\dagger}-\frac{1}{2}\hat{a}^{\dagger}\hat{a}\rho_{t}-\frac{1}{2}\rho_{t}\hat{a}^{\dagger}\hat{a}),$

where $\tau$ is the time spent by the atoms inside the cavity and $H^{ca}(\omega_{h})=\omega_{h}(\hat{a}^{\dagger}\hat{a}+\frac{1}{2})$ is the Hamiltonian of the single-mode cavity. The coefficients $r_{1}^{h}$ and $r_{2}^{h}$ are the arrival rates for the atoms in the excited and ground states, respectively, and these coefficients are thus associated with the probabilities of emission and absorption of a photon in the cavity, respectively. These two coefficients are determined using $r_{1}^{h}=\rho_{e}^{h}+\rho_{d}^{h}$ and $r_{2}^{h}=\rho_{g}^{h}+\rho_{d}^{h}$, where $\rho_{e}^{h}=\mathrm{exp}(-\beta_{h}\omega_{h})/[2\mathrm{cosh}(\beta_{h}\omega_{h})+2\mathrm{cosh}(\beta_{h}\xi)]$, $\rho_{g}^{h}=\mathrm{exp}(\beta_{h}\omega_{h})/[2\mathrm{cosh}(\beta_{h}\omega_{h})+2\mathrm{cosh}(\beta_{h}\xi)]$, and $\rho_{d}^{h}=\mathrm{cosh}(\beta_{h}\xi)/[2\mathrm{cosh}(\beta_{h}\omega_{h})+2\mathrm{cosh}(\beta_{h}\xi)]$.

The optical cavity, as the working substance, is allowed to relax to the stationary state at the end of the hot isochoric stroke, where the first term on the right hand side of Eq. (9) is vanishing because $[H^{ca}(\omega_{h}),{\rho_{t}}]=0$, and the time duration $\tau_{h}$ of this stroke is assumed to be much shorter than the duration of the unitary stroke. The asymptotic steady-state solution to Eq. (9), which describes the stationary state of the optical cavity, is then obtained as:

$$\rho^{ss}=(e^{\beta_{h}^{\mathrm{eff}}\omega_{h}/2}-e^{-\beta_{h}^{\mathrm{eff}}\omega_{h}/2})e^{-\beta_{h}^{\mathrm{eff}}H^{ca}(\omega_{h})},$$

(10)

where $\beta_{h}^{\mathrm{eff}}$ is introduced to denote the effective inverse temperature of the optical cavity. The detailed balance condition $\mathrm{exp}(-\beta_{h}^{\mathrm{eff}}\omega_{h})=r_{1}^{h}/r_{2}^{h}$ then produces

$${\beta_{h}^{\mathrm{eff}}}=\beta_{h}-\frac{1}{\omega_{h}}\mathrm{ln}\frac{1+e^{\beta_{h}\omega_{h}}\mathrm{cosh}(\beta_{h}\xi)}{e^{\beta_{h}\omega_{h}}+\mathrm{cosh}(\beta_{h}\xi)}.$$

(11)

Given the frequency $\omega_{h}$ and the interaction strength $\xi$, the effective inverse temperature $\beta_{h}^{\mathrm{eff}}$ will thus increase monotonically with increasing inverse temperature $\beta_{h}$, as illustrated in Fig. 3(a). The normalized temperature $\beta_{h}^{\mathrm{eff}}/\beta_{h}$ as a function of $\xi$ is plotted in Fig. 3(b), demonstrating that $\beta_{h}^{\mathrm{eff}}/\beta_{h}$ is a monotonically decreasing function of $\xi$ and $\beta_{h}^{\mathrm{eff}}/\beta=1$ for a vanishing $\xi$, as expected. The discord $\mathcal{Q}(\rho_{12})$ given by Eq. (6) is determined by the interaction strength $\xi$, which means that the effective temperature $\beta_{h}^{\mathrm{eff}}$ is closely dependent on the discord $\mathcal{Q}(\rho_{12})$. Because the discord is a monotonically increasing function of the interaction strength, the ratio $\beta_{h}^{\mathrm{eff}}/\beta_{h}$ decreases with increasing $\mathcal{Q}(\rho_{12})$ [see the inset of Fig. 3(b)]. Physically, a greater discord means that the system deviates further away from the thermal state of the inverse temperature $\beta$. We also note from Eq. (11) that the expression for the normalized temperature $\beta_{h}^{\mathrm{eff}}/\beta_{h}$ can be simplified to give ${\beta_{h}^{\mathrm{eff}}}/{\beta_{h}}\simeq 1-(\beta_{h}\xi)^{2}/4=1-2\mathcal{Q}(\rho_{12})\mathrm{ln}2$ in the high-temperature and/or weak-coupling limit .

(iii) Unitary expansion $C\rightarrow D$: The optical cavity is isolated again from the heat reservoir and undergoes a unitary expansion during the time period $\tau_{hc}$ with the driving function $\omega_{hc}(t)=\omega_{c}\omega_{h}\tau_{hc}/[(\omega_{h}-\omega_{c})(t-\tau_{ch})+\omega_{c}\tau_{hc}]$. Because this stroke can be accomplished by reversing the protocol used in the compression process described above, we set $\tau_{\mathrm{dri}}\equiv\tau_{hc}=\tau_{ch}$ to obtain $\omega_{hc}(t)=\omega_{ch}(2\tau_{\mathrm{dri}}-t)$. During this stroke, no heat is exchanged, and the change in the system’s internal energy is equal to the work done on the system. For the given quantum expansion work $w_{ch}$ and heat injection $q_{h}$, by using the two-point-measurement scheme again, the probability distribution of the stochastic work during the expansion process is given by

$$p({w_{hc}}|{w_{ch}},{q}_{h})=\mathop{\sum}\limits_{i,j}\delta[{w_{hc}}-(\varepsilon_{j}^{c}-\varepsilon_{i}^{h})]p_{i\rightarrow j}^{{hc}}p_{i}^{C},$$

(12)

where $p_{i}^{C}=\delta_{il}$ is the probability that the optical cavity is in the eigenstate $|l\rangle$ after the second projective measurement in the isochoric heating. In addition, $p_{i\rightarrow j}^{hc}=|\langle j|U_{hc}|i\rangle|^{2}$ is the transition probability between the instantaneous eigenstates $|i\rangle$ and $|j\rangle$, where $U_{hc}$ is the unitary operator.

(iv) Isochoric cooling $D\rightarrow A$: The system, which has the constant Hamiltonian $H^{ca}(\omega)=H^{ca}(\omega_{c})$, is coupled to a thermal reservoir of constant inverse temperature $\beta_{c}$ in the time duration $\tau_{c}$, which is much shorter than the adiabatic driving time $\tau_{\mathrm{dri}}$. The system reaches thermal equilibrium with the heat reservoir at the end point of the cold isochore, and the state of the system at this end point is given by $\rho_{0}=e^{-\beta_{c}H^{ca}(\omega_{c})}/\mathrm{Tr}[e^{-\beta_{c}H^{ca}(\omega_{c})}]$.

After a single cycle, the joint distribution of the total work ${w}$ and the heat injection ${q}_{h}$ can be calculated by combining Eqs. (7), (8), and (12) to give Lutz20 ; Hol121 ; Lutz21 ; jiao21

	$\displaystyle p({w},{q_{h}})$	$\displaystyle=$	$\displaystyle\mathop{\sum}\limits_{n,m,i,j}\delta({w}+\varepsilon_{n}^{c}-\varepsilon_{j}^{c}+\varepsilon_{i}^{h}-\varepsilon_{m}^{h})$		(13)
		$\displaystyle\times$	$\displaystyle\delta({q_{h}}-\varepsilon_{i}^{h}+\varepsilon_{m}^{h})|\langle m|U_{ch}|n\rangle|^{2}|\langle j|U_{hc}|i\rangle|^{2}$
		$\displaystyle\times$	$\displaystyle\frac{e^{-\beta_{c}\varepsilon_{n}^{c}}e^{-\beta_{h}^{\mathrm{eff}}\varepsilon_{i}^{h}}}{Z_{c}^{ca}Z_{h}^{ca}},$

where the partition functions of the optical cavity are $Z_{c}^{ca}=\mathrm{Tr}[e^{-\beta_{c}H^{ca}(\omega_{c})}]$ and $Z_{h}^{ca}=\mathrm{Tr}[e^{-\beta_{h}^{\mathrm{eff}}H^{ca}(\omega_{h})}]$.

To determine the statistics of the work and the heat, the characteristic function of the joint distribution function given in Eq. (13), denoted by $G(u,v)\equiv\int\int p(w,q_{h})e^{-ivq_{h}-iuw}dwdq_{h}$, can be obtained as follows Lutz21 ; Fei22 :

$$G(u,v)=\langle e^{-iv{q_{h}}-iuw}\rangle=\frac{G_{c}G_{h}}{1+\phi},$$

(14)

where

$$G_{c}=\frac{e^{\beta_{c}\omega_{c}}-1}{\sqrt{r_{\phi}\mathrm{cosh}(\beta_{c}\omega_{c}+iv_{0})+\mathrm{cosh}(\beta_{c}\omega_{c}+iu_{0})-2}},$$

(15)

and

$$G_{h}=\frac{e^{\beta_{h}^{\mathrm{eff}}\omega_{h}}-1}{\sqrt{r_{\phi}\mathrm{cos}(v_{0}-i\beta_{h}^{\mathrm{eff}}\omega_{h})+\mathrm{cos}(i\beta_{h}^{\mathrm{eff}}\omega_{h}+u_{0})-2}}.$$

(16)

Here, we have used $u_{0}=u(\omega_{c}-\omega_{h})+v\omega_{h}$, $v_{0}=u(\omega_{c}+\omega_{h})-v\omega_{h}$, and $r_{\phi}=(1-\phi)/(1+\phi)$. As shown in Appendix A, the parameter $\phi$ can be derived to be

$$\phi\equiv 1+\frac{1-\cosh\left(\sqrt{1-\zeta}\ln\frac{\omega_{h}}{\omega_{c}}\right)}{\zeta-1},$$

(17)

where $\zeta=[2\tau_{\mathrm{dri}}\omega_{c}\omega_{h}/(\omega_{h}-\omega_{c})]^{2}$, and this parameter is called the nonadiabatic factor. In the case where $\zeta>1$, Eq. (17) then becomes $\phi=1+\frac{1-\mathrm{cos}[\sqrt{\zeta-1}\mathrm{ln}(\omega_{h}/\omega_{c})]}{\zeta-1}$, indicating that the extreme limit of $\zeta\gg 1$ leads to $\phi\rightarrow 1$. Within the long time limit, $\phi$ approaches $1$, and thus the quantum adiabatic condition is satisfied .

The average and the variance of the work, and the average heat injection, can be obtained explicitly as

	$\displaystyle\langle{w}\rangle$	$\displaystyle=$	$\displaystyle-i\frac{\partial\mathrm{ln}G(u,v)}{\partial u}\bigg{|}_{u=v=0}$		(18)
		$\displaystyle=$	$\displaystyle\omega_{h}(\phi\langle n_{c}^{eq}\rangle-\langle n_{h}^{ss}\rangle)+\omega_{c}(\phi\langle n_{h}^{ss}\rangle-\langle n_{c}^{eq}\rangle),$

	$\displaystyle\delta w^{2}$	$\displaystyle=$	$\displaystyle\langle w^{2}\rangle-\langle w\rangle^{2}=-\frac{\partial^{2}\mathrm{ln}G(u,v)}{\partial u^{2}}\bigg{|}_{u=v=0}$
		$\displaystyle=$	$\displaystyle\omega_{h}^{2}[-\frac{1}{2}+(2\phi^{2}-1)\langle n_{c}^{eq}\rangle^{2}+\langle n_{h}^{ss}\rangle^{2}]$
		$\displaystyle+$	$\displaystyle\omega_{c}^{2}[-\frac{1}{2}+\langle n_{c}^{eq}\rangle^{2}+(2\phi^{2}-1)\langle n_{h}^{ss}\rangle^{2}]$
		$\displaystyle+$	$\displaystyle\omega_{h}\omega_{c}\phi(1-2\langle n_{c}^{eq}\rangle^{2}-2\langle n_{h}^{ss}\rangle^{2}),$

$\displaystyle\langle{q_{h}}\rangle=-i\frac{\partial\mathrm{ln}G(u,v)}{\partial v}\bigg{|}_{u=v=0}=\omega_{h}(\langle n_{h}^{ss}\rangle-\phi\langle n_{c}^{eq}\rangle),$

(20)

respectively, where we have used $\langle n_{c}^{eq}\rangle=\mathrm{coth}(\beta_{c}\omega_{c}/2)/2$ and $\langle n_{h}^{ss}\rangle={\mathrm{coth}(\beta_{h}^{\mathrm{eff}}\omega_{h}/2)}/{2}$ to denote the mean numbers of the cavity with the cold thermal reservoir and the hot nonthermal reservoir, respectively.

The average values of the work and the heat injection depend on the quantum discord, which enters into $\langle n_{h}^{ss}\rangle$, and under this condition, even the operation mode can change when the quantum discord is involved. Notably, the mean work given by Eq. (18) can be split into the following sum of two parts: $-\langle{w}\rangle=\langle{w_{\mathrm{adi}}}\rangle-\langle{w_{\mathrm{fric}}}\rangle$, where the first part $\langle w_{\mathrm{adi}}\rangle=(\omega_{h}-\omega_{c})(\langle n_{h}^{ss}\rangle-\langle n_{\mathrm{c}}^{eq}\rangle)$ is the mean work in the quantum adiabatic case, and the second part $\langle{w_{\mathrm{fric}}}\rangle=(\phi-1)(\omega_{h}\langle n_{c}^{eq}\rangle+\omega_{c}\langle n_{h}^{ss}\rangle)$ is the frictional work caused by diabatic transitions occurring along the two unitary driven processes.

As repeatedly emphasized in the work above, the quantum entanglement is vanishing in the case of multipartite separable states, where the quantum discord is, however, nonzero [see Fig. 1(b)]. The quantum discord is thus considered to be a more general quantum resource than the quantum concurrence in our engine model. To complete this picture, we have produced, as shown in Fig. 3(c), the average values of the work, the absorbed heat, the concurrence, and the discord as functions of the interaction strength $\xi$. We see that, for the given parameters, $\mathcal{Q}(\rho_{12})>0$ for $\xi>0$ but $C(\rho_{12})>0$ only for $\xi>2.922$; however, in the range where $C(\rho_{12})$ is finite, both $C(\rho_{12})$ and $\mathcal{Q}(\rho_{12})$ are monotonically increasing functions of $\xi$. Furthermore, the ranges of the quantities ($-\langle w\rangle$ and $\langle q_{c,h}\rangle)$ when parametrized using $\mathcal{Q}(\rho_{12})$ are much larger than those obtained when using $C(\rho_{12})$. Therefore, we use the quantum discord rather than the concurrence in the analysis of the quantum machine’s performance.

The thermodynamic efficiency, which acts as one of performance measures and is defined by $\eta_{th}=-\langle{w}\rangle/\langle q_{h}\rangle$, where $\langle{w}\rangle$ and $\langle{q_{h}}\rangle$ were given by Eqs. (18) and (20), respectively, can then be obtained as

$$\eta_{th}=1-\frac{\omega_{c}}{\omega_{h}}-\frac{\omega_{c}}{\omega_{h}}\frac{(\phi-1)(\langle n_{h}^{ss}\rangle+\langle n_{c}^{eq}\rangle)}{\langle n_{h}^{ss}\rangle-\phi\langle n_{c}^{eq}\rangle},$$

(21)

which reduces to the so-called Otto efficiency, $\eta_{O}=1-\omega_{\mathrm{c}}/\omega_{\mathrm{h}}$ if two unitary strokes are quantum adiabatic, irrespective of the existence of any correlations between the two atoms. The efficiency (21) can then be rewritten as $\eta_{th}=\eta_{C}^{\mathrm{gen}}-\langle\sigma\rangle/(\beta_{h}^{\mathrm{eff}}\langle q_{h}\rangle)$, where $\eta_{C}^{\mathrm{gen}}=1-\beta_{h}^{\mathrm{eff}}/\beta_{c}$ is the so-called generalized Carnot efficiency, and $\langle\sigma\rangle=-\beta_{h}^{\mathrm{eff}}\langle{q}_{h}\rangle-\beta_{c}\langle{q}_{c}\rangle\geq 0$ is the average entropy production of the heat engine during a single cycle ar22 . Although the efficiency $\eta_{\mathrm{th}}$ may surpass the Carnot efficiency $\eta_{C}=1-\beta_{h}/\beta_{c}$, it does satisfy the second law of thermodynamics because it must be bound by the generalized Carnot value $\eta_{C}^{\mathrm{gen}}$ as a result of positive entropy production.

Figure 4(a) shows that, for $\beta_{c}=0.6$ and $\beta_{h}=0.15$, as the quantum discord increases, both the mean heat injection $\langle{q_{h}}\rangle$ and the average work $-\langle w\rangle$ also increase, leading to higher thermodynamic efficiency $\eta_{th}$. Given the specifications selected for the control parameters, including $\omega_{c,h}$ and $\beta_{c,h}$, the machine always operates as a heat engine with $-\langle w\rangle>0$ and $\eta_{th}>0$ for finite or vanishing quantum discord. However, when different bath temperature values $\beta_{c,h}$ are selected, the machine mode may be changed by varying the discord. When the control parameters $\omega_{c,h}$ are given for appropriate values of the bath temperature, gradually increasing $\mathcal{Q}(\rho_{12})$ results in the occurrence of the following three consecutive operating modes [upper panels, Figs. 4(b) and 4(c) ]: (i) the refrigerator mode ($-\langle{w}\rangle<0,\langle{q_{h}}\rangle<0$, $\langle{q_{c}}\rangle>0$), (ii) the heater mode ($-\langle{w}\rangle<0,\langle{q_{c}}\rangle<0$), and (iii) the heat engine mode ($-\langle{w}\rangle>0$, $\langle{q_{h}}\rangle>0$, $\langle{q_{c}}\rangle<0$).

When the quantum discord $\mathcal{Q}(\rho_{12})$ increases gradually, the thermodynamic efficiency also increases and tends towards the Otto limit $\eta_{O}=1-\omega_{c}/\omega_{h}$, as shown in lower panels of Fig. 4(a)$-$4(c). For the heater mode, where the heat must flow into the cold reservoir, $\langle q_{c}\rangle<0$, the average heat injection $\langle q_{h}\rangle$ can be an arbitrary real number, and no useful work is extracted, i.e., $-\langle w\rangle\leq 0$. When operating as a refrigerator, the machine absorbs the heat from the cold reservoir, which entails $\langle q_{c}\rangle>0$, and releases this heat into the hot reservoir, such that $\langle q_{h}\rangle<0$, and work should be consumed, i.e., $-\langle w\rangle<0$. The refrigerator’s coefficient of performance is defined as the ratio of the cooling load to the work input, i.e., $\varepsilon=\langle q_{c}\rangle/\langle w\rangle$, which is positive [see the lower panels in Figs. 4(b) and 4(c)]. We also note from comparison of Fig. 4(b) with Fig. 4(c) that, even at a constant quantum discord, the machine mode can be changed when the bath temperature is modified.

The nonadiabatic factor (17) is dependent on the driving time $\tau_{\mathrm{dri}}$, and thus the efficiency $\eta_{th}$ is a function of both the quantum discord $\mathcal{Q}(\rho_{12})$ and the driving time $\tau_{\mathrm{dri}}$. When the driven stroke speeds downward, the nonadiabatic factor decreases, although not monotonically, and it approaches the lower bound of $\phi=1$ for large values of $\tau_{\mathrm{dri}}$ [see the inset of Fig. 5(a)]. Consequently, the efficiency increases very quickly at first and then levels off, becoming asymptotic to the Otto efficiency [see Fig. 5(a)]. In addition to the efficiency, the power output $P=-\langle{w}\rangle/\tau_{cyc}$ represents another important measure of the machine’s performance. The power output as a function of $\tau_{\mathrm{dri}}$ shows similar behavior to the efficiency [see Fig. 5(b)]. The average work output $-\langle w\rangle$ (18) is dependent on the driving time $\tau_{\mathrm{dri}}$ because of the time-dependent nonadiabatic factor $\phi=\phi(\tau_{\mathrm{dri}})$, and it reaches its maximum value $(-\langle w\rangle)_{\mathrm{max}}$ when $\phi=1$. Therefore, for different values of the quantum discord, the values of the dimensionless work $-\langle w\rangle/(-\langle w\rangle)_{\mathrm{max}}$ obtained as a function of $\tau_{\mathrm{dri}}$ fall on top of each other at the point $\phi(\tau_{\mathrm{dri}})=1$, as illustrated in the inset of Fig. 5(b). We also observe that, when using the given parameters $\beta_{h,c}$ and $\omega_{c,h}$, the driving time at the maximum work is independent of the interaction strength and thus remains fixed, but the driving time at the maximum power output changes if we change the interaction strength associated with the quantum discord. The power initially increases rapidly to approach the maximum value and then decreases nonmonotonically. In physical terms, when the driving time $\tau_{\mathrm{dri}}$ increases, the frictional work $\langle{w}_{\mathrm{fri}}\rangle$ drops rapidly, and then decreases very slowly to zero.

We also observe from Figs. 5(a) and 5(b) that both the power output and the efficiency are enhanced by the quantum discord. This enhancement behavior is caused by the dependences of the excitation number $\langle n_{h}^{ss}\rangle$ in Eqs. (18) and (21) on the quantum discord $\mathcal{Q}(\rho_{12})$. Physically, the quantum discord, as a type of quantum resource, contributes to the mean work extracted and thus enhances the machine’s performance by increasing its power output and efficiency. Furthermore, we note that the oscillations of the curves for the efficiency and the power diminish with increasing quantum discord. The final observation from Fig. 5 is that the quantum entanglement is vanishing at the fixed quantum discord $\mathcal{Q}(\rho_{12})=0.031$, and the engine is thus fueled by the reservoirs in separable states. The presence of the quantum discord enhances the machine’s performance by increasing its efficiency and power and also improves its stability by reducing the relative fluctuations in power (see also the inset of Fig. 6), even in the case where the reservoirs driving the quantum heat engines are in separable states, i.e., where the quantum entanglement is vanishing.

For quantum heat engines where the heat and the work are stochastic, the fluctuations in the power must be considered because they are associated with the machine’s stability. We consider the coefficient of variation of power $\sqrt{\delta P^{2}}/P$, which is equal to the square root of the relative work fluctuations, given by $-\sqrt{\delta{w}^{2}}/\langle{w}\rangle.$ This dimensionless coefficient is shown as a function of the quantum discord in Fig. 6, with the same coefficient being shown as a function of the driving time $\tau_{\mathrm{dri}}$ in the inset. The coefficient of variation of power decreases monotonically and its oscillations are damped as the quantum discord increases, thereby illustrating that the quantum discord can be used to damp the oscillations and the fluctuations in the power output, and thereby can improve the machine’s stability.

The total stochastic entropy production $\sigma$ is distributed according to the probability distribution $p(\sigma)$, where $\sigma(q_{h},w)=(\beta_{c}-\beta_{h}^{\mathrm{eff}})q_{h}+\beta_{c}w$. These distributions for the working system, which involve time-reversal symmetry and satisfy the fluctuation theorems, can be expressed as $p(\sigma)=p_{R}(-\sigma)e^{\sigma}$ Crooks98 , with $p_{R}(-\sigma)$ being the probability distribution of the time-reversed cycle [in the clockwise direction in Fig. 2(b)]. These theorems always imply the generalized thermodynamic uncertainty relationship for the stochastic work of the following form: $\delta w^{2}/\langle w\rangle^{2}\geq\mathrm{csch}^{2}[f(\langle\sigma\rangle)]$ Tim19 , where $f(x)$ is the inverse function of $x\tanh(x)$. Therefore, the coefficient of variation of power should satisfy the following relation: $\sqrt{\delta P^{2}}/P\geq\mathrm{csch}[f(\langle\sigma\rangle)],$ as shown in Fig. 6, where the function $\mathrm{csch}[f(\langle\sigma\rangle)]$ is indicated by the green dot-dot-dashed line.