Full statistics of nonequilibrium heat and work for many-body quantum Otto engines and universal bounds: A nonequilibrium Green’s function approach

Thermodynamic cycles are fundamental concepts in the field of thermodynamics that describe the exchange and conversion of energy and play a crucial role in understanding and optimizing the performance of various energy conversion devices Çengel and Boles (2004). Traditionally, thermodynamics provides a description of systems in terms of bulk or macroscopic properties by focusing primarily on averages or mean values of thermodynamic quantities Zemansky (1968); Callen et al. (1985). However, for small-scale systems consisting of a few degrees of freedom, thermal and quantum fluctuations about the mean values can be significant. Exploring the universal properties of these fluctuations for out-of-equilibrium systems has been of key interest in the field of stochastic and quantum thermodynamics Jarzynski (2011); Seifert (2012, 2008); Goold et al. (2016); Vinjanampathy and Anders (2016). Significant theoretical insight has been achieved in this regard over the past two decades, owing to the diverse kinds of universal fluctuation relations Jarzynski (2011); Holubec and Ryabov (2021); Esposito et al. (2009); Campisi et al. (2011).

In recent times, quantum engines have garnered significant interest as they provide a privileged platform to investigate the intricate relationship between quantum phenomena and nonequilibrium thermodynamics. Furthermore, they offer opportunities to harness the unique properties of the quantum working medium to boost performance and enhance stability Jaramillo et al. (2016); Vroylandt et al. (2018); Niedenzu and Kurizki (2018); Boubakour et al. (2023). Thanks to rapid technological developments, it has now become possible to fabricate small-scale thermal machines Peterson et al. (2019); Myers et al. (2022a); Blickle and Bechinger (2012); Martínez et al. (2016, 2017); Bouton et al. (2021); Maslennikov et al. (2019); Abah et al. (2012a); Roßnagel et al. (2016); Zhang et al. (2022); Bhattacharjee and Dutta (2021). For example, the smallest possible quantum Otto engine, with a qubit as working fluid has been realized in an NMR setup Peterson et al. (2019). Another paradigmatic model, a harmonic oscillator heat engine has also been realized in the ion-trap setup Roßnagel et al. (2016). More complex and exotic many-body working fluids executing a cycle has now been proposed Carollo et al. (2020); Chen et al. (2019); Sun et al. (2017); Altintas et al. (2014); Thomas and Johal (2011); Huang et al. (2013a); Myers et al. (2022b); Henrich et al. (2007a); Huang et al. (2013b); Singh and Rebari (2020); Jaramillo et al. (2016).

The study of fluctuation in the context of such small-scale thermal machines is of tremendous significance Watanabe and Minami (2022); Ito et al. (2019); Denzler and Lutz (2020); Pietzonka and Seifert (2018a); Jiao et al. (2021). Recently, there have been intense ongoing efforts to quantify the performance of quantum engines involving efficiency, power, and power fluctuation or so-called constancy Pietzonka and Seifert (2018b). In this regard, Thermodynamic uncertainty relations (TURs), establishing a trade-off relation between the relative fluctuation of thermodynamic currents and entropy production, have been put forward that rule out the possibility of achieving Carnot efficiency without allowing power fluctuation to diverge Barato and Seifert (2015); Pietzonka and Seifert (2018b); Gingrich et al. (2016); Falasco et al. (2020); Guarnieri et al. (2019); Brandner et al. (2018); Agarwalla and Segal (2018). Another direction of studies established universal upper and lower bounds on the ratio of power to input heat current fluctuations of generic continuously coupled autonomous thermal machines operating in the linear response regime Saryal et al. (2021). These studies have been further extended to discrete four-stroke Otto cycles with specific working mediums Saryal and Agarwalla (2021); Mohanta et al. (2023) and very recently for model-specific periodically driven continuous machines Das et al. (2023). However, a detailed understanding of the bounds on nonequilibrium fluctuations for non-autonomous machines with generic many-body working medium is still lacking.

In this work, we consider a generic four-stroke quantum Otto cycle Solfanelli et al. (2020); Kosloff (2019); Geva and Kosloff (1992); Feldmann and Kosloff (2000); Kieu (2004); Henrich et al. (2007b); Scully (2002); Watanabe et al. (2017); Rezek and Kosloff (2006); Lin and Chen (2003); Abah et al. (2012b); Kosloff and Rezek (2017); Uzdin et al. (2015) consisting of an arbitrary many-body working medium and derive an analytical expression for the joint cumulant generating function (CGF) of total work and input heat, valid up to the second order of the driving amplitude, by employing the rigorous Schwinger-Keldysh nonequilibrium Green’s function (NEGF) formalism Fei and Quan (2020); Cavina et al. (2023); Rammer and Smith (1986); Schwinger (1961). The CGF obtained for an arbitrary many-body working medium satisfies the nonequilibrium fluctuation relation Campisi (2014). Furthermore, we demonstrate the linear response limit characterized by a small driving amplitude and small temperature difference and obtain the expressions for the Onsager’s transport coefficients. Notably, we reveal a breakdown of the traditional fluctuation-dissipation relation (FDR) Kubo (1966); Weber (1956); Felderhof (1978); Pathria and Beale (2021) for the total work, while, the traditional FDR remains intact for heat. Our findings have remarkable implications for the generic Otto cycle: when operating as an engine, the ratio of output (work) to input (heat from the hot reservoir) fluctuations is universally bounded from below, while in the refrigerator regime, a similar quantity is universally bounded from above. Additionally, we also establish connections to the TURs. Finally, we compare the obtained universal bounds for the discrete Otto cycle with those of autonomous steady-state machines Saryal et al. (2021); Mohanta et al. (2023); Saryal and Agarwalla (2021); Das et al. (2023).

We organize the paper as follows: Section II briefly describes the quantum Otto cycle and the two-point measurement protocol to construct the CGF corresponding to the joint probability distribution of total work and heat. We then demonstrate how to map the CGF on the modified Keldysh contour. In Section III, we obtain the expression for the CGF valid up to the second order of the external driving amplitude. We validate that the CGF respects the fluctuation symmetry. In Section IV, we discuss the linear response limit, obtain the Onsager transport coefficients, and derive universal bounds on nonequilibrium fluctuations. We illustrate the obtained results with a paradigmatic model in Section V. Finally, we summarize our main findings and provide an outlook in Section VI.

We consider a generic many-body working medium for the quantum Otto cycle and treat the many-body interaction part as a time-dependent perturbation. The total Hamiltonian of the working fluid is given by

$$H(t)=H_{0}+\lambda(t)H_{1},$$

(1)

where $H_{0}$ is the non-interacting part, often referred to here as the free Hamiltonian, and $H_{1}$ represents an arbitrary many-body interaction term. $\lambda(t)$ is an arbitrary driving protocol that drives the system away from equilibrium. Importantly, $\lambda(t)$ serves the purpose of a perturbation parameter, and the values it can take are significantly smaller than those of other dimensionally comparable free parameters.

The working medium performs a four-stroke Otto cycle consisting of two unitary and two thermalization strokes with respect to the cold and hot thermal reservoirs with fixed inverse temperatures $\beta_{c}$ and $\beta_{h}$, respectively (see Fig. 1). Below we briefly describe these four strokes: (1) Unitary expansion stroke ($A\!\to\!B$). The system starts at time $t\!=\!0$ from the initial thermal state $A$ characterized by the cold inverse temperature $\beta_{c}$ and Hamiltonian $H(0)\!=\!H_{0}+\lambda(0)H_{1}$, where $\lambda(0)\!=\!\lambda_{0}$. The system is then decoupled from the reservoir and $\lambda(t)$ drives the system out of equilibrium. The first unitary stroke ends at time $t\!=\!\tau_{1}$ with $\lambda(\tau_{1})\!=\!\lambda_{1}$. At the end of this stroke, the state of the system $B$ is a non-thermal one. (2) Hot isochore ($B\!\to\!C$). During the first heat exchange stroke, the system is brought into weak contact with the hot reservoir while keeping $\lambda_{1}$ fixed. The system then waits until it reaches a thermal state $C$, parameterized by ($\lambda_{1},\beta_{h}$) at time $t\!=\!\tau_{1}+\tau_{h}$. (3) Unitary compression ($C\!\to\!D$). The system is then decoupled from the hot reservoir and the interaction strength $\lambda(t)$ is decreased unitarily from $\lambda_{1}$ to $\lambda_{0}$, and consequently, the system ends up in a non-thermal state $D$ at time $t\!=\!\tau_{1}+\tau_{h}+\tau_{2}$. (4) Cold isochore ($D\!\to\!A$): In the final step, the system is placed in weak contact with the cold reservoir while keeping $\lambda_{0}$ fixed until it reaches the thermal state parametrized by ($\lambda_{0},\beta_{c}$) at time $t=\tau_{1}+\tau_{h}+\tau_{2}+\tau_{c}=\tau_{cyc}$. Here, we assume that the system achieves full thermalization during both hot and cold isochores. This assumption is valid when the heat exchange stroke times $\tau_{h}$ and $\tau_{c}$ are much larger than the system relaxation time. Additionally, we only consider here the symmetric driving case, meaning the unitary evolution operators, represented by $\mathcal{U}_{\mathrm{e}}$ and $\mathcal{U}_{\mathrm{c}}$, governing the expansion and compression strokes, respectively, are mirror images of each other, i.e., $\mathcal{U}_{\mathrm{c}}\!=\!\Theta\mathcal{U}_{\mathrm{e}}^{\dagger}\Theta^{\dagger}$, where $\Theta$ is the anti-unitary time-reversal operator, and as a result, the two unitary stroke duration are equal, $\tau_{1}\!=\!\tau_{2}$, which we set as $\tau$.

Under the complete thermalization assumption, we construct the joint probability distribution (PD) of total work and heat exchange in the hot isochore $P(w,q_{h})$ by employing the projective measurements of the energy of the system at the end points of each stroke i.e., at $A$, $B$, $C$, and $D$, respectively Denzler and Lutz (2020). The characteristic function (CF) corresponding to the joint distribution is given by the Fourier transform of the PD,

	$\displaystyle\chi($	$\displaystyle\gamma_{w},\gamma_{h})=\int\!\!\!\int d{w}\,d{q_{h}}\,P({w},{q_{h}})\,e^{i\gamma_{w}{w}}\,e^{i\gamma_{h}{q_{h}}}$
		$\displaystyle=\frac{\mathrm{Tr}\big{[}{\cal U}_{\mathrm{e}}^{\dagger}\,e^{(i\gamma_{w}-i\gamma_{h})H(\tau)}\,{\cal U}_{\mathrm{e}}\,e^{(-i\gamma_{w}-\beta_{c})H(0)}\big{]}}{{\cal Z}_{c}[H(0)]}$
		$\displaystyle\times\frac{\mathrm{Tr}\big{[}{\cal U}_{\mathrm{c}}^{\dagger}\,e^{i\gamma_{w}H(0)}\,{\cal U}_{\mathrm{c}}\,e^{(-i\gamma_{w}+i\gamma_{h}-\beta_{h})H(\tau)}\big{]}}{{\cal Z}_{h}[H(\tau)]},$		(2)

where $\gamma_{w},\gamma_{h}$ are the counting variables associated to the total work $w\!=\!w_{1}+w_{3}$ and the heat exchange with the hot reservoir $q_{h}$, respectively. In Eq. (2), $H(0)\!=\!H_{0}+\lambda_{0}H_{1}$ and $H(\tau)\!=\!H_{0}+\lambda_{1}H_{1}$. ${\cal Z}_{c}[H(0)]\!=\!\operatorname{Tr}[e^{-\beta_{c}H(0)}]$ and ${\cal Z}_{h}[H(\tau)]\!=\!\operatorname{Tr}[e^{-\beta_{h}H(\tau)}]$ are the canonical partition functions with respect to the cold and hot inverse temperatures $\beta_{c}$ and $\beta_{h}$, and Hamiltonians $H(0)$ and $H(\tau)$, respectively. We observe that, due to the complete thermalization assumption, the joint PD in Eq. (2) takes a product form. Furthermore, an intriguing implication of this assumption is that the joint PD of total work and heat from the hot reservoir alone suffices to determine the statistics of heat exchange with the cold reservoir Mohanta et al. (2023).

In order to obtain a perturbative expression of the CF given in Eq. (2), we utilize the Schwinger-Keldysh formal- ism, a systematic approach for calculating nonequilibrium Green’s functions and correlation functions under the influence of time-dependent perturbations by casting the problem on a closed time-ordered contour in the complex time plane, known as the Keldysh contour. For the problem considered here, all the unitary evolution operators and the exponential operators, e.g., $e^{(i\gamma_{w}-i\gamma_{h})H(\tau)}$, appearing in Eq. (2) are mapped on two modified Keldysh contours. The first trace in Eq. (2) corresponds to the modified contour $C$, as shown in Fig. 2(a), and the second trace gives rise to another modified contour $C^{\prime}$, as shown in Fig. 2(b). It is important to mention that, in the Schwinger-Keldysh formalism one typically expects closed contours while dealing with time-dependent perturbations. However, here, none of our modified contours are closed, only articulating the fact that both the traces that are mapped contain partial information about the non-unitary heat exchange stroke. The tails represent the presence of interaction $H_{1}$ in the thermal states at $A$ and $C$ (see Fig. 1). We have demonstrated the details of this mapping in Appendix A and here we only provide a compact expression for the joint CF on the Keldysh contour,

	$\displaystyle\chi(\gamma_{w},\gamma_{h})$	$\displaystyle=\frac{\big{\langle}{\cal T}_{C}\,e^{-\frac{i}{\hbar}\int_{C}\!\lambda_{C}(s_{1})H_{1}^{I}(s_{1})ds_{1}}\big{\rangle}_{\widetilde{\beta}_{c}}}{\big{\langle}{\cal T}_{C}\,e^{-\frac{i}{\hbar}\int_{\hbar\gamma_{h}}^{\hbar\gamma_{h}\!-\!i\hbar\beta_{c}}\!\lambda_{0}H_{1}^{I}(s_{1})ds_{1}}\big{\rangle}_{{\beta}_{c}}}$
		$\displaystyle\times\,\frac{\big{\langle}{\cal T}_{C^{\prime}}\,e^{-\frac{i}{\hbar}\int_{C}\!{\lambda}_{C^{\prime}}(s_{2})H_{1}^{I}(s_{2})ds_{2}}\big{\rangle}_{\widetilde{\beta}_{h}}}{\big{\langle}{\cal T}_{C^{\prime}}\,e^{-\frac{i}{\hbar}\int_{\tau\!-\!\hbar\gamma_{h}}^{\tau\!-\!\hbar\gamma_{h}\!-\!i\hbar\beta_{h}}\!\lambda_{1}H_{1}^{I}(s_{2})ds_{2}}\big{\rangle}_{{\beta}_{h}}}$
		$\displaystyle\times\,\langle e^{-i\gamma_{h}H_{0}}\rangle_{\beta_{c}}\langle e^{i\gamma_{h}H_{0}}\rangle_{\beta_{h}}.$		(3)

Here, crucially, the averages in both the numerators are taken with respect to shifted inverse temperatures $\widetilde{\beta}_{c}\!=\!\beta_{c}+i\gamma_{h}$ and $\widetilde{\beta}_{h}\!=\!{\beta}_{h}-i\gamma_{h}$ which are counting-field dependent. ${\cal T}_{C}$ and ${\cal T}_{C^{\prime}}$ are the contour-time ordered operators defined on the contours $C$ and $C^{\prime}$, respectively, which orders the time-dependent operators according to their contour-time argument; i.e., sorting the operators from left to right with their contour-time arguments decreasing. In Appendix B we have expanded the exponential operators in Eq. (3) up to the second order in driving protocol $\lambda(s)$ and obtained the expression of the joint cumulant generating function (CGF) $\ln{\chi(\gamma_{w},\gamma_{h})}$ for arbitrary many-body working medium. The CGF is a useful quantity as it enables us to calculate all the cumulants (denoted by double angular brackets) of total work and heat exchange with the hot reservoir by taking partial derivatives with respect to $\gamma_{w}$ and $\gamma_{h}$, respectively,

$\displaystyle\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}w^{k}q_{h}^{l}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}=\frac{\partial^{k}\partial^{l}}{\partial(i\gamma_{w})^{k}\partial(i\gamma_{h})^{l}}\,\ln{\chi(\gamma_{w},\gamma_{h})}\Big{|}_{\gamma_{w},\gamma_{h}=0}.$

(4)

In this section, we write down the final expression for the CGF valid up to the second order of the driving protocol $\lambda(s)$ (see Appendix B for the details). We obtain

	$\displaystyle\ln\chi($	$\displaystyle\gamma_{w},\gamma_{h})\!=\!\big{[}\!-\!i\gamma_{w}(\lambda_{1}\!-\!\lambda_{0})+i\gamma_{h}\lambda_{1}\big{]}\,\big{[}\langle H_{1}\rangle_{\widetilde{\beta}_{h}}\!\!-\!\langle H_{1}\rangle_{\widetilde{\beta}_{c}}\big{]}+\big{[}\!-\!i\hbar\gamma_{w}(\lambda_{1}^{2}\!-\!\lambda_{0}^{2})+i\hbar\gamma_{h}\lambda_{1}^{2}\big{]}\!\int\frac{d\omega}{2\pi}\frac{\widetilde{G}_{h}^{>}(\omega)\!-\!\widetilde{G}_{c}^{>}(\omega)}{\omega}$
		$\displaystyle-\!\beta_{c}\lambda_{0}\Big{[}\langle H_{1}\rangle_{\widetilde{\beta}_{c}}\!\!-\!\langle H_{1}\rangle_{\beta_{c}}+\hbar\lambda_{0}\!\int\frac{d\omega}{2\pi}\frac{\widetilde{G}_{c}^{>}(\omega)\!-\!{G}_{c}^{>}(\omega)}{\omega}\Big{]}\!-\!\beta_{h}\lambda_{1}\Big{[}\langle H_{1}\rangle_{\widetilde{\beta}_{h}}\!\!-\!\langle H_{1}\rangle_{\beta_{h}}+\hbar\lambda_{1}\!\int\frac{d\omega}{2\pi}\frac{\widetilde{G}_{h}^{>}(\omega)\!-\!{G}_{h}^{>}(\omega)}{\omega}\Big{]}$
		$\displaystyle+\!\int\frac{d\omega}{2\pi}\frac{1\!-\!e^{i\hbar\omega(\gamma_{w}\!-\!\gamma_{h})}}{\omega^{2}}A(\omega)\widetilde{G}_{c}^{>}(\omega)+\!\int\frac{d\omega}{2\pi}\frac{1\!-\!e^{i\hbar\omega\gamma_{w}}}{\omega^{2}}{A}(\omega)\widetilde{G}_{h}^{>}(\omega)+\ln{\langle e^{-i\gamma_{h}H_{0}}\rangle_{\beta_{c}}}+\ln{\langle e^{i\gamma_{h}H_{0}}\rangle_{\beta_{h}}}.$		(5)

Equation (5) is the first central result of this paper. It describes the statistical properties of total work and input heat fluctuations in the many-body Otto cycle in terms of one-point and two-point correlators of the interaction Hamiltonian $H_{1}$. In this expression, we suppress the integration limits for $\omega$, which ranges from $-\infty$ to $+\infty$. Here, the nonequilibrium Green’s functions are defined as connected two-point correlators of $H_{1}$ in the interaction picture,

	$\displaystyle\widetilde{G}_{j}^{>}(s)$	$\displaystyle=\left(\!-\frac{i}{\hbar}\right)^{2}\!\!\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{1}^{I}(s)H_{1}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\widetilde{\beta}_{j}},$		(6)
	$\displaystyle{G}_{j}^{>}(s)$	$\displaystyle=\left(\!-\frac{i}{\hbar}\right)^{2}\!\!\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{1}^{I}(s)H_{1}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{{\beta}_{j}},$		(7)

where $j\!=\!c,h$. In this context, $\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\cdot\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$ implies a connected correlation function, also known as a cumulant correlation function,

$$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}X^{I}(s_{1})Y^{I}(s_{2})\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\beta}=\langle X^{I}(s_{1})Y^{I}(s_{2})\rangle_{\beta}\!-\!\langle X\rangle_{\beta}\langle Y\rangle_{\beta}.$$

(8)

The shifted inverse temperatures $\widetilde{\beta}_{j}$, $j\!=\!c,h$, appearing in Eq. (5) play a crucial role in calculating the cumulants of the heat exchange. In this paper, we adopt the Fourier transformation convention as

$${G}_{j}^{>}(\omega)=\!\int_{-\infty}^{+\infty}\!\!ds\,{G}_{j}^{>}(s)\,e^{i\omega s}.$$

(9)

Crucially, all the information about the time-dependent driving protocol is encapsulated in the quantity $A(\omega)$, defined as

$\displaystyle A(\omega)=\left|\int_{0}^{\tau}\!\dot{\lambda}(t)e^{i\omega t}dt\right|^{2}.$

(10)

It is worth nothing that the CGF in Eq. (5) satisfies the normalization condition, $\ln{\chi(0,0)}\!=\!0$. More importantly, it satisfies the universal heat engine fluctuation symmetry Campisi (2014); Campisi et al. (2015); Jarzynski and Wójcik (2004),

$$\ln\chi\big{[}i\beta_{c},\,i(\beta_{c}\!-\!\beta_{h})\big{]}=0.$$

(11)

Notably, this fluctuation symmetry is a direct consequence of the Kubo-Martin-Schwinger (KMS) boundary condition satisfied by the Green’s functions introduced in Eqs. (6) and (7) Breuer and Petruccione (2002); Carmichael (1993); Kubo (1957); Martin and Schwinger (1959),

	$\displaystyle\widetilde{G}^{>}(s\!-\!i\hbar\widetilde{\beta}_{c})=\widetilde{G}^{>}(-s),$		(12)
	$\displaystyle{G}^{>}(s\!-\!i\hbar{\beta}_{c})={G}^{>}(-s).$		(13)

Here, $\widetilde{\beta}_{c}$ and ${\beta}_{c}$ are the shifted and original inverse temperatures of the cold reservoir, respectively. Furthermore, the fluctuation symmetry in Eq. (11) implies that $\langle e^{-\Sigma}\rangle\!=\!1$, where $\Sigma\!=\!\beta_{c}w+(\beta_{c}\!-\!\beta_{h})q_{h}$ is the stochastic entropy production in a single cycle. It is worth highlighting that utilizing Jensen’s inequality, $\langle e^{-\Sigma}\rangle\geq e^{-\langle\Sigma\rangle}$, we immediately recover the second law of thermodynamics, $\langle\Sigma\rangle\geq 0$ Landi and Paternostro (2021).

In the following, we present the formal expressions of the first and second cumulants of total work and heat exchange with the hot reservoir. The average total work is given by

	$\displaystyle\langle w\rangle\!=$	$\displaystyle-\!(\lambda_{1}\!-\!\lambda_{0})\left[\langle H_{1}\rangle_{\beta_{h}}\!\!-\!\langle H_{1}\rangle_{\beta_{c}}\right]$
		$\displaystyle-\hbar(\lambda_{1}^{2}\!-\!\lambda_{0}^{2})\!\int\!\frac{d\omega}{2\pi}\frac{1}{\omega}\left[{G}_{h}^{>}(\omega)\!-\!{G}_{c}^{>}(\omega)\right]$
		$\displaystyle-\!\hbar\!\int\!\frac{d\omega}{2\pi}\frac{A(\omega)}{\omega}\left[G^{>}_{c}(\omega)+G^{>}_{h}(\omega)\right].$		(14)

The first two terms on the right-hand side represent the boundary terms that depend on the end values of the driving protocol and are associated with the total change in free energy during the two unitary work strokes.

The corresponding variance (second cumulant) of total work is given by

$$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}w^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}=-\hbar^{2}\!\int\!\frac{d\omega}{2\pi}A(\omega)\left[G^{>}_{c}(\omega)+G^{>}_{h}(\omega)\right].$$

(15)

Notably, we find that calculating the cumulants of heat exchange $q_{h}$ requires more involved calculations due to the presence of shifted inverse temperatures $\widetilde{\beta}_{c}$ and $\widetilde{\beta}_{h}$. This results in connected correlators of the form $\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{1}H_{0}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$, $\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{1}^{I}(s)H_{1}H_{0}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$, and so on. We note that, for an arbitrary operator $X$,

	$$\displaystyle\frac{d}{d(i\gamma_{h})}\langle X\rangle_{\widetilde{\beta}_{j}}\!=\!\mp\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}XH_{0}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\widetilde{\beta}_{j}},$$		(16)
	$$\displaystyle\frac{d^{2}}{d(i\gamma_{h})^{2}}\langle X\rangle_{\widetilde{\beta}_{j}}\!=\!\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}XH_{0}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\widetilde{\beta}_{j}},$$		(17)

where $(\!-\!)$ corresponds to $j\!=\!c$, and $(+)$ corresponds to $j\!=\!h$, respectively. Below we list the expressions for the average heat exchange with the hot reservoir and its variance,

	$\displaystyle\langle q_{h}\rangle\!=$	$\displaystyle\langle H_{0}\rangle_{\beta_{h}}\!\!-\!\langle H_{0}\rangle_{\beta_{c}}+\lambda_{1}\big{[}\langle H_{1}\rangle_{\beta_{h}}\!\!-\!\langle H_{1}\rangle_{\beta_{c}}\big{]}+\hbar\lambda_{1}^{2}\!\int\!{\frac{d\omega}{2\pi}\frac{{G}_{h}^{>}(\omega)\!-\!{G}_{c}^{>}(\omega)}{\omega}}+\beta_{c}\lambda_{0}\Big{[}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{0}H_{1}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\beta_{c}}\!+\hbar\lambda_{0}\!\int\!\frac{d\omega}{2\pi}\frac{G^{{}^{\prime}}_{c}(\omega)}{\omega}\Big{]}$
		$\displaystyle\quad\quad\quad\quad-\!\beta_{h}\lambda_{1}\Big{[}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{0}H_{1}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\beta_{h}}\!+\hbar\lambda_{1}\!\int\!\frac{d\omega}{2\pi}\frac{G^{{}^{\prime}}_{h}(\omega)}{\omega}\Big{]}+\hbar\!\int\!\frac{d\omega}{2\pi}\frac{A(\omega)G^{>}_{c}(\omega)}{\omega},$		(18)
	$\displaystyle\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}q_{h}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}\!=$	$\displaystyle\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{0}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\beta_{c}}\!+\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{0}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\beta_{h}}\!+2\lambda_{1}\big{[}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{0}H_{1}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\beta_{c}}\!+\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{0}H_{1}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\beta_{h}}\big{]}+2\hbar\lambda_{1}^{2}\int\!\frac{d\omega}{\pi}\frac{G_{c}^{{}^{\prime}}(\omega)+G_{h}^{{}^{\prime}}(\omega)}{\omega}\!-\!\beta_{c}\lambda_{0}\Big{[}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{0}^{2}H_{1}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\beta_{c}}$
		$\displaystyle\!\!+\hbar\lambda_{0}\!\int\!\frac{d\omega}{2\pi}\frac{G_{c}^{{}^{\prime\prime}}(\omega)}{\omega}\Big{]}\!-\!\beta_{h}\lambda_{1}\Big{[}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}H_{0}^{2}H_{1}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\beta_{h}}\!+\hbar\lambda_{1}\!\int\!\frac{d\omega}{2\pi}\frac{G_{h}^{{}^{\prime\prime}}(\omega)}{\omega}\Big{]}+2\hbar\!\int\!\frac{d\omega}{2\pi}\frac{G^{{}^{\prime}}_{c}(\omega)A(\omega)}{\omega}+\hbar^{2}\!\int\!\frac{d\omega}{2\pi}{A(\omega)G_{c}^{>}(\omega)},$		(19)

where we have introduced a compact notation,

	$\displaystyle G_{j}^{{}^{\prime}}(\omega)$	$\displaystyle=$	$\displaystyle-\frac{\partial}{\partial\beta_{j}}G_{j}^{>}(\omega),$		(20)
	$\displaystyle G_{j}^{{}^{\prime\prime}}(\omega)$	$\displaystyle=$	$\displaystyle\frac{\partial^{2}}{\partial\beta_{j}^{2}}G_{j}^{>}(\omega),$		(21)

with $j=c,h$.

Now that we have obtained the expressions for the mean and the variance for total work and heat exchange with the hot reservoir from the joint CGF, in the following section, we focus on establishing the linear response regime for the generic quantum Otto cycle. We derive the expressions for Onsager’s transport coefficients and explore the fate of traditional fluctuation-dissipation relations for work and heat.

Let us first parameterize the driving protocol $\lambda(t)$ as

$\displaystyle\lambda(t)=\lambda+\Delta\lambda\,f(t),$

(22)

where, $\lambda\!=\!\lambda_{0}$ and $\Delta\lambda\!=\!\lambda_{1}\!-\!\lambda_{0}$. The function $f(t)$ is defined such that $f(0)\!=\!0$ and $f(\tau)\!=\!1$ to ensure that $\lambda_{0}$ and $\lambda_{1}$ are the boundary values of the protocol $\lambda(t)$ at time $t\!=\!0$ and time $t\!=\!\tau$, respectively. Recall that $\tau$ is the duration of each work stroke. To identify the proper thermodynamic affinities and the corresponding conjugate fluxes, we begin with the definition of the stochastic entropy production in a single cycle:

	$\displaystyle\Sigma=$	$\displaystyle\beta_{c}w+(\beta_{c}\!-\!\beta_{h})\,q_{h}$
	$\displaystyle=$	$\displaystyle\bm{(}\beta_{c}\Delta\lambda\bm{)}\frac{w}{\Delta\lambda}+\bm{(}\beta_{c}\!-\!\beta_{h}\bm{)}q_{h}.$		(23)

From this expression, we can identify the thermodynamic affinities as ${\cal F}_{w}\!=\!\beta_{c}\Delta\lambda$ and ${\cal F}_{q}\!=\!\beta_{c}\!-\!\beta_{h}$, and the corresponding integrated work and heat fluxes as $w/\Delta\lambda$ and $q_{h}$, respectively. One can further define the heat and work currents by dividing the integrated heat and work fluxes with total cycle time $\tau_{cyc}\!=\!2\tau+\tau_{h}+\tau_{c}$, as follows:

$$\displaystyle\;\;\mathcal{J}_{w}\!=\!\frac{1}{\tau_{cyc}}\frac{w}{\Delta\lambda}\;\;\mathrm{and}\;\;\mathcal{J}_{h}\!=\!\frac{q_{h}}{\tau_{cyc}}.$$

(24)

Consequently, the average entropy production rate in one complete cycle $\langle\sigma\rangle\!=\!\langle\Sigma\rangle/\tau_{cyc}$ can be expressed as

$\displaystyle\langle\sigma\rangle=\sum_{i}\langle\mathcal{J}_{i}\rangle\mathcal{F}_{i},$

(25)

where $i\!=\!w,h$. From this point onwards, we adopt $\beta_{c}\!=\!\beta$ and $\beta_{c}\!-\!\beta_{h}\!=\!\Delta\beta$, and focus on the linear response regime where $\Delta\beta\ll\beta$. In this regime, we can effectively utilize the following expansion,

	$\displaystyle\langle X\rangle_{\beta\!-\!\Delta\beta}=$	$\displaystyle\langle X\rangle_{\beta}\!-\!\Delta\beta\,\frac{\partial}{\partial\beta}\langle X\rangle_{\beta}$
	$\displaystyle=$	$\displaystyle\langle X\rangle_{\beta}+\Delta\beta\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}XH_{0}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\beta}.$		(26)

Note that the appearance of the free Hamiltonian $H_{0}$ in the second term arises from the derivative with respect to $\beta$ over a canonical distribution characterized by $e^{-\beta H_{0}}/\operatorname{Tr}{[e^{-\beta H_{0}}]}$. To proceed with our linear response analysis for the generic Otto cycle, we calculate the average work and heat fluxes, $\langle\mathcal{J}_{w}\rangle$ and $\langle\mathcal{J}_{h}\rangle$, using Eqs. (14) and (18), respectively, which are accurate up to the second order of the driving protocol $\lambda(t)$. Finally, carrying out the Taylor expansion given in Eq. (26) and keeping terms up to the linear order of the affinities, $\mathcal{F}_{w}$ and $\mathcal{F}_{q}$, we derive the following relation between the currents and affinities, connected via the Onsager matrix:

$\displaystyle\begin{pmatrix}\langle\mathcal{J}_{w}\rangle\\ \langle\mathcal{J}_{h}\rangle\end{pmatrix}=\begin{pmatrix}{\cal L}_{ww}&{\cal L}_{wq}\\ {\cal L}_{qw}&{\cal L}_{qq}\end{pmatrix}\begin{pmatrix}{\cal F}_{w}\\ {\cal F}_{q}\end{pmatrix},$

(27)

where the Onsager transport coefficients ${\cal L}_{ij}$, with $i,j=w,q$, are given by

	$\displaystyle{\cal L}_{ww}\!$	$\displaystyle=\!-\frac{1}{\tau_{cyc}}\hbar\!\int\!\frac{d\omega}{2\pi}\frac{2F(\omega)G^{>}(\omega)}{\beta\,\omega},$		(28)
	$\displaystyle{\cal L}_{wq}\!$	$\displaystyle=\!{\cal L}_{qw}\!=\!\frac{1}{\tau_{cyc}}\frac{\partial}{\partial\beta}\Big{[}\langle H_{1}\rangle_{\beta}\!+\!2\hbar\lambda\!\!\int\!\frac{d\omega}{2\pi}\frac{G^{>}(\omega)}{\omega}\Big{]},$		(29)
	$\displaystyle\mathcal{L}_{qq}\!$	$\displaystyle=\!-\frac{\partial}{\partial\beta}\Big{[}\langle H_{0}\rangle_{\beta}\!+\!\lambda\frac{\partial}{\partial\beta}\Big{(}\beta\langle H_{1}\rangle_{\beta}\!+\hbar\lambda\beta\!\!\int\!\frac{d\omega}{2\pi}\frac{G^{>}(\omega)}{\omega}\Big{)}\Big{]}.$		(30)

Here, note that

$\displaystyle F(\omega)=\left|\int_{0}^{\tau}\!\dot{f}(t)e^{i\omega t}dt\right|^{2}\!=\frac{A(\omega)}{{\Delta\lambda}^{2}}$

(31)

encodes the information of driving where $\Delta\lambda$ and $f(t)$ are defined in Eq. (22). The Green’s function $G^{>}(\omega)$ is evaluated at the equilibrium temperature $\beta$. Importantly, we recover the Onsager reciprocity relation Onsager (1931), i.e. $\mathcal{L}_{wq}\!=\!\mathcal{L}_{qw}$, as shown in Eq. (29). The derived expressions for the Onsager transport coefficients, characterizing the behavior of an arbitrary working medium executing the Otto cycle in the linear response regime, constitute another central result of our work. With these results in hand, we proceed to provide some important insights and remarks.

The traditional “Fluctuation-Dissipation relations” (FDRs) typically connect the diagonal elements of the Onsager matrix to equilibrium fluctuations Kubo (1966). Remarkably, in our analysis, we precisely recover the traditional FDR for the heat current $\langle\mathcal{J}_{h}\rangle$ (detailed derivation in Appendix E),

$\displaystyle{{\cal L}_{qq}}=\frac{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{h}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{2}\bigg{|}_{\mathcal{F}_{w},\mathcal{F}_{q}=0},$

(32)

where we recall that $\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{h}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}\!=\!\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}q_{h}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}/\tau_{cyc}$. However, importantly, the same is not true for the case of work current. In fact, for the work current fluctuation, following Eq. (15), we obtain

$\displaystyle\frac{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{w}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{2}\bigg{|}_{\mathcal{F}_{w},\mathcal{F}_{q}=0}=-\frac{1}{\tau_{cyc}}\hbar^{2}\!\int\!\frac{d\omega}{2\pi}F(\omega)G^{>}(\omega),$

(33)

which does not match the diagonal Onsager coefficient $\mathcal{L}_{ww}$. To investigate this discrepancy further, let us rewrite the expressions for $\mathcal{L}_{ww}$ and $\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{w}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$ as follows:

	$$\displaystyle{\cal L}_{ww}\!=\!-\frac{1}{\tau_{cyc}}\frac{\hbar^{2}}{2}\!\!\int\!\frac{d\omega}{2\pi}F(\omega)\frac{\tanh{\big{(}\frac{\beta\hbar\omega}{2}\big{)}}}{\frac{\beta\hbar\omega}{2}}\big{[}G^{>}(\omega)+G^{>}(\!-\omega)\big{]}.$$		(34)
	$$\displaystyle\frac{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{w}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{2}\bigg{|}_{\mathcal{F}_{w},\mathcal{F}_{q}=0}\!=\!\!-\frac{1}{\tau_{cyc}}\frac{\hbar^{2}}{2}\!\int\!\frac{d\omega}{2\pi}F(\omega)\big{[}G^{>}(\omega)+G^{>}(\!-\omega)\big{]},$$		(35)

Here, we have utilized an important identity involving the Green’s functions,

$$G^{>}(\omega)\!-\!G^{>}(\!-\omega)=\tanh{\Big{(}\frac{\beta\hbar\omega}{2}\Big{)}}\big{[}G^{>}(\omega)+G^{>}(\!-\omega)\big{]}.$$

(36)

which arises from the KMS condition in the Fourier space, $e^{-\beta\hbar\omega}G^{>}(\omega)\!=\!G^{>}(-\omega)$. Now, considering the fact that ${\tanh{(x)}}/{x}\!\leq\!1$, we deduce from Eqs. (35) and (34) that

$\displaystyle{\cal L}_{ww}\leq\frac{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{w}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{2}\bigg{|}_{\mathcal{F}_{w},\mathcal{F}_{q}=0},$

(37)

highlighting a breakdown of the traditional FDR for work current. Importantly, this breakdown is purely of quantum origin, and the equality is restored in the quantum-adiabatic limit, i.e., when $[H_{0},H_{1}]\!=\!0$. This observation aligns with the recent findings of Ref. Miller et al., 2019; Scandi et al., 2020 where a similar breakdown was reported.

The obtained results for the Onsager coefficients and the FDR breakdown for work current have remarkable implications in determining the universal bounds on the performance of the quantum Otto cycle operating as an engine or refrigerator. Particularly, it is noteworthy that the thermodynamic uncertainty relation Sacchi (2021a, b) holds for the currents $\langle\mathcal{J}_{w}\rangle$ and $\langle\mathcal{J}_{h}\rangle$, as shown in Appendix D,

$\displaystyle\langle\sigma\rangle\frac{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{w}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{\langle\mathcal{J}_{w}\rangle^{2}}\geq 2,\quad\langle\sigma\rangle\frac{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{h}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{\langle\mathcal{J}_{h}\rangle^{2}}\geq 2,$

(38)

where $\langle\sigma\rangle=\langle\Sigma\rangle/\tau_{\textrm{cyc}}$ represents average entropy production rate. Notably, these TURs remain valid independently of the nature of the Otto cycle’s operation, i.e., whether it acts as an engine, refrigerator, heater, or accelerator Solfanelli et al. (2020).

Let us now focus on understanding bounds on nonequilibrium fluctuations when the Otto cycle works as an engine, i.e., converting a fraction of spontaneous heat flux flowing from hot to cold reservoir into useful work. To quantify the engine’s performance, we introduce the power output as

$$\langle\mathcal{P}\rangle=-\langle w\rangle/\tau_{cyc}=\!-\frac{{\cal F}_{w}}{\beta}\langle\mathcal{J}_{w}\rangle$$

(39)

and the input heat current is given by $\langle\mathcal{J}_{h}\rangle$. The engine operational regime of the Otto cycle is characterized by $\langle\mathcal{P}\rangle\!>\!0$, and $\langle\mathcal{J}_{h}\rangle\!>\!0$, as per our sign convention where energy flowing into the working fluid is considered positive. To achieve this, the necessary requirement is $\mathcal{F}_{q}\!>\!\mathcal{F}_{w}$. Following the second law of thermodynamics, $\langle\sigma\rangle\!\geq\!0$, it is straightforward to show that the average efficiency of the Otto engine, $\langle\eta\rangle\!\coloneqq\!\langle\mathcal{P}\rangle/\langle\mathcal{J}_{h}\rangle$, is universally upper bounded by the the Carnot bound $\eta_{c}\!=\!\Delta\beta/\beta$ Callen et al. (1985); Lieb and Yngvason (1999), $\langle\eta\rangle\!\leq\!\eta_{c}$. In order to derive bounds on nonequilibrium fluctuations, following Ref. Saryal et al. (2021), we construct the quantity $\eta^{(2)}$, which is the ratio of power fluctuation to the fluctuation of input heat current,

$\displaystyle\eta^{(2)}=\frac{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{P}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{h}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}.$

(40)

Recently, it was shown that in the case of a continuously coupled autonomous steady-state heat engine operating in the linear response regime, where both the work and heat FDRs are valid, the ratio of output power to input heat current fluctuations is subject to both universal lower and upper bounds Saryal et al. (2021),

$\displaystyle\langle\eta\rangle^{2}\leq\big{[}\eta^{(2)}\big{]}_{\mathrm{auto}}\leq\eta_{c}^{2},$

(41)

where the subscript ‘$\mathrm{auto}$’ stands for the autonomous steady-state heat engine. However, in our case of a generic finite-time Otto engine operating in the linear response limit, we obtain the following inequality (see Appendix E for details),

$$\eta_{c}^{2}\geq\frac{{\cal F}_{w}^{2}{\cal L}_{ww}}{\beta^{2}{\cal L}_{qq}}\geq\langle\eta\rangle^{2}.$$

(42)

A key point to note here is that, contrary to the autonomous steady-state engine, we observe a violation of the traditional work FDR in the discrete Otto cycle, as given in Eq. (37). This FDR breakdown for the work current leads to the following inequality:

$\displaystyle\eta^{(2)}\geq\frac{{\cal F}_{w}^{2}{\cal L}_{ww}}{\beta^{2}{\cal L}_{qq}}.$

(43)

Importantly, the equality is restored when the driving Hamiltonian $H_{1}$ commutes with the free Hamiltonian $H_{0}$, i.e., in the quantum-adiabatic driving limit. As an immediate consequence of the inequality in Eq. (43), the lower bound on $\eta^{(2)}$ remains robust,

$\displaystyle\eta^{(2)}\geq\langle\eta\rangle^{2}.$

(44)

However, the upper bound on $\eta^{(2)}$ may not always hold true, i.e., $\eta^{(2)}\nleq\eta_{c}^{2}$, the reason being purely of quantum origin. This distinction is a significant departure from the autonomous case, where both the lower and upper bounds remain intact. In the Otto cycle, due to the violation of the traditional work FDR, the upper bound is not guaranteed, highlighting the impact of nonequilibrium and quantum effects on the engine’s performance. This constitutes another central result of this work. As a direct consequence of the lower bound for $\eta^{(2)}$ in Eq. (44), a hierarchical relation between the TURs is observed. This follows from the fact that

$$\eta^{(2)}\!\geq\!\langle\eta\rangle^{2}\Rightarrow\frac{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{w}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{\langle\mathcal{J}_{w}\rangle^{2}}\!\geq\!\frac{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{h}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{\langle\mathcal{J}_{h}\rangle^{2}}.$$

(45)

As a result, in the engine regime, a strict hierarchy between the TURs of work (output) and heat (input) currents is obtained:

$\displaystyle\langle\sigma\rangle\frac{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{w}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{\langle\mathcal{J}_{w}\rangle^{2}}\geq\langle\sigma\rangle\frac{\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{h}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{\langle\mathcal{J}_{h}\rangle^{2}}\!\geq\!2.$

(46)

After our exploration of the engine regime, we now shift our focus to the refrigerator regime. In this mode of operation, the primary objective is to extract heat from the cold reservoir by utilizing external work. Following our sign convention, i.e., energy entering the working fluid is considered positive, the refrigerator regime is characterized by a positive current flowing out of the cold reservoir, $\langle\mathcal{J}_{c}\rangle\!>\!0$, and work is performed on the system, $\langle\mathcal{J}_{w}\rangle\!>\!0$. To achieve refrigeration, it is necessary that the affinities satisfy $\mathcal{F}_{w}\!>\!\mathcal{F}_{q}$. The performance of the refrigerator is quantified by its coefficient of performance (COP), denoted as $\langle\varepsilon\rangle$. Notably, the COP of a refrigerator is universally upper bounded by the Carnot COP Callen et al. (1985); Abah and Lutz (2016), i.e.,

$\displaystyle\langle\varepsilon\rangle\coloneqq\frac{\beta\langle\mathcal{J}_{c}\rangle}{\mathcal{F}_{w}\langle\mathcal{J}_{w}\rangle}\leq\varepsilon_{c},$

(47)

where $\varepsilon_{c}\!=\!(1-\eta_{c})/\eta_{c}\approx\beta/\Delta\beta$. Let us now investigate the quantity $\varepsilon^{(2)}$ defined as

$\displaystyle\varepsilon^{(2)}\coloneqq\frac{\beta^{2}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{c}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}{\mathcal{F}_{w}^{2}\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{w}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}}.$

(48)

In Appendix F we have derived the expression of CGF for the heat extracted from the cold reservoir, from which we can obtain the expressions for the average heat extracted and its variance. It is worth mentioning that for steady-state autonomous refrigerators operating in the linear response regime, the quantity $\varepsilon^{(2)}$ was reported to possess both universal upper and lower bounds Mohanta et al. (2022),

$$\varepsilon_{c}^{2}\geq\big{[}\varepsilon^{(2)}\big{]}_{\mathrm{auto}}\geq\langle\varepsilon\rangle^{2}.$$

(49)

where, recall that the subscript ‘$\mathrm{auto}$’ denotes the autonomous steady-state refrigerator. In our case of a generic Otto cycle operating as a refrigerator, we derived the following inequality (see Appendix G for details):

$$\varepsilon_{c}^{2}\geq\frac{\beta^{2}{\cal L}_{qq}}{{\cal F}_{w}^{2}{\cal L}_{ww}}\geq\langle\varepsilon\rangle^{2}.$$

(50)

Interestingly, the breakdown of the traditional work FDR leads to the following inequality for the quantity $\varepsilon^{(2)}$:

$\displaystyle\varepsilon^{(2)}\leq\frac{\beta^{2}{\cal L}_{qq}}{{\cal F}_{w}^{2}{\cal L}_{ww}},$

(51)

where it is worth noting that the opposite inequality arises due to the presence of $\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}\mathcal{J}_{w}^{2}\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$ in the denominator of the definition of $\varepsilon^{(2)}$. As a consequence, we observe a contrasting trend compared to the engine regime. Specifically, in the refrigeration regime, the lower bound may not always be valid, i.e., $\varepsilon^{(2)}\!\ngeq\!\langle\varepsilon\rangle^{2}$, whereas the upper bound remains intact,

$$\varepsilon^{(2)}\!\leq\!\varepsilon_{c}^{2}.$$

(52)

Furthermore, it is essential to mention that the potential violation of the lower bound for $\varepsilon^{(2)}$ leads to the loss of the hierarchy in the TURs within the refrigerator regime, unlike what we observed in the engine regime of the Otto cycle as given in Eq. (46). Intriguingly, a similar observation was recently reported for a specific model-dependent periodically driven continuous machine Das et al. (2023). To summarize, the violation of the traditional FDR for the work current in the discrete Otto cycle gives rise to distinct and remarkable consequences for the engine and refrigerator regimes, thereby revealing the impact of quantum effects on the discrete Otto cycle. It is only in the quantum-adiabatic driving limit, i.e., when $[H_{0},H_{1}]\!=\!0$, that both the upper and lower bounds on $\eta^{(2)}$ and $\varepsilon^{(2)}$ are restored, akin to what is observed in continuously coupled autonomous thermal machines.

To illustrate our findings, we now examine a specific model example. Our working medium consists of $N$ non-interacting bosons confined in a one-dimensional parabolic trap with frequency $\omega_{0}$. During the unitary strokes, the potential is perturbed with time. The unperturbed and the perturbing Hamiltonians of the working medium are given by

	$\displaystyle H_{0}$	$\displaystyle=\sum_{i=1}^{N}\frac{p_{i}^{2}}{2m}+\frac{1}{2}m\omega_{0}^{2}x_{i}^{2},$		(53)
	$\displaystyle H_{1}$	$\displaystyle=\sum_{i=1}^{N}m\omega_{0}^{2}x_{i}^{2},$		(54)

respectively, such that $H(t)\!=\!H_{0}+\lambda(t)H_{1},$ where, $\lambda(0)\!=\!\lambda$ and $\lambda(\tau)\!=\!\lambda+\Delta\lambda$. In the second-quantized formulation, we have

	$\displaystyle H_{0}$	$\displaystyle=\sum\nolimits_{m}\epsilon_{m}a_{m}^{\dagger}a_{m},$		(55)
	$\displaystyle H_{1}$	$\displaystyle=\frac{\hbar\omega_{0}}{2}\sum\nolimits_{m}\sqrt{(m+1)(m+2)}\big{(}a_{m}^{\dagger}a_{m+2}+a_{m+2}^{\dagger}a_{m}\big{)}$
		$\displaystyle\quad\quad\quad\quad+(1+2m)a_{m}^{\dagger}a_{m},$		(56)

where $m\!\in\!\mathbb{N}$. Here $\epsilon_{m}\!=\!\hbar\omega_{0}(m+\frac{1}{2})$ represents the single-particle energy-spectrum of $H_{0}$ and $a_{m}$ ($a_{m}^{\dagger}$) is the bosonic annihilation (creation) operator for the $m^{\mathrm{th}}$ energy level. Due to the number conservation, i.e., $\sum_{m}a^{\dagger}_{m}a_{m}\!=\!N$, we are restricted to work with the canonical ensemble description. The central repercussion of this constraint is that it induces correlations between the occupation numbers of different single-particle occupation states. However, one can still compute the canonical partition function $Z_{N}$ for $N$ bosonic particles via a set of recursion relations starting from $N\!=\!1$ Barghathi et al. (2020); Magnus et al. (2017); Mullin and Fernández (2003); Giraud et al. (2018); Borrmann et al. (1999); Schönhammer (2017); Grabsch et al. (2018),

$\displaystyle Z_{N}=$

$\displaystyle\frac{1}{N}\sum_{k=1}^{N}Z_{1}(k\beta)Z_{N-k}.$

(57)

Here, $Z_{1}(k\beta)\!=\!\sum\nolimits_{m}e^{-k\beta\epsilon_{m}}$ represents the partition function for a single-particle state. For our model example, the Green’s function $G^{>}(\omega)$ depends on the two-point correlation involving the occupation number operator. Let $n_{m}\!=\!a_{m}^{\dagger}a_{m}$ be the occupation number operator corresponding to the $m^{\mathrm{th}}$ energy level, and $\overline{\cdot\cdot}^{(\!N\!)}$ represents thermal averaging with respect to the canonical $N$-particle density matrix

$\displaystyle\rho_{N}=\frac{1}{Z_{N}}e^{-\beta\sum\nolimits_{m}\epsilon_{m}n_{m}}\delta_{\sum\nolimits_{m}n_{m},N}.$

(58)

The average occupation of $m^{\mathrm{th}}$ level can be computed using the following recursion relations Barghathi et al. (2020)

	$\displaystyle\overline{n}^{(\!N\!)}_{m}=$	$\displaystyle\frac{1}{Z_{N}}\sum_{k=1}^{N}e^{-\beta\epsilon_{m}k}Z_{N-k},$		(59)
	$\displaystyle\overline{n}^{(\!N+1\!)}_{m}=$	$\displaystyle\frac{Z_{N}}{Z_{N+1}}e^{-\beta\epsilon_{m}}\big{(}1+\overline{n}^{(\!N\!)}_{m}\big{)}.$		(60)

starting from $Z_{0}\!=\!1$ and $\overline{n}^{(0)}_{m}\!=\!0$.