# Monitoring quantum Otto engines

Jeongrak Son  Department of Physics and Astronomy, Seoul National University,
Seoul 08826, Republic of Korea.  Center for Theoretical Physics of Complex
Systems, Institute for Basic Science (IBS), Daejeon 34126, Republic of Korea.
Peter Talkner<EMAIL_ADDRESS>Institut für Physik,
Universität Augsburg, Universitätsstraße 1, D-86135 Augsburg, Germany.  Juzar
Thingna<EMAIL_ADDRESS>Center for Theoretical Physics of Complex
Systems, Institute for Basic Science (IBS), Daejeon 34126, Republic of Korea.
Basic Science Program, Korea University of Science and Technology, Daejeon
34113, Republic of Korea.

###### Abstract

Unlike classical systems, a measurement performed on a quantum system always
alters its state. In this work, the impacts of two diagnostic schemes to
determine the performance of quantum Otto heat engines are compared: In one
scheme the energy of the engine’s working substance is measured after each
stroke (repeated measurements) and in the other one the energies after each
stroke are recorded in one or two pointer states and measured only after the
completion of a prescribed number of cycles (repeated contacts). A single
pointer state suffices if one is only interested in either work or heat. For a
joint work and heat diagnostics two pointers are needed. These schemes are
applied to Otto engines, whose working substance consists of a two-level
system. Depending on the engine protocol, the duration of a single cycle may
be infinite or finite. Because in the repeated contact scheme the number of
measurements is drastically reduced compared to the repeated measurement
scheme, the quantum coherence after, but also during the contact diagnostics
is much better maintained as compared to repeated measurements which destroy
any coherence at the end of each stroke. We demonstrate that maximum power,
reliability, and efficiency of the engine in presence of repeated contacts
typically outperform these figures of merit of repeated measurements. Due to
the improved coherence persistence, heat engines with a finite cycle duration
require a larger number of cycles to reach a periodically asymptotic state.
Overall, our results document the importance of taking into account the
particular nature of diagnostic tools for monitoring and testing purposes but
also for feedback control, both in theory and experiment.

##  I Introduction

Heat engines have not only been important devices from a technological
perspective, but they also played a significant role in establishing and
examining the laws of thermodynamics [1]. Remarkably, with the extension of
thermodynamics from macroscopic to ever smaller systems during the past
decades, the interest in the same paradigmatic model systems such as Carnot
and Otto engines has revived, with efficiency and power as figures of merit. A
particular aspect of small engines results from the generic randomness of
their dynamics due to the coupling of the microscopically small working
substance with heat baths at different temperatures [2]. Already shortly after
the invention of the maser the similarity of its action principle with a heat
engine was noted [3]. Apart from scattered early considerations [4, 5, 6], it
took about forty years until this concept gained broader interest and the idea
of quantum Carnot and Otto engines was further elaborated [7, 8]. Since then
basically all of the distinguishing features of quantum systems have been
investigated in view of their potentially beneficial or possibly detrimental
aspects for an engine, such as coherence [9, 10, 11, 12, 13], discreteness and
spacing of energy levels [14, 15, 16], quantum correlations and entanglement
[17, 18], Fermi versus Bose statistics [19], and projective measurements as an
energy source [20, 21]. Experimental investigations of quantum heat engines
were suggested in Ref. [22] and realized using trapped ions [23, 24],
optomechanical systems [25], nitrogen-vacancy centers [26], and NMR techniques
[27].

So far, the monitoring of engines has found only relatively little attention,
be it as a merely diagnostic tool or as a basic element of a feedback control
device [28]. The functioning of four stroke Otto-like engines during a single
cycle was tracked in Refs. [29, 21, 30] by projective energy measurements of
the working substance immediately after the completion of each stroke. From
these four measured values the energy changes during work strokes and
thermalization strokes corresponding to work and heat, respectively, were
characterized by their joint statistics. The reduction of the engine’s working
substance state with respect to the instantaneous energy basis specifying each
measurement, in general alters the subsequent dynamics and hence may influence
the resulting energy differences between subsequent measurements. This total
suppression of coherence is specific for the two point projective measurement
scheme of work [31, 32, 33] and has often been interpreted as a deficiency
[34, 35]. On the other hand, the traditional assignment of average work done
on a thermally closed system by a control parameter change as the difference
of the average energies calculated after and before the forcing [4, 36] cannot
be extended to the definition of a proper random variable [37, 38, 39, 40].
Only if the interference of the system with a particular monitoring device is
specified, e.g. in terms of two projective energy measurements in the two
point measurement scheme or by two subsequent interactions of the system with
the same pointer state of a measurement apparatus [38], or by any other
generalized measurement scheme [41], work can be properly specified as a
random variable at the price that its probability distribution depends on the
chosen measurement procedure. This principle restriction still leaves the
possibility to look for detection schemes that optimize particular features
such as the accuracy or the amount of back-action on the system, to name just
two possibly complementary criteria.

For the diagnostics of a reciprocating engine such as an Otto engine, the
monitoring of the energy of the working substance several times during a
single cycle is crucial in order to be able to infer the total performed work
and the consumed heat during a single, and also over many cycles. With
projective energy measurements any coherence possibly built up during a work
stroke [11] or after a thermalization stroke [42] would be completely removed
with an impact on the performance of the subsequent strokes. The coherencies
of the bath have been identified as a potential boost to the performance of an
engine [9], while those built up in the working substance have been found as
potentially detrimental [8, 43, 44, 11].

In the present paper we adapt a recently proposed method to determine the sum
of an observable by means of repeated contacts [45] to find the total work and
the heat supplied to an engine in a prescribed number of cycles. We compare
the obtained results with those from repeated measurements and demonstrate
that the repeated contact method is in general less invasive and also leads to
a better performance of the engine.

We organize the paper as follows: In Sec. II we introduce the general quantum
Otto engine and formulate the repeated contact and repeated measurement
schemes. General expressions for the joint work and heat probability density
functions (PDFs) are presented for both schemes. In Sec. III, a two-level
system, i.e. a single qubit, is considered as the working substance for which
general expressions for the work strokes are presented. For the heat strokes
we consider perfect thermalization beyond weak system-bath coupling and
alternatively finite time thermalization at weak system-bath coupling. We
compare our analytic expressions of work and heat PDFs as well as the second
moments resulting from the two measurement schemes for a single engine cycle.
In Sec. IV numerical results for several cycles are presented. Finally, in
Sec. V, we draw our conclusions.

##  II Quantum Otto engines

###  II.1 Mode of operation

Figure 1: Illustration for repeated measurements and two-pointer repeated
contacts applied to quantum Otto engines running for $N$ cycles. The top panel
illustrates repeated measurements with $4N$ pointers, five pointers of which
are shown. Each pointer is initialized in the same state $\sigma$. The
measurements are taken at the beginnings and the end of the work strokes when
the working substance is thermally isolated. The results of the measurements
are symbolically presented at the right margin. The registered pointer states
can be used to calculate heat and work. The lower panel illustrates the
repeated contact scheme with two pointers, the upper one for the work and the
lower one for the heat. Again both pointers are initialized in the same state
$\sigma$. Here the actions of the contacts on the two pointers are accumulated
with appropriate signs, see the Sec. II.2 for the details, and read out by a
projective position measurement after a prescribed number of engine cycles is
completed. Between subsequent contacts the states of the two pointers do not
evolve.

For the sake of concreteness we consider an Otto engine with a quantum
mechanical working substance, which later will be specified as a two-level
system. The subsequent four strokes constituting a full cycle of the engine
consist in compressing, heating, expanding, and cooling the working substance,
see Fig. 1. The compression and expansion strokes are caused by changing a
control parameter of the working substance in thermal isolation. Therefore,
the energy change of the working substance during these strokes can be
referred to as work. The parameter change in the compression stroke is
designed such that the energy level distances of the working substance
Hamiltonian increase. Hence, the time evolution of the density matrix caused
by the compression is determined by a unitary operator $U$ with
$\rho_{2}=U\rho_{1}U^{\dagger}$, with $\rho_{k}=\rho(t_{k})$, $k=1,2,\cdots$,
denoting the density matrices at the times indicated in Fig. 1, whereby times
$t_{k}$ with $k=1,\cdots,4$ specify instants within the first cycle, with
$k=5,\cdots,8$ within the second cycle, etc.. For the sake of simplicity, the
expansion is assumed to be the time-reversed protocol of the compression.
Therefore, the unitary operator giving the corresponding transformation can be
written as $\tilde{U}=K^{*}U^{\dagger}K$ with $K^{*}$ being the adjoint of the
antiunitary operator $K$ ($K^{*}K=KK^{*}=\mathbb{1}$) such that
$\rho_{4}=\tilde{U}\rho_{3}\tilde{U}^{\dagger}$. The relation between
$\tilde{U}$ and $U^{\dagger}$ holds as long as $K^{*}H(t)K=H(t)$ for all time
$t$ during the work strokes with $H(t)$ being the time-dependent engine
Hamiltonian. During the second and the fourth stroke the working substance is
brought into contact with a hot and a cold heat bath, respectively, while its
Hamiltonian remains unchanged. These two heat strokes are characterized by
trace preserving, positive, linear maps $\Phi_{h}$ and $\Phi_{c}$ relating the
respective density matrices of the working substance at the beginnings and the
ends of these strokes according to $\rho_{3}=\Phi_{h}(\rho_{2})$ and
$\rho_{5}=\Phi_{c}(\rho_{4})$. Later, we shall give explicit expressions
specifying the unitary operator $U$, $\tilde{U}$ and the maps $\Phi_{h}$ and
$\Phi_{c}$ for the case of a two-level system as working substance.

###  II.2 Diagnostic means I: Repeated measurements

The central quantities characterizing the performance of an engine are the
total work $W$ performed on and the total heat $Q$ supplied by the hot bath to
the working substance during $N$ cycles. Both quantities follow from the
sequence of energies $\mathbf{e}=(e_{1},e_{2},\cdots,e_{4N})$ of the working
substance taken at the beginning and the end of each work stroke when the
working substance is decoupled from either heat bath yielding

$\displaystyle W$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{4N}(-1)^{k}e_{k},$ (1) $\displaystyle Q$
$\displaystyle=$ $\displaystyle\sum_{k=1}^{4N}\mu_{k}e_{k},$ (2)

where $\mu_{1+4l}=\mu_{4+4l}=0$, $\mu_{2+4l}=-1$, $\mu_{3+4l}=1$ with
$l=0,\cdots,N-1$.

In the repeated measurement scheme the energies $e_{k}$ are measured
individually by bringing the working substance for a short time into contact
with the pointer of a measurement apparatus [46]. This contact causes a
unitary transformation of the joint state of the pointer and the working
substance given by

$V_{k}=e^{-i\kappa H_{k}\otimes P_{k}},$ (3)

where $H_{k}$ is the working substance Hamiltonian at the start of the $k$th
stroke, i.e. $H_{1}=H_{4}=H_{5}=\cdots=H_{c}$ and $H_{2}=H_{3}=\cdots=H_{h}$,
see also Fig. 1; furthermore, $P_{k}$ denotes the momentum operator conjugate
to the pointer coordinate of the $k$th measurement apparatus. The parameter
$\kappa$ results as the product of the strength and the duration of the
contact which is assumed to be so short that the free motion of the working
substance appears as frozen. Above, and throughout this work we will set
$\hbar=k_{B}=1$. Before the interaction, the working substance and the pointer
are uncorrelated. The latter resides in a Gaussian pure state $\sigma_{k}$
with vanishing mean value and variance $\Sigma^{2}$. Its position matrix
element $\sigma(x,y)=\langle x|\sigma_{k}|y\rangle$ is independent of $k$ and
is given by

$\sigma(x,y)=\frac{1}{\sqrt{2\pi\Sigma^{2}}}\exp\left[-\frac{x^{2}+y^{2}}{4\Sigma^{2}}\right].$
(4)

If, after the contact has taken place, the pointer position is measured
projectively at a position $x$, the non-normalized density matrix of the
working substance becomes

$\displaystyle\phi_{x}(\rho)$ $\displaystyle=$ $\displaystyle\langle
x|V_{k}\rho\otimes\sigma_{k}V^{\dagger}_{k}|x\rangle$ (5) $\displaystyle=$
$\displaystyle\sum_{m,m^{\prime}}\mathcal{P}^{m}_{k}\rho\mathcal{P}^{m^{\prime}}_{k}\sigma_{k}\left(x-e^{m}_{k},x-e^{m^{\prime}}_{k}\right),$

where $e^{m}_{k}$ and $\mathcal{P}^{m}_{k}$ denote the eigenvalues and eigen-
projectors of the Hamiltonian $H_{k}=\sum_{m}e^{m}_{k}\mathcal{P}^{m}_{k}$ and
$\rho$ denotes the density matrix of the working substance before the contact.
For the sake of simplicity and without restriction of generality, the pointer
position is measured in units of energy such that $\kappa=1$. The subsequent
measurements of the energies $e_{k}$ during $N$ cycles, giving rise to the
pointer positions $\vec{x}_{N}=(x_{1},\cdots,x_{4N})$, lead to a mapping of
the initial working substance density matrix $\rho$ onto the non-normalized
density matrix $\phi^{\mathrm{RM}}_{\vec{x}_{N}}(\rho)$ given by

$\phi^{\mathrm{RM}}_{\vec{x}_{N}}(\rho)=\sum_{\vec{m}_{N},\vec{m}^{\prime}_{N}}D^{\vec{m}_{N},\vec{m}^{\prime}_{N}}(\rho)\prod_{k=1}^{4N}\sigma\left(x_{k}-e_{k}^{m_{k}},x_{k}-e_{k}^{m_{k}^{\prime}}\right),$
(6)

where the sum runs over all possible sequences
$\vec{m}^{(\prime)}_{N}=(m^{(\prime)}_{1},m^{(\prime)}_{2},\cdots,m^{(\prime)}_{4N})$
of quantum numbers $m^{(\prime)}_{k}$ labelling the eigenstates of the
Hamiltonian $H_{k}$. The operator $D^{\vec{m}_{N},\vec{m}^{\prime}_{N}}(\rho)$
describes the action of the alternating time-evolution operators
characterizing the subsequent strokes and the interrupting projections due to
the energy measurements on the initial density matrix $\rho$. They are
recursively given by

$D^{\vec{m}_{N+1},\vec{m}^{\prime}_{N+1}}(\rho)=S^{\mathbf{m},\mathbf{m}^{\prime}}\left(D^{\vec{m}_{N},\vec{m}^{\prime}_{N}}(\rho)\right),$
(7)

where the sequences
$\vec{m}^{(\prime)}_{N+1}=\vec{m}^{(\prime)}_{N}\oplus\mathbf{m}^{(\prime)}$
denote the concatenation of $\vec{m}^{(\prime)}_{N}$ and
$\mathbf{m}^{(\prime)}=(m_{4N+1}^{(\prime)},\cdots,m_{4(N+1)}^{(\prime)})$.
Furthermore, the map $S^{\mathbf{m},\mathbf{m}^{\prime}}$ acts as

$\displaystyle S^{\mathbf{m},\mathbf{m}^{\prime}}(\rho)$ $\displaystyle=$
$\displaystyle\Phi_{c}\left(\mathcal{P}^{m_{4(N+1)}}_{1}\tilde{U}\mathcal{P}^{m_{4N+3}}_{2}\Phi_{h}\left(\mathcal{P}^{m_{4N+2}}_{2}U\mathcal{P}^{m_{4N+1}}_{1}\rho\mathcal{P}^{m_{4N+1}^{\prime}}_{1}U^{\dagger}\mathcal{P}^{m_{4N+2}^{\prime}}_{2}\right)\mathcal{P}^{m_{4N+3}^{\prime}}_{2}\tilde{U}^{\dagger}\mathcal{P}^{m_{4(N+1)}^{\prime}}_{1}\right).$
(8)

Taking the trace of $\phi^{\mathrm{RM}}_{\vec{x}}(\rho)$ one obtains for the
joint PDF $p(\vec{x})$to find the pointers at the values given by the
components of the vector $\vec{x}$ the formal expression

$p(\vec{x})=\mathrm{Tr}\left[\phi^{\mathrm{RM}}_{\vec{x}}(\rho)\right].$ (9)

For Gaussian pointers as specified in Eq. (4) this PDF can be written as

$p(\vec{x})=\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}\sum_{k=1}^{4N}\left(e^{m_{k}}_{k}-e^{m_{k}^{\prime}}_{k}\right)^{2}}\prod_{k=1}^{4N}g_{\Sigma^{2}}\left(x_{k}-\frac{1}{2}\left(e^{m_{k}}_{k}+e^{m_{k}^{\prime}}_{k}\right)\right),$
(10)

where $g_{\Sigma^{2}}(x)=(2\pi\Sigma^{2})^{-1/2}\exp[-x^{2}/(2\Sigma^{2})]$
denotes a Gaussian PDF with vanishing mean value and variance $\Sigma^{2}$.
The coefficient $D^{\vec{m},\vec{m}^{\prime}}$ is defined as

$D^{\vec{m},\vec{m}^{\prime}}=\mathrm{Tr}[D^{\vec{m},\vec{m}^{\prime}}(\rho)].$
(11)

In order not to overburden the notation we omitted the index $N$ in the
vectors $\vec{m}^{(\prime)}$ and $\vec{x}$.

Identifying the components $x_{k}$ as the measured energies, one obtains for
the joint probability of work and heat, $P^{\mathrm{RM}}(W,Q)$, the expression

$\displaystyle P^{\mathrm{RM}}(W,Q)$ $\displaystyle=$ $\displaystyle\int
d^{4N}x\,\delta\left(W-\sum_{k=1}^{4N}(-1)^{k}x_{k}\right)$ (12)
$\displaystyle\times\delta\left(Q-\sum_{k=1}^{4N}\mu_{k}x_{k}\right)p(\vec{x}).$

Going to the characteristic function $G^{\mathrm{RM}}(u,v)=\int
dWdQe^{iuW}e^{ivQ}P^{\mathrm{RM}}(W,Q)$ one can perform the $4N$-fold
$x$-integration to find

$\displaystyle G^{\mathrm{RM}}(u,v)$
$\displaystyle=\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{\tilde{m}}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}\sum_{k=1}^{4N}\left(e^{m_{k}}_{k}-e^{m_{k}^{\prime}}_{k}\right)^{2}}$
$\displaystyle\quad\times
e^{iu\mathcal{W}^{\vec{m},\vec{m}^{\prime}}}e^{iv\mathcal{Q}^{\vec{m},\vec{m}^{\prime}}}e^{-N\Sigma^{2}(2u^{2}-2uv+v^{2})},$

where

$\displaystyle\mathcal{W}^{\vec{m},\vec{m}^{\prime}}$
$\displaystyle=\frac{1}{2}\sum_{k=1}^{4N}(-1)^{k}\left(e^{m_{k}}_{k}+e^{m_{k}^{\prime}}_{k}\right),$
(14) $\displaystyle\mathcal{Q}^{\vec{m},\vec{m}^{\prime}}$
$\displaystyle=\frac{1}{2}\sum_{k=1}^{4N}\mu_{k}\left(e^{m_{k}}_{k}+e^{m_{k}^{\prime}}_{k}\right),$
(15)

with $\mu_{k}$ as defined below Eq. (2). This characteristic function
corresponds to a linear superposition of Gaussian PDFs
$g_{\mathcal{M}}(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}},Q-\mathcal{Q}^{\vec{m},\vec{m}^{\prime}})$
with mean values $\langle W\rangle=\mathcal{W}^{\vec{m},\vec{m}^{\prime}}$ and
$\langle Q\rangle=\mathcal{Q}^{\vec{m},\vec{m}^{\prime}}$ and with the
covariance matrix $\mathcal{M}$ given by

$\mathcal{M}=\left(\begin{array}[]{cc}\langle(\delta
W)^{2}\rangle&\langle\delta W\delta Q\rangle\\\ \langle\delta W\delta
Q\rangle&\langle(\delta
Q)^{2}\rangle\end{array}\right)=2N\Sigma^{2}\left(\begin{array}[]{cc}2&-1\\\
-1&1\end{array}\right),$ (16)

with $\delta X=X-\langle X\rangle$ ($X=W,Q$). Hence, the PDF is given by

$\displaystyle P^{\mathrm{RM}}(W,Q)$
$\displaystyle=\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}{\sum_{k=1}^{4N}}(e^{m_{k}}_{k}-e^{m_{k}\prime}_{k})^{2}}$
(17) $\displaystyle\quad\times
g_{\mathcal{M}}(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}},Q-\mathcal{Q}^{\vec{m},\vec{m}^{\prime}}).$

All diagonal coefficients with $\vec{m}=\vec{m}^{\prime}$ are non-negative
whereas nondiagonal contributions with $\vec{m}\neq\vec{m}^{\prime}$ may
assume complex values. Yet, by construction, the total sum is guaranteed to be
always non-negative. Moreover, the nondiagonal contributions are suppressed by
exponential factors
$\exp[-\sum_{k=1}^{4N}\left(e^{m_{k}}_{k}-e^{m_{k}^{\prime}}_{k}\right)^{2}/(8\Sigma^{2})]$.
When the pointer state variance $\Sigma^{2}$ is sufficiently small compared to
the minimum squared energy level distance of either Hamiltonian $H_{c}$ or
$H_{h}$, all non-diagonal elements are negligibly small and hence only
Gaussians having a maximum at
$\mathcal{W}^{\vec{m},\vec{m}},\mathcal{Q}^{\vec{m},\vec{m}}$ contribute. Due
to the $N$-dependence of the covariance matrix the widths of each of the
contributions grow proportionally to the square root of the number of cycles.
The marginal work and heat PDFs result as

$\displaystyle P^{\mathrm{RM}}(W)$
$\displaystyle=\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}{\sum_{k=1}^{4N}}(e^{m_{k}}_{k}-e^{m_{k}\prime}_{k})^{2}}g_{4N\Sigma^{2}}(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}}),$
(18) $\displaystyle P^{\mathrm{RM}}(Q)$
$\displaystyle=\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}{\sum_{k=1}^{4N}}(e^{m_{k}}_{k}-e^{m_{k}\prime}_{k})^{2}}g_{2N\Sigma^{2}}(Q-\mathcal{Q}^{\vec{m},\vec{m}^{\prime}}).$
(19)

The more pronounced broadening of the Gaussian contributions in the marginal
work PDF compared to the marginal heat PDF is caused by the twofold number of
energy measurements required for the work compared to the heat.

###  II.3 Diagnostic means II: Repeated contacts

Instead of registering each of the $4N$ energy values in a separate pointer
state one may follow the idea of Ref. [47] and subsequently transfer the
actual energy value at the beginning of each stroke with the appropriate sign
according to the Eqs. (1) and (2) to two pointer states, one for work and the
other one for heat. The transfer of information from the working substance to
the pointer states is mediated in a similar way as for repeated measurements
by a sufficiently short interaction causing a unitary transformation $V_{k}$
of the form

$V_{k}=e^{i\kappa H_{k}\otimes((-1)^{k}P_{W}+\mu_{k}P_{Q})},$ (20)

acting on the joint Hilbert space of the working substance and the two pointer
states. Here $P_{W}$ and $P_{Q}$ denote the momentum operators that are
canonically conjugate to the position operators of the work and heat pointers,
respectively. In contrast to the repeated measurement scheme, the pointers are
only read out by a projective position measurements after the completion of
$N$ cycles. Again, the pointers are assumed to be gauged in units of energy.
Between two subsequent contacts the working substance does not interact with
the two pointers but their states remain correlated as a result of the
previous contacts. To simplify the further analysis we assume that the
pointers have no own dynamics, i.e. they change their states only due to the
interactions as specified by the unitary operators defined in Eq. (20). The
subsequent strokes and contacts, resulting after the completion of $N$ cycles
in pointer state positions at specific values $W$ and $Q$, lead to a non-
normalised reduced density matrix of the working substance, which is given by

$\displaystyle\phi^{\mathrm{RC}}_{W,Q}(\rho)$ $\displaystyle=$
$\displaystyle\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}(\rho)\sigma(W-\mathcal{W}^{\vec{m}},W-\mathcal{W}^{\vec{m}^{\prime}})$
(21)
$\displaystyle\times\sigma(Q-\mathcal{Q}^{\vec{m}},Q-\mathcal{Q}^{\vec{m}^{\prime}}),$

where $\rho$ denotes the initial density matrix of the working substance and
individual work and heat outcomes are denoted as

$\displaystyle\mathcal{W}^{\vec{m}}$
$\displaystyle=\sum_{k=1}^{4N}(-1)^{k}e^{m_{k}}_{k},$ (22)
$\displaystyle\mathcal{Q}^{\vec{m}}$
$\displaystyle=\sum_{k=1}^{4N}\mu_{k}e^{m_{k}}_{k},$ (23)

with $\mu_{k}$ as defined below Eq. (2). The joint PDF $P^{\mathrm{RC}}(W,Q)$
of finding these work and heat values is given by the trace of
$\phi^{\mathrm{RC}}_{W,Q}(\rho)$ yielding

$P^{\mathrm{RC}}(W,Q)=\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{W}^{\vec{m}}-\mathcal{W}^{\vec{m}\prime})^{2}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{Q}^{\vec{m}}-\mathcal{Q}^{\vec{m}\prime})^{2}}g_{\Sigma^{2}}(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}})g_{\Sigma^{2}}(Q-\mathcal{Q}^{\vec{m},\vec{m}^{\prime}}).$
(24)

As for repeated measurements it is a linear superposition of Gaussian PDFs
with weights $D^{\vec{m},\vec{m}^{\prime}}$, however, with the following
differences:

  1. 1.

for repeated contacts the contributing Gaussians for work and heat factorize
while for repeated measurements they are correlated;

  2. 2.

the widths of the contributing Gaussians are independent of $N$ given by
$\Sigma$ whereas those for repeated measurements are proportional to
$\sqrt{N}\Sigma$;

  3. 3.

for repeated contacts nondiagonal contributions are exponentially suppressed
by rates
$[(\mathcal{W}^{\vec{m}}-\mathcal{W}^{\vec{m}^{\prime}})^{2}+(\mathcal{Q}^{\vec{m}}-\mathcal{Q}^{\vec{m}^{\prime}})^{2}]/(8\Sigma^{2})$
while for repeated measurements the suppression is determined by rates of the
form
$\sum_{k=1}^{4N}(e^{m_{k}}_{k}-e^{m_{k}^{\prime}}_{k})^{2}/(8\Sigma^{2})$.

Therefore, for repeated contacts only contributions with different sequences
$\vec{m}$ and $\vec{m}^{\prime}$ leading to different values of work and heat
are penalized while in the repeated measurement approach deviations of
individual components lead to a suppression of coherences due to the
monitoring.

For repeated contacts with a non-selective heat measurement the marginal PDF
$P^{\mathrm{RC,2P}}(W)$ defined as

$P^{\mathrm{RC,2P}}(W)=\int
dQP^{\mathrm{RC}}(W,Q)=\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{W}^{\vec{m}}-\mathcal{W}^{\vec{m}\prime})^{2}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{Q}^{\vec{m}}-\mathcal{Q}^{\vec{m}\prime})^{2}}g_{\Sigma^{2}}(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}})$
(25)

differs by the extra factor of
$\exp[-(\mathcal{Q}^{\vec{m}}-\mathcal{Q}^{\vec{m}\prime})^{2}/(8\Sigma^{2})]$
from the PDF resulting from a mere work measurement performed with a single
pointer state. In the latter case one obtains

$\displaystyle P^{\mathrm{RC,1P}}(W)$ $\displaystyle=$
$\displaystyle\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{W}^{\vec{m}}-\mathcal{W}^{\vec{m}\prime})^{2}}$
(26) $\displaystyle\times
g_{\Sigma^{2}}(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}}).$

For repeated contacts the same type of contextuality also holds for the
marginal heat PDF. In both cases it is more pronounced for rather imprecise
registration, i.e. for pointer states with relatively large variances
$\Sigma^{2}$. For repeated measurements the joint work and heat PDF [see Eq.
(17)] results according to Eq. (12) in a classical way from the joint PDF
$p(\vec{x})$ of all individual energies. Therefore repeated measurements do
not present this particular contextuality feature.

Once the work and heat PDFs are known for a specified diagnostic method, the
performance of the engine can be characterized by its efficiency $\eta$ and
reliability $R$ which are defined as

$\displaystyle\eta$ $\displaystyle=$ $\displaystyle-\frac{\langle
W\rangle}{\langle Q\rangle},$ (27) $\displaystyle R$ $\displaystyle=$
$\displaystyle-\frac{\langle W\rangle}{\sqrt{\langle W^{2}\rangle-\langle
W\rangle^{2}}},$ (28)

where $\langle\cdot\rangle$ denotes the average of the indicated quantity with
respect to the PDF of the considered diagnostic method. Higher order
efficiencies as introduced in Ref. [48] can be obtained from higher order
moments of work and heat. Also a fluctuating efficiency, $W/Q$, can in
principle be determined [49], but will not be considered here.

##  III Two-level system as working substance

###  III.1 Work strokes

The general formulation outlined in Sec. II can be applied to any working
substance undergoing adiabatic or non-adiabatic protocols. For a numerically
feasible but yet representative example, we will use in this and subsequent
sections a two-level system as the working substance and allow for non-
adiabatic work strokes, i.e. strokes comprising shifts of the instantaneous
energy eigenstates and transitions between them during the unitary evolution.
Without specifying the time-dependent Hamiltonian in detail, the system
Hamiltonians $H_{c}$ and $H_{h}$ before and after the first work stroke,
respectively, are written as

$\displaystyle H_{c}$ $\displaystyle=$
$\displaystyle\epsilon_{c}\left(|+_{c}\rangle\langle+_{c}|-|-_{c}\rangle\langle-_{c}|\right),$
$\displaystyle H_{h}$ $\displaystyle=$
$\displaystyle\epsilon_{h}\left(|+_{h}\rangle\langle+_{h}|-|-_{h}\rangle\langle-_{h}|\right),$
(29)

where $|-_{u}\rangle$ and $|+_{u}\rangle$ denote the ground and the excited
states, respectively, of the Hamiltonian $H_{u}$ with $u=h,c$. The most
general form of the unitary evolution, during the work stroke, then reads

$\displaystyle U=$
$\displaystyle\sqrt{1-\alpha}\left(e^{-i\varphi}|+_{h}\rangle\langle+_{c}|+e^{i\varphi}|-_{h}\rangle\langle-_{c}|\right)$
(30)
$\displaystyle-\sqrt{\alpha}\Big{(}|+_{h}\rangle\langle-_{c}|-|-_{h}\rangle\langle+_{c}|\Big{)},$
$\displaystyle\tilde{U}=$
$\displaystyle\sqrt{1-\alpha}\left(e^{-i\varphi}|+_{c}\rangle\langle+_{h}|+e^{i\varphi}|-_{c}\rangle\langle-_{h}|\right)$
(31)
$\displaystyle+\sqrt{\alpha}\Big{(}|+_{c}\rangle\langle-_{h}|-|-_{c}\rangle\langle+_{h}|\Big{)}.$

Above $\tilde{U}\equiv C^{*}U^{\dagger}C$ is the evolution operator for the
time-reversed protocol with the antiunitary operator $K$ chosen to be the
complex conjugation operator $C$. Here $\alpha\in[0,1]$ is the transition
probability between the ground and the excited states and $\varphi$ is the
corresponding phase (both parameters depend on the stroke protocol and in
particular on its duration). The limiting case $\alpha=0$ refers to an
adiabatic process while $\alpha=1$ corresponds to a swap. For example, a
Landau-Zener protocol is described by the time-dependent Hamiltonian

$H(t)=-\frac{vt}{2}\sigma_{z}-\epsilon_{c}\sigma_{x}\>,$ (32)

where the velocity of the drive,
$v=4\sqrt{\epsilon_{h}^{2}-\epsilon_{c}^{2}}/T_{1}$, is chosen such that the
eigenvalues of the Hamiltonian change from $\pm\epsilon_{c}$ at $t=0$ to
$\pm\epsilon_{h}$ at $t=T_{1}/2$, the time corresponding to the duration of a
work stroke. The parameters specifying the unitary operator $U$ are given in
Ref. [50] by

$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle
e^{-2\pi\delta},\quad\delta=\frac{\epsilon_{c}T_{1}}{4\sqrt{(\frac{\epsilon_{h}}{\epsilon_{c}})^{2}-1}},$
$\displaystyle\varphi$ $\displaystyle=$
$\displaystyle\frac{\pi}{4}-\delta(\log{\delta}-1)-\mathrm{arg}\Gamma(1-i\delta).$
(33)

where $\Gamma(z)$ denotes the Gamma function. For other protocols
characterized by a single parameter see Ref. [51] and for more general cases
Ref. [52, 53].

In the first instance we shall not refer to a particular work protocol and use
the general form of a unitary time evolution $U$ as specified by Eq. (30).
Before we do so, we consider the heating and cooling strokes of the two-level
working substance that are achieved by a coupling to the heat reservoirs at
temperatures $T_{h}$ and $T_{c}$, respectively. Depending on the duration and
the strength of the respective contacts perfect or imperfect thermalization of
the two-level system can be reached.

###  III.2 Perfect thermalization

Provided that the contact between the working substance and the heat bath is
sufficiently long, the finally reached reduced state of the two-level system
is always the same, dependent on the temperature of the heat bath and the
strength of the interaction but independent of its state $\rho$ at the
beginning of the contact. Hence, the operation $\Phi_{u}$, $u=h,c$ describing
a perfect thermalization stroke is given by the projection onto the final
state $\rho_{\beta}$ and can be written as

$\Phi_{u}(\rho)=\rho_{u}\mathrm{Tr}[\rho].$ (34)

Only for weak coupling $\rho_{u}$, $u=h,c$ coincides with the Gibbs state
$\rho_{u}\propto\exp[-\beta_{u}H_{u}]$ defined by the Hamiltonian of the two-
level system at the inverse temperature $\beta_{u}$. While these states lack
coherence the latter may be found in generalized Gibbs states resulting from
the partial trace of the total system’s Gibbs state with respect to the bath
[40, 42]. For a perturbative treatment in the coupling strength up to second
order see Append. A.

In the asymptotic limit of an initially sharp pointer state with $\Sigma=0$
the joint work and heat PDFs characterizing a single cycle are identical for
repeated measurements and repeated contacts (two-pointer scheme), i.e.,
$P^{\mathrm{RM}}_{1}(W,Q)=P^{\mathrm{RC}}_{1}(W,Q)$. The first two moments can
be expressed as

$\displaystyle\langle W\rangle$ $\displaystyle=$
$\displaystyle\mathcal{A}_{c,h}\,\epsilon_{h}+\mathcal{A}_{h,c}\,\epsilon_{c},\quad\quad\langle
W^{2}\rangle=2\mathcal{B}_{c,h}\left(\epsilon_{c}^{2}+\epsilon_{h}^{2}\right)+\left[\frac{1-\mathcal{B}_{c,h}}{1-2\alpha}-(1-2\alpha)\left(1+\mathcal{B}_{c,h}\right)\right]\epsilon_{c}\epsilon_{h},$
(35) $\displaystyle\langle Q\rangle$ $\displaystyle=$
$\displaystyle-\mathcal{A}_{c,h}\,\epsilon_{h},\quad\quad\quad\quad\quad~{}~{}\langle
Q^{2}\rangle=2\mathcal{B}_{c,h}\,\epsilon_{h}^{2},$ (36) $\displaystyle\langle
WQ\rangle$ $\displaystyle=$
$\displaystyle-2\mathcal{B}_{c,h}\,\epsilon_{h}^{2}-\frac{1}{2}\left[\frac{1-\mathcal{B}_{c,h}}{1-2\alpha}-(1-2\alpha)\left(1+\mathcal{B}_{c,h}\right)\right]\epsilon_{c}\epsilon_{h},$
(37)

with $\mathcal{A}_{c,h}=2\left(\alpha+d_{c}-d_{h}-2\alpha d_{c}\right)$ and
$\mathcal{B}_{c,h}=1-(1-2\alpha)(1-2d_{c})(1-2d_{h})$. Here, $d_{c}$ and
$d_{h}$ denote the populations of the exited states of the cold and the hot
density matrices, respectively, see the Eqs. (53) and (LABEL:d) in Append. A.
Even when the Gaussian pointer has a finite width $\Sigma\neq 0$, the
differences between repeated measurements and contacts in the average values
of work and heat are of the order $\exp[-\epsilon^{2}_{c(h)}/(2\Sigma^{2})]$
which could be extremely small for a precise measurement apparatus, i.e. if
$\epsilon_{c(h)}<\Sigma$. On the other hand, the second moments differ by
$3\Sigma^{2}+\mathcal{O}(\exp[-\epsilon_{c(h)}^{2}/(2\Sigma^{2})])$ which is
dominated by the $3\Sigma^{2}$ term arising due to the difference in the
number of measurements. Accordingly, also the reliability [see Eq. (28)]
differs for repeated contacts and repeated measurements. Overall, however,
even though the two measurement schemes are quite different they do not show
major significant differences for a perfectly thermalizing quantum Otto
engine.

###  III.3 Imperfect thermalization

Perfect thermalization requires long times during which the working substance
is coupled to one of the two heat baths, leading to long cycle periods and
correspondingly small power output of the engine. The dynamics of the
relaxation process in general is very complicated. A substantial
simplification occurs only at week coupling in the Markovian limit in which
the dynamics is governed by a Lindblad equation

$\dot{\rho_{u}}(t)=-i[H_{u},\rho_{u}]+\Gamma_{u}^{+}\mathcal{D}_{u}^{+}[\rho_{u}]+\Gamma_{u}^{-}\mathcal{D}_{u}^{-}[\rho_{u}],$
(38)

with the subscript $u=c,h$, the dissipators
$\mathcal{D}_{u}^{\pm}[\rho_{u}]=2L^{\mp}_{u}\rho_{u}L^{\pm}_{u}-\\{L^{\pm}_{u}L^{\mp}_{u},\rho_{u}\\}$,
the rates $\Gamma_{u}^{\pm}=\gamma\epsilon_{u}/(1+\exp[\mp
2\beta_{u}\epsilon_{u}])$, and the Lindblad jump operators
$L^{\pm}_{u}=|\pm_{u}\rangle\langle\mp_{u}|$. The Hamiltonians $H_{u}$ are
given by Eq. (29). Even though the stationary solutions of the master
equations for the hot and the cold heat baths are given by the Gibbs states,
corresponding to the Hamiltonians $H_{u}$ at the respective temperatures
$\beta_{u}$ and hence do not possess any coherence, the solutions of the
master equations at any finite time in general do contain coherences. This
leads to differences between repeated measurements and repeated contacts even
in the limit of infinite precision ($\Sigma\rightarrow 0$). For a single
cycle, starting at some initial state
$\rho=d|+_{\mathrm{c}}\rangle\langle+_{\mathrm{c}}|+(1-d)|-_{\mathrm{c}}\rangle\langle-_{\mathrm{c}}|+[q|+_{\mathrm{c}}\rangle\langle-_{\mathrm{c}}|+\mathrm{h.c.}]$
one finds the following expressions for the first two moments of work and heat
for a small pointer state variance including leading order
$\mathcal{O}(\Sigma^{2})$ corrections:

$\displaystyle\langle W\rangle^{\mathrm{RM}}$ $\displaystyle=$ $\displaystyle
2\alpha(1-\alpha)\left[\epsilon_{h}\tanh\left(\beta_{h}\epsilon_{h}\right)(1-e^{-2\gamma\theta})+\epsilon_{c}(1-2d)(1+e^{-2\gamma\theta})\right]$
(39)
$\displaystyle+\sum_{j=0,1}\left[(-1)^{j}\epsilon_{h}+\epsilon_{c}\right]\left[1+(-1)^{j}\left(2\alpha-1\right)\right]^{2}\left(\frac{e^{(-1)^{j}\beta_{h}\epsilon_{h}}}{e^{\beta_{h}\epsilon_{h}}+e^{-\beta_{h}\epsilon_{h}}}-d\right)\left(1-e^{-2\gamma\theta}\right),$
$\displaystyle\langle W\rangle^{\mathrm{RC}}$ $\displaystyle=$
$\displaystyle\langle
W\rangle^{\mathrm{RM}}-4\epsilon_{c}\alpha(1-\alpha)(1-2d)e^{-\gamma\theta}\cos\left(2(\theta+\varphi)\right),$
(40) $\displaystyle\langle W^{2}\rangle^{\mathrm{RM}}$ $\displaystyle=$
$\displaystyle
4\left(\epsilon_{h}^{2}+\epsilon_{c}^{2}\right)\bigg{\\{}(\alpha+d-2\alpha
d)+\left[(1-2d)(1-2\alpha)-2\alpha^{2}d\right]\frac{e^{-\beta_{h}\epsilon_{h}}}{e^{\beta_{h}\epsilon_{h}}+e^{-\beta_{h}\epsilon_{h}}}\bigg{\\}}\left(1-e^{-2\gamma\theta}\right)$
(41)
$\displaystyle+8\epsilon_{h}\epsilon_{c}\bigg{\\{}\alpha^{2}(1-2d)-(1-2\alpha)d-\left[(1-2d)(1-2\alpha)-2\alpha^{2}d\right]\frac{e^{-\beta_{h}\epsilon_{h}}}{e^{\beta_{h}\epsilon_{h}}+e^{-\beta_{h}\epsilon_{h}}}\bigg{\\}}\left(1-e^{-2\gamma\theta}\right)$
$\displaystyle+4\epsilon_{c}^{2}\alpha(1-\alpha)e^{-2\gamma\theta}+4\Sigma^{2},$
$\displaystyle\langle W^{2}\rangle^{\mathrm{RC}}$ $\displaystyle=$
$\displaystyle\langle
W^{2}\rangle^{\mathrm{RM}}-8\epsilon_{c}^{2}\alpha(1-\alpha)e^{-\gamma\theta}\cos\left(2(\theta+\varphi)\right)-3\Sigma^{2},$
(42) $\displaystyle\langle WQ\rangle^{\mathrm{RC}}$ $\displaystyle=$
$\displaystyle-4\epsilon_{h}^{2}\bigg{\\{}(\alpha+d-2\alpha
d)+\left[(1-2d)(1-2\alpha)-2\alpha^{2}d\right]\frac{e^{-\beta_{h}\epsilon_{h}}}{e^{\beta_{h}\epsilon_{h}}+e^{-\beta_{h}\epsilon_{h}}}\bigg{\\}}\left(1-e^{-2\gamma\theta}\right)$
(43)
$\displaystyle-4\epsilon_{h}\epsilon_{c}\bigg{\\{}\alpha^{2}(1-2d)-(1-2\alpha)d-\left[(1-2d)(1-2\alpha)-2\alpha^{2}d\right]\frac{e^{-\beta_{h}\epsilon_{h}}}{e^{\beta_{h}\epsilon_{h}}+e^{-\beta_{h}\epsilon_{h}}}\bigg{\\}}\left(1-e^{-2\gamma\theta}\right),$
$\displaystyle\langle WQ\rangle^{\mathrm{RM}}$ $\displaystyle=$
$\displaystyle\langle WQ\rangle^{\mathrm{RC}}-2\Sigma^{2},$ (44)

where $\theta=\epsilon_{h}\tau_{h}$ with $\tau_{h}$ being the duration of the
hot bath stroke. A dependence of the moments on the initial coherence
parameter $q$ is only seen in higher order $\Sigma^{2}$ corrections. From the
expressions of these moments the according variances and the covariance of
work and heat can be obtained. Instead of these bulky expressions we present
the variances and covariances for the two monitoring schemes in Fig. 2a, b as
functions of $\theta$. The difference between repeated measurements and
repeated contacts is governed by damped oscillations stemming from the
$e^{-\gamma\theta}\cos\left(2(\theta+\varphi)\right)$-term in the deviation of
the first two work moments, see Eqs. (40), (42) and Fig. 2c, d. Yet, the
maximal amplitude of these oscillations is rather small rendering the
(co)variances for repeated measurements and repeated contacts almost
identical.

Figure 2: The variance of work $\langle\delta W^{2}\rangle=\langle
W^{2}\rangle-\langle W\rangle^{2}$ and the covariance of work and heat,
$\langle\delta W\delta Q\rangle=\langle WQ\rangle-\langle W\rangle\langle
Q\rangle$, for repeated contacts and repeated measurements for a single cycle
as a function of the dimensionless time of contact with the heat baths,
$\theta=\epsilon_{h}\tau_{h}=\epsilon_{c}\tau_{c}$ are drawn in panels a and b
respectively. The dashed horizontal lines represent the asymptotic values in
the perfect thermalization limit $\theta\rightarrow\infty$. In both of the
plots, the differences, which are almost imperceptible in the panels a and b,
can be divided into two parts. The first part is the constant shift
proportional to $\Sigma^{2}$ in Eqs. (42) and (44). The remaining nontrivial
differences (NTD), which are plotted in panels c (for $\langle\delta
W^{2}\rangle$) and d (for $\langle\delta W\delta Q\rangle$), result from the
terms proportional to $\cos(2(\theta+\varphi))$ in the Eqs. (40) and (42). The
engine parameters are as follows: The initial half energy gap
$\epsilon_{c}=1$, the final half energy gap $\epsilon_{h}=3.7$, and the
Gaussian width of the measurement is $\Sigma=0.2$. The work strokes are
characterized by $\alpha=0.1$ and $\varphi=0$. The imperfect thermalization
strokes follow a master equation dynamics which asymptotically approaches
Gibbs states at the temperatures $\beta_{c}=0.25$, $\beta_{h}=0.025$ at the
relaxation rate $\gamma=0.025$. The engine is started in the invariant state
$\rho^{*}$ [Eq. (45)] of the unobserved cycle dynamics given by the Eq. (74).

The moments of heat can be obtained by setting $\epsilon_{c}=0$ and
multiplying appropriate global signs (negative for odd order moments and
positive for even ones) in the expressions for the moments of work. As seen
from the above equations, this also implies that the average heat for two-
pointer repeated contacts and repeated measurements are the same, i.e.,
$\langle Q\rangle=\langle Q\rangle^{\mathrm{RM}}=\langle
Q\rangle^{\mathrm{RC}}$. The variances of the heat are not displayed. They
behave in a similar way as the those of the work as presented in Fig. 2 with
solely a constant difference of $3\Sigma^{2}$ between repeated measurements
and repeated contacts and therefore are not extra displayed. Note that the
engine may starts at an arbitrary initial state specified by the parameters
$d$ and $q$. Here and in the following examples, however, we will initiate the
monitored sequence of $N$ cycles with the invariant state $\rho^{*}$ of the
map $\mathcal{F}$, i.e.

$\rho^{*}=\mathcal{F}(\rho^{*})$ (45)

where

$\mathcal{F}(\rho)=\Phi_{c}(\tilde{U}\Phi_{h}(U\rho
U^{\dagger})\tilde{U}^{\dagger})$ (46)

describes the dynamics of a complete cycle in the absence of any monitoring.

Furthermore we note that the master equation (38), which governs the
thermalization strokes of the engine, leads to decoupled equations of motion
for the populations and coherences. Hence, the thermalization maps $\Phi_{u}$
resulting as the solutions of the master equation (38) obey the conditions

$\mathrm{Tr}[L^{\pm}_{u}\Phi_{u}(\mathcal{P}^{k}_{u})]=\mathrm{Tr}[\mathcal{P}^{k}_{u}\Phi_{u}(L^{\pm}_{u})]=\mathrm{Tr}[L^{\pm}_{u}\Phi_{u}(L^{\mp}_{u})]=0,$
(47)

with swap-operators $L^{\pm}_{u}=|\pm_{u}\rangle\langle\mp_{u}|$, as defined
below Eq. (38), and $\mathcal{P}^{k}_{u}=|k_{u}\rangle\langle k_{u}|$ with
$k=\pm$, denoting the projection operators onto the eigenstates of the
Hamiltonian $H_{u}$. Due to this decoupling, the work PDFs for one- and two-
pointer repeated contacts become identical and we denote the average value as
$\langle W\rangle^{\mathrm{RC}}$ (see Append. B for more details). In the
limit $\theta\rightarrow\infty$, $\Sigma=0$, and $\rho=\rho_{c}$, which leads
to perfect thermalization, Eqs. (39)–(44) match Eqs. (35)–(37). For a highly
non-adiabatic work stroke ($\alpha$ large) the cosine terms in the expressions
for the first two moments play a significant role. Then the sign of
$\cos(2\theta)$ determines which one of the two monitoring schemes leads to a
larger work on average for a single cycle. This affects also the efficiency of
the engine because the average heat is the same for either scheme.

##  IV Numerical results

For a working substance with a $d$-dimensional Hilbert space, the number of
possible realizations of $(\vec{m},\vec{m}^{\prime})$ increases exponentially
as $d^{8N}$, and so do the coefficients $D^{\vec{m},\vec{m}^{\prime}}$. Hence,
the computational resources (CPU time and RAM) required to track the
individual coefficients become inaccessible even for $d=2$ already at a
moderately large number of cycles. To avoid this numerical impossibility, we
device a scheme to reduce the computational complexity to $N^{2}$ by grouping
terms according to their work values. The scheme (see Append. C for details)
is restricted to a two-level working substance undergoing thermalization
strokes that are governed by decoupled population and coherence dynamics as
specified by Eq. (47). Under these conditions the method is exact and allows
us to obtain the work distributions for a large number of cycles such that the
asymptote can be reached.

###  IV.1 One cycle, perfect thermalization

Figure 3: Panels a–b display the probability density functions (PDFs) of work
for adiabatic work strokes ($\alpha=0$) and different measurement schemes (a)
as well as for different non-adiabatic work strokes (b). Solid lines represent
two-pointer repeated contacts PDFs, and the dashed line presents the PDF for
repeated measurements. One-pointer results are not drawn since they are
exactly the same as two-pointer scheme for $\alpha=0$ and have a negligible
difference for $\alpha\neq 0$. Panel c displays the joint PDF of work and heat
calculated for repeated measurements, and d is the joint PDF calculated for
two-pointer repeated contacts. The PDFs are computed for perfect
thermalization with thermal baths at inverse temperatures $\beta_{c}=0.25$,
$\beta_{h}=0.025$, both characterized by the the same Ohmic spectal density
with Lorentz-Drude cut-off frequency $\omega_{D}=0.2$, system-operator
$S=\sigma_{z}+\sigma_{x}$, and damping coefficient $\gamma=0.5$. This results
in density matrices with an excited-state probability $d_{c}=0.37759$ and
$d_{h}=0.45388$, and coherence $q_{u}=\langle+_{u}|\rho_{u}|-_{u}\rangle$
given by $q_{h}=-1.9205\times 10^{-6}$ and $q_{c}=5.0813\times 10^{-5}$ (see
Append. A for details). The phase during the non-adiabatic evolution is set to
$\varphi=0$ in all cases, and the joint PDF is shown for $\alpha=0.1$. The
engine is started in the generalized Gibbs state of the cold bath with
parameters $d_{c}$ and $q_{c}$. All other parameters are the same as in Fig.
2.

Figure 3 presents the joint and marginal work and heat PDFs of an engine with
perfect thermalization for a single cycle ($N=1$). Panel a displays the work
PDFs for repeated measurements and repeated contacts for adiabatic work
strokes, i.e., for $\alpha=0$. In this case, there is no difference between
the marginal two-pointer work PDF given by Eq. (25) and the one-pointer work
PDF (26) as explained in Append. D. In accordance with the Eq. (26) the PDF
consists of three peaks located at $W=0$ and $W=\pm
2(\epsilon_{c}-\epsilon_{h})$ because transitions between different energy
levels are impossible for adiabatic work strokes. These peaks are broader for
repeated measurements compared to repeated contacts. This feature prevails for
non-adiabatic work strokes. In panel b we present the results for moderate and
strong non-adiabatic work strokes. In order not to clutter the plot we present
only the work PDFs resulting from repeated contacts. With decreasing
adiabaticity (increasing $\alpha$) six further peaks appear in accordance with
Eq. (26) accompanied by a decrease of the weights of the three adiabatic
peaks. At the same time the weights of those peaks at positive work values may
increase to such an extent that the average work becomes positive. That means
that the engine only consumes energy rather than to perform work. As for
classical engines any non-adiabaticity is detrimental with respect to the
efficiency but must be accepted to achieve finite power. The same features as
described for repeated contacts also hold for the repeated measurement PDFs
with the only difference that their peaks are broader.

Figures 3c and d present the joint work and heat PDFs for repeated
measurements and contacts, respectively. In the case of repeated measurements,
each individual Gaussian contribution contains a negative covariance
indicating an anti-correlation of the according work and heat contributions.
In contrast, for repeated contacts the individual Gaussian contributions are
isotropic. In both cases the main contribution to the covariance of work and
heat results from the linear superposition of the individual Gaussian
distributions.

###  IV.2 $N$ cycles, imperfect thermalization

Figure 4: Panel a presents a plot of the work PDFs for repeated measurements
(red) and repeated contacts (black) for $N=5$ cycles and panel b does so for
$N=100$. Two-pointer and one-pointer repeated contacts give exactly the same
PDFs of work due to the dynamical decoupling between populations and
coherences during the thermalization process (see text for details). While for
repeated contacts separate sharp lines with $N$-independent widths remain
visible also for large numbers of cycles, the according terms for repeated
measurements contribute with a width $\propto\sqrt{N}$ and hence merge in a
broad, practically Gaussian PDF. Panel c and d display the average work output
per cycle and the engine’s reliability, respectively. The horizontal black
dashed lines in c correspond to the asymptotic values of the average work per
cycle obtained from the analytic maps described in Append. B. The reliability
scales as the square root of the cycle number $N$ for large $N$. The repeated
measurement scheme has a detrimental impact on the engine in terms of both
average work output and reliability. The work strokes are characterized by
non-adiabatic parameters $\alpha=0.05$, $\varphi=0$. The energy gaps are
chosen as $\epsilon_{c}=1$ and $\epsilon_{h}=3.7$. The bath parameters are
$\beta_{c}=0.25$, $\beta_{h}=0.025$, the coupling constant $\gamma=0.025$, and
the thermalization times is taken such that
$\theta=\epsilon_{c}\tau_{c}=\epsilon_{h}\tau_{h}=8$. All pointer states
initially are Gaussian with a width of $\Sigma=0.2$. The engine is started in
the according invariant state $\rho^{*}$ [Eq. (45)] of the unobserved cycle
dynamics given by the Eq. (74).

Numerical results for adiabatic and non-adiabatic work strokes and imperfect
thermalization are presented in Fig. 4. The panels a and b display the work
PDFs for repeated measurements and repeated contacts for two different numbers
of cycles $N$. For the smaller number of cycles, $N=5$, still several peaks
are visible in the repeated measurements work PDF. Because their widths grow
as $\sqrt{N}$ this multiply peaked structure disappears with increasing number
of cycles to finally approach a Gaussian shape, which in the present case is
virtually indistinguishable from the exact PDF for $N=30$. In contrast, the
fine-structure of the repeated contact work PDF remains visible also for large
numbers of cycles due to the fact that the widths of the individual Gaussian
contributions are solely determined by the initial pointer state and hence are
$N$ independent.

In Fig. 4c the average work per cycle is displayed as it results as a function
of the number of cycles both for repeated measurements (red closed circles)
and repeated contacts (blue crosses). The initial state of the working
substance is chosen as the asymptotic state of an unobserved engine reached
after sufficiently many cycles. Only for a single cycle ($N=1$) the average
work for repeated measurements is larger than that for repeated contacts. For
any larger number of cycles the repeated contact scenario outperforms the
repeated measurements both with respect to the average work (see panel c) and
the reliability (see panel d). As a function of the number of cycles, the work
behaves like a biased random walk in that its first and second moments
asymptotically grow in proportion to the number of cycles yielding a constant
work per cycle and a reliability that increases as the square root of the
number of cycles.

For repeated measurements a map can be constructed that assigns the density
matrix of the working substance after a cycle to any initial density matrix.
The map possesses a unique invariant state describing the asymptotic regime of
a large number of cycles. This state can be used to calculate the asymptotic
average work per cycle, see Append. B. Even though such a map of the reduced
working substance density matrix does not exist for repeated contacts, in the
special cases of the Lindblad equation (38) that dynamically decouples
populations and coherences, the asymptotic value of the average work can be
determined from the asymptotic density matrix in the absence of any contacts
or measurements, see Append. B. This unobserved asymptotic density matrix of
the working substance can again be obtained as the invariant state of a
quantum map in the absence of contacts as considered in Ref. [54]. The
resulting average works are presented in panel c as dashed black lines.

###  IV.3 Power of a Landau-Zener engine

So far we have characterized the two work strokes by the Hamiltonians of the
expanded and the compressed working substances [see Eqs. (32) and (33)] and
unitary operators specifying the changes of its state, however, without
specifying a particular forcing protocol. As a consequence, the duration of
the work strokes is not specified. In general, determining the unitary time
evolution operator for a particular forcing protocol of given duration
presents a nontrivial task [51, 53]. In order to also have full control over
the duration of an engine cycle we consider a forcing protocol with a linear
driving as specified in Eq. (32). The unitary time evolution for a work stroke
of the duration $T_{1}/2$ is then determined by the parameters as given in the
Eq. (33). The duration of the thermalization strokes, $\tau_{h}$ and
$\tau_{c}$ are assumed to be fixed by a the single parameter
$\theta=\epsilon_{c}\tau_{c}=\epsilon_{h}\tau_{h}$, resulting in the total
thermalization time

$T_{2}=\tau_{c}+\tau_{h}=\theta(1/\epsilon_{h}+1/\epsilon_{c})$ (48)

Therefore, a single cycle takes the time $T_{1}+T_{2}$. The power of the
engine is given by

$P=-\frac{\langle W\rangle}{T_{1}+T_{2}}.$ (49)

where $\langle W\rangle$ denotes the average work performed on the working
substance per cycle.

Figure 5: Power per cycle (a-f) and average efficiency (g-l) as a function of
work and heat stroke durations $T_{1}$ and $T_{2}$ with cycle number $N=1$ (a,
d, g, and j), $N=20$ (b, e, h, and k), and $N=\infty$ (c, f, i, and l; see
Append. B). Panels a-c and g-i are for the engine measured via repeated
contacts whereas the other panels are for repeated measurements. The white
region in all panels correspond to the dud regime of the engine. The extended
dark red regions in g–l indicate regimes with efficiencies close to the ideal
Otto engine value $1-\epsilon_{c}/\epsilon_{h}\approx 0.73$. Fixed parameters
are $\epsilon_{c}=1$, $\epsilon_{h}=3.7$,
$\theta=\epsilon_{c}\tau_{c}=\epsilon_{h}\tau_{h}=8$, $\beta_{c}=0.25$,
$\beta_{h}=0.025$, $\gamma=0.025$, and $\Sigma=0.2$. The engine is always
started in the invariant state $\rho^{*}$ [Eq. (45)] of the according
unobserved cycle dynamics given by the Eq. (74).

In Fig. 5 the power and the efficiency for repeated contacts and repeated
measurements as well as for different numbers of cycles are presented as
functions of the duration of the work and thermalization strokes. The variance
of the pointer states is chosen relatively small such that coherences that may
build up during the non-adiabatic work strokes and also partially during an
imperfect thermalization, are almost completely suppressed by repeated
measurements while they can partly survive in the presence of repeated
contacts. Because the generation of coherences is further suppressed in the
limit of adiabatic work processes as well as for perfect thermalization, one
finds substantial differences between repeated measurements and repeated
contacts for small values of $T_{1}$ and small to intermediate values of
$T_{2}$. As already observed in Sec. IV.1, the total average work performed on
the working substance in a single cycle becomes positive for perfect
thermalization and work strokes with a sufficiently large parameter $\alpha$
and hence for a sufficiently short duration of the work strokes. This feature
of a dud engine with formally negative power and efficiency persists for
imperfect thermalization and also for any numbers of cycles at sufficiently
small values of the work stroke duration $T_{1}$ indicated as white regions in
Fig. 5. The border between the regions of a properly working and a faulty
engine displays well pronounced oscillations as a function of $T_{1}$ for any
number of cycles in the case of repeated contacts. The fact that these
oscillations are visible in the case of repeated measurements only for small
number of cycles indicates that they are caused by coherences which are
suppressed by repeated measurements.

Figure 6: Maximum power $P^{*}$, with respect to durations $T_{1}$ and
$T_{2}$, of the engine as a function of cycle number $N$ for both repeated
contacts and repeated measurements, blue and red symbols, respectively (panel
a). The dashed lines represent the asymptotic values of maximum power of the
engines in their periodic state asymptotically reached after many cycles, see
Append. B. The deviations of the engine from this limiting behavior decay
exponentially as the cycle number increases. Panel b represents the slowest
decaying factor $\Lambda_{2}$ (see B last paragraph for definition) in engines
monitored by repeated contacts, while c is the corresponding plot of engines
being repeatedly measured. The initial condition for the engine in panel a is
again given by the according invariant state $\rho^{*}$ as in the previous
figures. The other parameters $\epsilon_{c}=1$, $\epsilon_{h}=3.7$,
$\theta=\epsilon_{c}\tau_{c}=\epsilon_{h}\tau_{h}=8$, $\beta_{c}=0.25$,
$\beta_{h}=0.025$, and $\gamma=0.025$ are the same as in Fig. 5

The dependence of the maximal power $P^{*}=\max_{T_{1},T_{2}}P$ on the number
of cycles for repeated measurements and contacts is displayed in Fig. 6a. This
maximum is chosen from the set of all engines with variable stroke times
$T_{1}$ and $T_{2}$ and with all other parameters fixed. In both cases, the
power approaches an asymptotic value at large numbers of cycles. While the
asymptotic limit is reached from above already after a few cycles for repeated
measurements, it is approached considerably slower from below for repeated
contacts. The asymptotic values are indicated as broken horizontal lines. In
both cases, these values can be determined from the invariant maps that
propagate the density matrix cycle by cycle. For the repeated contact
scenario, the map describes the dynamics in the absence of any monitoring,
while for the repeated measurement scenario it describes the single-cycle
engine dynamics in the presence of non-selective energy measurements at the
times indicated in Fig. 1. For a more detailed description we refer to Append.
B.

For the considered parameter values of the engine, the position where the
maximum power is assumed in the $T_{1}$-$T_{2}$ plane turns out to be almost
independent of the number of cycles. Therefore, the convergence speed of the
map as it is iterated, mainly determines the convergence of the power with
increasing number of cycles as it can be seen in Fig. 6a. For sufficiently
many cycles the convergence speed is ruled by the second largest eigenvalue of
the respective map, while its largest eigenvalue is given by $\Lambda_{1}=1$.
It is non-degenerate for all engines with a positive thermalization time
$T_{2}>0$. The dependence of the second largest eigenvalue $\Lambda_{2}$ as a
function of $T_{1}$ and $T_{2}$ is exemplified in the panels Fig. 6b,c. In
both scenarios the second eigenvalue converges for large thermalization times
to zero, however much faster for repeated measurements than for repeated
contacts in agreement with the faster convergence of the maximal power towards
its asymptotic limit, i.e., $\Lambda^{RM}_{2}<\Lambda^{RC}_{2}$. Outside the
regime of short work strokes, the second eigenvalue is almost independent of
$T_{1}$ as it is primarily determined by the duration and strength of the
thermalization strokes.

##  V Conclusions

In this paper, we compared two different monitoring schemes, both designed to
record the statistics of work performed on and heat supplied to a quantum Otto
engine over a specified number of engine cycles. Based on these statistics, we
investigated the performance of such engines running subject to different
modes of operation, such as different work strokes ranging from adiabatic to
sudden compression and expansion of the working substance, and different types
of heating and cooling strokes. The method was formulated for general working
substances and exemplified for a two-level system acting as working substance.

Apart from the joint and marginal work and heat PDFs, we studied performance
measures such as efficiency, power and reliability that are all determined by
the first two moments of work and heat. For an engine that performs $N$
cycles, the two monitoring schemes are based on contacts of the working
substance with one, two, or $4N$ pointers, depending on the particular scheme.
These engine-pointer contacts are supposed to be so short that the dynamics of
the working substance is frozen during these moments but effective enough to
induce a change of the pointer state depending on the actual energy of the
system.

For the repeated measurement scheme each contact is immediately followed by a
projective measurement of the pointer state whose value corresponds to the
energy of the working substance at the instant of the contact, hence
representing a generalized energy measurement according to the von Neumann
scheme. The precision of the measurement can be controlled by the position
variance $\Sigma^{2}$ of the initial pointer states which we supposed to be
Gaussian and identical for all employed pointers. From the joint statistics of
all measured energy values, the joint and marginal PDFs of work and heat can
be found. These PDFs result as linear superpositions of Gaussian contributions
with centers at work and heat values given by proper sums and differences of
the $4N$ energy values but also fictitious values resulting from a double set
of energies. The latter though are exponentially suppressed provided the
variance $\Sigma^{2}$ is sufficiently small. All contributing Gaussians are of
the same shape with a variance growing proportionally to the number of cycles
leading to a smeared-out, almost Gaussian PDF for sufficiently large numbers
of cycles.

For repeated contacts there are just two pointers one being associated to work
and the other one to heat. These pointers are “kept on” during the whole
monitored period of cycles. At each contact the pointer state is modified
according to the present energetic state of the working substance, left
unchanged during the strokes and read out by projective position measurements
only after the completion of $N$ cycles. The energetic states at the
subsequent contacts are inscribed with alternating signs in the work pointer.
For the heat pointer only those contacts before and after a heat stroke are
registered with appropriate signs. It is interesting to note that the marginal
distributions of work obtained from the joint work and heat PDF in general
differ from the work PDF of a repeated contact scheme with only the work
pointer but without the heat pointer. The same holds for the heat PDFs
obtained either as the marginal of the joint work and heat PDF or as the
outcome of a single heat pointer repeated contact scheme. This aspect of
quantum contextuality is absent in the repeated measurement scheme. It
indicates a higher sensitivity of the repeated contact scheme to the
quantumness of a process as compared to the repeated measurement scheme.

The higher sensitivity of the repeated contact scheme is closely related to
the lesser back-action of the sole contacts as compared to the individual
read-out of the pointer-state position for repeated measurements leading to
different suppression factors making long-lived coherences possible. In the
example of an Otto engine with a two-level system working substance and
imperfect thermalization strokes this manifests itself in an average work per
cycle that first increases with the number of cycles and finally saturates at
a constant value for repeated contacts. In contrast, for repeated measurements
the average work per cycle decreases towards an asymptotic value that is
smaller than that for repeated contacts. Accordingly, the power but also the
efficiency as well as the reliability of the work, defined as the ratio of the
average work and its root-mean-square deviation, are larger for repeated
contacts than for repeated measurements provided that the engine protocol
allows coherences being still present after a cycle. Due to the less strong
back-action imposed by repeated contacts in comparison to repeated
measurements these coherences may subsist leading to a better performance of a
repeated-contacts-monitored engine. We expect that this behaviour is not
restricted to the special case of a two-level working substance but that it
prevails in general.

In a recent work, Watanabe et. al. [55] studied the amount of energy a quantum
Otto engine may deposit into an external system over several engine cycles.
The comparison of the amount of out-coupled energy into a “flywheel” with the
energy produced by the free engine according to a repeated contact diagnostic
scheme may provide further insight into the effectiveness of charging and
discharging strategies of quantum batteries.

## Acknowledgements

This research was supported by the Institute for Basic Science in Korea
(IBS-R024-Y2).

##  Appendix A Generalized Gibbs state

The generalized Gibbs state can be obtained up to second order in system-bath
interaction using canonical perturbation theory [42]. In this Appendix we omit
the subscripts $c$ or $h$ used in the main text for the cold and hot reservoir
for notational simplicity. The generalized Gibbs state

$\rho=\frac{\mathrm{Tr}_{B}\left[e^{-\beta
H_{\mathrm{tot}}}\right]}{\mathrm{Tr}\left[e^{-\beta
H_{\mathrm{tot}}}\right]},$ (50)

can be interpreted as the reduced density matrix of an extended system
governed by the total Hamiltonian $H_{\mathrm{tot}}=H_{S}+H_{B}+V$ in a
canonical state at the inverse temperature $\beta$. Here, $H_{S}$ and $H_{B}$,
denote the Hamiltonians of the proper system and the bath, respectively and
$V$ denotes the interaction. We assume that the bath is composed of
independent Harmonic oscillators of frequency $\omega_{n}$ and the system-bath
interaction takes the form $V=-S\otimes\sum_{n}c_{n}x_{n}$, with $S$ being a
system operator and $x_{n}$ being the position operator of the $n$th mode of
the bath coupling to the system via strength $c_{n}$. Performing a
perturbative expansion upto second-order in $V$ we obtain,

$\rho=\frac{e^{-\beta H_{S}}}{Z_{S}}+\frac{D}{Z_{S}}-\frac{e^{-\beta
H_{S}}\mathrm{Tr}[D]}{Z_{S}^{2}}.$ (51)

Above $Z_{S}=\mathrm{Tr}[\exp[-\beta H_{S}]]$ is the system partition function
of the bare system and the operator $D$ encodes information about finite
system-bath coupling, i.e.,

$\displaystyle D$ $\displaystyle=$ $\displaystyle e^{-\beta
H_{S}}\bigg{[}\int_{0}^{\beta}\,d\beta_{1}\int_{0}^{\beta_{1}}\,d\beta_{2}\,\tilde{S}_{\beta_{1}}\tilde{S}_{\beta_{2}}C(-i(\beta_{1}-\beta_{2}))\bigg{]}.$

The operator $\tilde{S}_{\mathrm{x}}=e^{\mathrm{x}H_{S}}Se^{-\mathrm{x}H_{S}}$
is the imaginary time free evolution of the system operator and the force-
force correlator
$C(\tau)=\sum_{n,n^{\prime}}c_{n}c_{n^{\prime}}\mathrm{Tr}[e^{-\beta
H_{B}}\tilde{x}_{n}(\tau)x_{n^{\prime}}]/\mathrm{Tr}[e^{-\beta H_{B}}]$ with
$\tilde{x}_{n}(\tau)=e^{i\tau H_{B}}x_{n}e^{-i\tau H_{B}}$. The force-force
correlator can be expressed in terms of the spectral density of the bath,

$J(\omega)=\pi\sum_{n}\frac{c_{n}^{2}}{2m_{n}\omega_{n}}\delta(\omega-\omega_{n}),$
(52)

as

$\displaystyle C(\tau)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\,d\omega\frac{J(\omega)}{\omega}\left[\coth\left(\frac{\beta\omega}{2}\right)\cos(\omega\tau)-i\sin(\omega\tau)\right].$

For our perfectly thermalizing two-level Otto engine discussed in Sec. III we
use the system operator $S=\sigma_{z}+\sigma_{x}$ such that the density matrix

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle
d|+\rangle\langle+|+(1-d)|-\rangle\langle-|+[q|+\rangle\langle-|+\mathrm{h.c.}],$
(53)

with

$\displaystyle d$ $\displaystyle=$
$\displaystyle\frac{e^{-\beta\epsilon}}{Z_{S}}\left(1+4\int_{0}^{\beta}\,d\lambda\frac{(\beta-\lambda)\sinh(2\epsilon\lambda)}{(1+e^{-2\beta\epsilon})}C(-i\lambda)\right),$
$\displaystyle q$ $\displaystyle=$
$\displaystyle\frac{e^{-\beta\epsilon}}{Z_{S}}\int_{0}^{\beta}\,d\lambda\frac{e^{2\beta\epsilon}(e^{-2\epsilon\lambda}-1)+(e^{2\epsilon\lambda}-1)}{\epsilon}C(-i\lambda),$

with $Z_{S}=e^{\beta\epsilon}+e^{-\beta\epsilon}$. Throughout this work we
have used an Ohmic spectral density with Lorentz-Drude cut-off as
$J(\omega)=\gamma\omega/(1+\omega^{2}/\omega_{D}^{2})$ for all our numerical
simulations.

##  Appendix B Reduced maps for measured engines

In this section, we present the methods to determine the map for the average
work performed in the $(N+1)$th cycle based on the knowledge of the average
work performed during $N$ cycles. In case of repeated measurements, the map
fully captures the distribution of work performed at any $(N+1)$th cycle,
whereas for repeated contacts we introduce a _fictitious_ work distribution
whose average is equivalent to the asymptotic average work per cycle. The
method is explained for work performed in the last cycle but can likewise be
applied for the heat absorbed in the last cycle.

###  B.1 Repeated measurements

First, we consider the case of repeated measurements. As described in Sec.
II.2, in this case after every stroke a new Gaussian pointer is adopted,
brought into contact with the engine, and measured. We begin with the non-
normalised density matrix for $(N+1)$ cycles, using Eq. (6),

$\phi^{\mathrm{RM}}_{\vec{x}_{N+1}}(\rho)=\sum_{\vec{m}_{N+1},\vec{m}^{\prime}_{N+1}}D^{\vec{m}_{N+1},\vec{m}^{\prime}_{N+1}}(\rho)\prod_{k=1}^{4(N+1)}\sigma\left(x_{k}-e_{k}^{m_{k}},x_{k}-e_{k}^{m_{k}^{\prime}}\right).$
(55)

Utilising the relation in Eq. (7), splitting the product above into one for
$N$ cycles and the other for $(N+1)$th cycle, and using the non-normalized
density matrix for $N$ cycles [Eq. (6)] we can rewrite the above equation as,

$\displaystyle\phi_{\vec{x}_{N+1}}^{\mathrm{RM}}(\rho)$ $\displaystyle=$
$\displaystyle\sum_{\mathbf{m},\mathbf{m}^{\prime}}S^{\mathbf{m},\mathbf{m}^{\prime}}\left(\phi_{\vec{x}_{N}}^{\mathrm{RM}}(\rho)\right)$
(56)
$\displaystyle\times\prod_{k=4N+1}^{4(N+1)}\sigma\left(x_{k}-e_{k}^{m_{k}},x_{k}-e_{k}^{m_{k}^{\prime}}\right).$

The work performed in the $(N+1)$th cycle is given by

$W^{(N+1)}=\sum_{k=4N+1}^{4(N+1)}(-1)^{k}e_{k},$ (57)

and its corresponding PDF is

$\displaystyle P^{\mathrm{RM}}(W^{(N+1)})$ $\displaystyle=$ $\displaystyle\int
d^{4(N+1)}x\,\delta\left(W^{(N+1)}-\sum_{k=4N+1}^{4(N+1)}(-1)^{k}x_{k}\right)p(\vec{x}_{N+1}),$
(58)

with $p(\vec{x}_{N+1})=\mathrm{Tr}[\phi_{\vec{x}_{N+1}}^{\mathrm{RM}}(\rho)]$
according to Eq. (9). Noting that the integration variables within the delta
function depend only on the $(N+1)$th cycle pointer positions $x_{k}$, the
integral can be split as,

$P^{\mathrm{RM}}(W^{(N+1)})=\int
d\mathbf{x}\,\delta\left(W^{(N+1)}-\sum_{k=4N+1}^{4(N+1)}(-1)^{k}x_{k}\right)\int
d^{4N}x\,p(\vec{x}_{N+1}),$ (59)

with $\mathbf{x}$ defined by $\vec{x}_{N+1}=\vec{x}_{N}\oplus\mathbf{x}$. All
information about the previous $N$ cycles is contained in the rightmost
integral above. Hence, using Eq. (56) this integral can be expressed as,

$\displaystyle\int d^{4N}x\,p(\vec{x}_{N+1})$ $\displaystyle=$
$\displaystyle\sum_{\mathbf{m},\mathbf{m}^{\prime}}\mathrm{Tr}\left[S^{\mathbf{m},\mathbf{m}^{\prime}}\left(\int
d^{4N}x\,\phi_{\vec{x}_{N}}^{\mathrm{RM}}(\rho)\right)\right]$
$\displaystyle\times\prod_{k=4N+1}^{4(N+1)}\sigma\left(x_{k}-e_{k}^{m_{k}},x_{k}-e_{k}^{m_{k}^{\prime}}\right).$

The simplification on the r.h.s. implies that as long as we know the map
$S^{\mathbf{m},\mathbf{m}^{\prime}}$ and the state it acts on, i.e., $\int
d^{4N}x\,\phi_{\vec{x}_{N}}^{\mathrm{RM}}(\rho)$, we can evaluate the work
distribution $P^{\mathrm{RM}}(W^{(N+1)})$. Thus, the task reduces to finding
the map $F^{\mathrm{RM}}$, such that,

$\int
d^{4N}x\,\phi_{\vec{x}_{N}}^{\mathrm{RM}}(\rho)=F^{\mathrm{RM}}\left(\int
d^{4N-4}x\,\phi_{\vec{x}_{N-1}}^{\mathrm{RM}}(\rho)\right).$ (60)

Using Eq. (56) with $N\rightarrow N-1$ and integrating both sides we can
easily obtain the map,

$F^{\mathrm{RM}}(\varrho)=\sum_{\mathbf{\tilde{m}},\mathbf{\tilde{m}}^{\prime}}S^{\mathbf{\tilde{m}},\mathbf{\tilde{m}}^{\prime}}(\varrho)\prod_{k=4N-3}^{4N}e^{-\frac{1}{8\Sigma^{2}}\left(e_{k}^{m_{k}}-e_{k}^{m_{k}^{\prime}}\right)^{2}},$
(61)

where $\varrho=\int d^{4N-4}x\,\phi_{\vec{x}_{N-1}}^{\mathrm{RM}}(\rho)$,
${\bf\tilde{m}^{(\prime)}}$ are defined by the sequence
$\vec{m}^{(\prime)}_{N}=\vec{m}^{(\prime)}_{N-1}\oplus{\bf\tilde{m}^{(\prime)}}$.
Note the difference between ${\bf\tilde{m}^{(\prime)}}$ and ${\bf
m^{(\prime)}}$ is that $N$ is replaced by $(N+1)$ in the concatenation. The
map $F^{\mathrm{RM}}$ is _Markovian_ since it gives us the state at the $N$th
cycle given the state at $N-1$ [see Eq. (60)]. Overall, in case of repeated
measurements since the map $F^{\mathrm{RM}}$ is known it is straightforward to
obtain the work PDF for the work performed in the $(N+1)$th cycle using Eqs.
(59)–(61).

###  B.2 Repeated contacts

In case of repeated contacts, one cannot assign a value of work to each
individual cycle as we did with $W^{(j)}$ for repeated measurements by using
the classically stored information about the outcomes of the sequence of
energy measurements for $N$ cycles. To circumvent this problem we restrict our
consideration to engines with a working substance consisting of a two-level
system and moreover being subject to thermalization strokes obeying the
decoupled population and coherences dynamics specified by the Eq. (47). In
this case we are able to introduce a fictitious work PDF that can be computed
recursively and that determines the average work per cycle in the limit of
large cycle numbers when the engine has reached a steady mode of operation.

In order to proceed, we begin with the work PDFs defined in Eqs. (25) and
(26),

$\displaystyle P^{\mathrm{RC,2P}}(W)$ $\displaystyle=$
$\displaystyle\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{W}^{\vec{m}}-\mathcal{W}^{\vec{m}\prime})^{2}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{Q}^{\vec{m}}-\mathcal{Q}^{\vec{m}\prime})^{2}}g_{\Sigma^{2}}(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}}),$
$\displaystyle P^{\mathrm{RC,1P}}(W)$ $\displaystyle=$
$\displaystyle\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{W}^{\vec{m}}-\mathcal{W}^{\vec{m}\prime})^{2}}g_{\Sigma^{2}}(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}}).$
(62)

In case of dynamical decoupling between populations and coherences for a two-
level working substance, [Eq. (47)], only the terms with

$e_{2i}^{m_{2i}}-e_{2i}^{m_{2i}^{\prime}}=e_{2i+1}^{m_{2i+1}}-e_{2i+1}^{m_{2i+1}^{\prime}}$
(63)

$\forall i=1,\cdots,2N-1$ contribute to the work distributions (see Append.
B.3 for details). This simplifies the exponents [see Eq. (22) and (23)]
appearing in the work distributions, namely,

$\displaystyle\mathcal{W}^{\vec{m}}-\mathcal{W}^{\vec{m}^{\prime}}$
$\displaystyle=$ $\displaystyle-
e_{1}^{m_{1}}+e_{1}^{m_{1}^{\prime}}+e_{4N}^{m_{4N}}-e_{4N}^{m_{4N}^{\prime}},$
(64) $\displaystyle\mathcal{Q}^{\vec{m}}-\mathcal{Q}^{\vec{m}\prime}$
$\displaystyle=$ $\displaystyle 0.$ (65)

Thus, for the models considered in Secs. IV.2–IV.3 the one- and two-pointer
work PDFs always coincide. Moreover, using the above simplifications one can
easily identify the non-normalized reduced density matrix for the working
substance as,

$\displaystyle\phi^{\mathrm{RC}}_{W}(\rho)$ $\displaystyle=$
$\displaystyle\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}(\rho)\,e^{-\frac{1}{8\Sigma^{2}}\left(e_{1}^{m_{1}}-e_{1}^{m_{1}^{\prime}}\right)^{2}}$
(66) $\displaystyle\times
g_{\Sigma^{2}}\left(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}}\right)$
$\displaystyle=$
$\displaystyle\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}(\tilde{\rho})g_{\Sigma^{2}}\left(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}}\right).$

Here we also omitted the only remaining contribution
$e_{4N}^{m_{4N}}-e_{4N}^{m_{4N}^{\prime}}$ to the exponential factor
$e^{-(\mathcal{W}^{\vec{m}}-\mathcal{W}^{\vec{m}^{\prime}})^{2}/(8\Sigma^{2})}$
because the terms $D^{\vec{m},\vec{m}^{\prime}}$ with $e_{4N}^{m_{4N}}\neq
e_{4N}^{m_{4N}^{\prime}}$ vanish in the final trace operation applied to
obtain the work PDFs from the non-normalized reduced density matrices. In
addition, we made use of the fact that the innermost part of the telescopic
action of the map $D^{\vec{m},\vec{m}^{\prime}}$ on $\rho$ is given by
$\mathcal{P}^{m_{1}}_{1}\rho\mathcal{P}^{m^{\prime}_{1}}_{1}$, see Eq. (7).
Splitting the density matrix $\rho$ into its diagonal part
$\rho_{d}=\sum_{m=1,2}\mathcal{P}^{m}_{1}\rho\mathcal{P}^{m}_{1}$ and the
remaining coherence contribution $\rho_{c}=\rho-\rho_{d}$ one obtains
$\exp\left[-(e_{1}^{m_{1}}-e_{1}^{m_{1}^{\prime}})^{2}/(8\Sigma^{2})\right]\mathcal{P}^{m}_{1}\rho\mathcal{P}^{m^{\prime}}_{1}=\mathcal{P}^{m}_{1}\tilde{\rho}\mathcal{P}^{m^{\prime}}_{1}$
where

$\tilde{\rho}\equiv\rho_{d}+e^{-\frac{\epsilon^{2}_{c}}{2\Sigma^{2}}}\rho_{s}.$
(67)

Note that the map $\rho\to\tilde{\rho}$ preserves the trace and the
positivity.

Next we introduce the fictitious work distribution

$\displaystyle\tilde{P}^{\mathrm{RC}}(W^{(N+1)})$ $\displaystyle=$
$\displaystyle\int
d\mathbf{x}\,\delta\left(W^{(N+1)}-\sum_{k=4N+1}^{4(N+1)}(-1)^{k}x_{k}\right)$
(68) $\displaystyle\times\int d^{4N}x\,\tilde{p}(\vec{x}_{N+1}),$

which exactly resembles Eq. (59) with,

$\displaystyle\tilde{p}(\vec{x}_{N})$ $\displaystyle=$
$\displaystyle\sum_{\vec{m}_{N},\vec{m}^{\prime}_{N}}\mathrm{Tr}\left[D^{\vec{m}_{N},\vec{m}^{\prime}_{N}}(\tilde{\rho})\right]$
(69)
$\displaystyle\times\prod_{k=1}^{4N}g_{\Sigma^{2}}\left(x_{k}-\frac{1}{2}\left(e^{m_{k}}_{k}+e^{m_{k}^{\prime}}_{k}\right)\right).$

The average work from the fictitious distribution,

$\langle{W}^{(N+1)}\rangle=\int
dW^{(N+1)}\,W^{(N+1)}\tilde{P}^{\mathrm{RC}}(W^{(N+1)}),$ (70)

represents the difference of the average work performed by the engine in the
$(N+1)$th and $N$th cycle, i.e., $\langle{W}^{(N+1)}\rangle=\langle
W\rangle_{N+1}-\langle W\rangle_{N}$ ( this relation follows from the
definitions of the individual quantities as explained below). As the monitored
engine reaches its asymptotic state, $\langle W^{(\infty)}\rangle$ represents
the average work performed by the engine per cycle, similarly to the repeated
measurements case, because the monitoring process on average has the same
effect on the engine state in the $N$th and $(N+1)$th cycle after sufficiently
many cycles.

The average work performed in $(N+1)$ and in $N$ cycles can be written as

$\displaystyle\langle W\rangle_{N+1}$ $\displaystyle=$ $\displaystyle\int
dW\,W\sum_{\vec{m}_{N+1},\vec{m}_{N+1}^{\prime}}\mathrm{Tr}\left[D^{\vec{m}_{N+1},\vec{m}_{N+1}^{\prime}}(\tilde{\rho})\right]g_{\Sigma^{2}}(W-\mathcal{W}^{\vec{m}_{N+1},\vec{m}_{N+1}^{\prime}})$
$\displaystyle=$ $\displaystyle\int
dW\,W\sum_{\mathbf{m},\mathbf{m}^{\prime}}\mathrm{Tr}\left[S^{\mathbf{m},\mathbf{m}^{\prime}}(\sum_{\vec{m}_{N},\vec{m}_{N}^{\prime}}D^{\vec{m}_{N},\vec{m}_{N}^{\prime}}(\tilde{\rho})])\right]g_{\Sigma^{2}}(W-\mathcal{W}^{\vec{m}_{N},\vec{m}_{N}^{\prime}}-\mathcal{W}^{\mathbf{m},\mathbf{m}^{\prime}}),$
$\displaystyle\langle W\rangle_{N}$ $\displaystyle=$ $\displaystyle\int
dW\,W\sum_{\mathbf{m},\mathbf{m}^{\prime}}\mathrm{Tr}\left[S^{\mathbf{m},\mathbf{m}^{\prime}}(\sum_{\vec{m}_{N},\vec{m}_{N}^{\prime}}D^{\vec{m}_{N},\vec{m}_{N}^{\prime}}(\tilde{\rho})])\right]g_{\Sigma^{2}}(W-\mathcal{W}^{\vec{m}_{N},\vec{m}_{N}^{\prime}}).$

The third line holds because the non-selective map $F^{\mathrm{RC}}$ that is
defined as

$F^{\mathrm{RC}}(\varrho)=\sum_{\mathbf{\tilde{m}},\mathbf{\tilde{m}}^{\prime}}S^{\mathbf{\tilde{m}},\mathbf{\tilde{m}}^{\prime}}(\varrho),$
(71)

preserves the trace so that

$\displaystyle\sum_{\mathbf{m},\mathbf{m}^{\prime}}\mathrm{Tr}\left[S^{\mathbf{m},\mathbf{m}^{\prime}}(\sum_{\vec{m}_{N},\vec{m}_{N}^{\prime}}D^{\vec{m}_{N},\vec{m}_{N}^{\prime}}(\tilde{\rho})])\right]$
$\displaystyle=$
$\displaystyle\sum_{\vec{m}_{N},\vec{m}_{N}^{\prime}}\mathrm{Tr}\left[D^{\vec{m}_{N},\vec{m}_{N}^{\prime}}(\tilde{\rho})]\right].$

Now the difference $\Delta W=\langle W\rangle_{N+1}-\langle
W\rangle_{N}$becomes

$\displaystyle\Delta W$ $\displaystyle=$
$\displaystyle\sum_{\mathbf{m},\mathbf{m}^{\prime}}\mathrm{Tr}\left[S^{\mathbf{m},\mathbf{m}^{\prime}}(\sum_{\vec{m}_{N},\vec{m}_{N}^{\prime}}D^{\vec{m}_{N},\vec{m}_{N}^{\prime}}(\tilde{\rho})])\right]$
$\displaystyle\int
dW\,W\left(g_{\Sigma^{2}}(W-\mathcal{W}^{\mathbf{m},\mathbf{m}^{\prime}})-g_{\Sigma^{2}}(W)\right)$
$\displaystyle=$
$\displaystyle\sum_{\mathbf{m},\mathbf{m}^{\prime}}\mathrm{Tr}\left[S^{\mathbf{m},\mathbf{m}^{\prime}}(\sum_{\vec{m}_{N},\vec{m}_{N}^{\prime}}D^{\vec{m}_{N},\vec{m}_{N}^{\prime}}(\tilde{\rho})])\right]\mathcal{W}^{\mathbf{m},\mathbf{m}^{\prime}},$

where we shifted the integrand by the amount of
$\mathcal{W}^{\vec{m},\vec{m}^{\prime}}$ and used the normalization of the
Gaussian distributions. On the other hand,

$\displaystyle\langle{W}^{(N+1)}\rangle$ $\displaystyle=$ $\displaystyle\int
dW^{(N+1)}\,W^{(N+1)}\tilde{P}^{\mathrm{RC}}(W^{(N+1)})$ $\displaystyle=$
$\displaystyle\int dW^{(N+1)}\,\int
d\mathbf{x}\,W^{(N+1)}\delta\left(W^{(N+1)}-\sum_{k=4N+1}^{4(N+1)}(-1)^{k}x_{k}\right)\sum_{\vec{m}_{N+1},\vec{m}_{N+1}^{\prime}}\mathrm{Tr}[D^{\vec{m}_{N+1},\vec{m}_{N+1}^{\prime}}(\tilde{\rho})]$
$\displaystyle\times\prod_{k=4N+1}^{4(N+1)}g_{\Sigma^{2}}\left(x_{k}-\frac{1}{2}(e_{k}^{m_{k}}+e_{k}^{m_{k}^{\prime}})\right).$

Integrating $\mathbf{x}$ first, we get

$\displaystyle\langle{W}^{(N+1)}\rangle$ $\displaystyle=$ $\displaystyle\int
dW^{(N+1)}\,W^{(N+1)}\sum_{\mathbf{m},\mathbf{m}^{\prime}}\mathrm{Tr}\left[S^{\mathbf{m},\mathbf{m}^{\prime}}(\sum_{\vec{m}_{N},\vec{m}_{N}^{\prime}}D^{\vec{m}_{N},\vec{m}_{N}^{\prime}}(\tilde{\rho})])\right]g(W^{(N+1)}-\mathcal{W}^{\mathbf{m},\mathbf{m}^{\prime}})$
$\displaystyle=$
$\displaystyle\sum_{\mathbf{m},\mathbf{m}^{\prime}}\mathrm{Tr}\left[S^{\mathbf{m},\mathbf{m}^{\prime}}\left(\sum_{\vec{m}_{N},\vec{m}_{N}^{\prime}}D^{\vec{m}_{N},\vec{m}_{N}^{\prime}}(\tilde{\rho})]\right)\right]\mathcal{W}^{\mathbf{m},\mathbf{m}^{\prime}}=\Delta
W.$

Furthermore, using a fictitious non-normalised density matrix, similar to Eq.
(55),

$\displaystyle\tilde{\phi}^{\mathrm{RC}}_{\vec{x}_{N+1}}(\tilde{\rho})$
$\displaystyle=$
$\displaystyle\sum_{\vec{m}_{N+1},\vec{m}_{N+1}^{\prime}}D^{\vec{m}_{N+1},\vec{m}_{N+1}^{\prime}}(\tilde{\rho})$
(72)
$\displaystyle\prod_{k=1}^{4(N+1)}\tilde{\sigma}\left(x_{k}-e_{k}^{m_{k}},x_{k}-e_{k}^{m_{k}^{\prime}}\right),$

with the fictitious pointer state

$\tilde{\sigma}\left(x_{k}-e_{k}^{m_{k}},x_{k}-e_{k}^{m_{k}^{\prime}}\right)=g_{\Sigma^{2}}\left(x_{k}-\frac{e_{k}^{m_{k}}+e_{k}^{m_{k}^{\prime}}}{2}\right),$
(73)

we can replicate the proof presented for repeated measurements in Append. B.1,
Eqs. (59)–(61) leading to the conclusion that the fictitious energy PDF
$\tilde{p}(\vec{x}^{N})$ [Eq. (69)] entering the fictitious work PDF
$\tilde{P}^{RC}(W^{(N+1)})$ [Eq. (68)] is determined by the invariant state
$\rho^{*}$ of the map $F^{RC}$, i.e. by the solution of the equation

$F^{RC}(\rho^{*})=\rho^{*},$ (74)

where $F^{RC}$ is defined in Eq. (71). Comparing the above equation with Eq.
(61) we notice the missing exponential terms in the map. This is simply
because the exponential terms in case of repeated contacts get simplified, see
Eq. (66), and they only affect the initial state of the engine. In other
words, in case of repeated measurements to obtain the state of the engine
after $N$ cycles we apply the map $F^{\mathrm{RM}}$ $N$ times on the initial
density matrix $\rho$ whereas in case of repeated contacts the map
$F^{\mathrm{RC}}$ is applied $N$ times on the fictitious initial density
matrix $\tilde{\rho}$.

As trace and positivity preserving maps, $F^{RC}$ and $F^{RM}$ possess at
least one invariant solution. The uniqueness of the invariant solution is
guaranteed by the dissipative nature of the thermalization strokes. As a
conseqence, the invariant densities can be iteratively reached as
$\rho^{*}=\lim_{N\to\infty}(F^{X})^{N}(\rho)$, where $X=RC,RM$. Accordingly,
the spectra of these maps contain one non-degenerate eigenvalue
$\Lambda_{1}=1$ and the absolute values of all other eigenvales are smaller
than one.

In Fig. 6b-c, the second-largest eigenvalues $\Lambda_{2}$ of the maps
$F^{\mathrm{RC(RM)}}$ are displayed as functions of the durations $T_{1}$ and
$T_{2}$. When $\Lambda_{2}$ is close to 1, the engine needs a long time to
converge to the asymptotic cycle and vice versa.

###  B.3 Simplification due to dynamical decoupling

In order to prove Eq. (63) we consider the coefficient,

$D^{\vec{m},\vec{m}^{\prime}}=\mathrm{Tr}[\mathcal{P}^{m_{4N}}_{4N}\cdots\mathcal{P}^{m_{2i+1}}_{2i+1}\Phi_{u}\left(\mathcal{P}^{m_{2i}}_{2i}\cdots\rho\cdots\mathcal{P}^{m_{2i}^{\prime}}_{2i}\right)\,\mathcal{P}^{m_{2i+1}^{\prime}}_{2i+1}\cdots\mathcal{P}^{m_{4N}^{\prime}}_{4N}],$
(75)

and focus on the terms involving $2i+1$ and $2i$ subscripts. Now,two
alternative cases exist:

  1. 1.

$m_{2i}=m_{2i}^{\prime}$ and

  2. 2.

$m_{2i}\neq m_{2i}^{\prime}$,

for which we study the central term in Eq. (75) within the trace, i.e.,

$\langle
m_{2i+1}|\Phi_{u}\left(\mathcal{P}^{m_{2i}}_{2i}\cdots\rho\cdots\mathcal{P}^{m_{2i}^{\prime}}_{2i}\right)|m_{2i+1}^{\prime}\rangle\propto
A,$ (76)

with $A=\langle m_{2i+1}|\Phi_{u}\left(|m_{2i}\rangle\langle
m_{2i}^{\prime}|\right)|m_{2i+1}^{\prime}\rangle$.

For case (i), we could have either $m_{2i+1}=m_{2i+1}^{\prime}$ or
$m_{2i+1}\neq m_{2i+1}^{\prime}$. For $m_{2i+1}=m_{2i+1}^{\prime}$, Eq. (63)
is satisfied. Whereas for $m_{2i+1}\neq m_{2i+1}^{\prime}$,
$A=\mathrm{Tr}[L^{m_{2i+1}}_{u}\Phi_{u}(\mathcal{P}^{m_{2i}}_{2i})]$ which is
zero due to Eq. (47). Hence, for case (i) only when Eq. (63) is satisfied we
have a contribution to the work distributions.

For case (ii), we divide into three sub-cases: (iia)
$m_{2i+1}=m_{2i+1}^{\prime}$ for which
$A=\mathrm{Tr}[\mathcal{P}^{m_{2i+1}}_{2i+1}\Phi_{u}(L_{u}^{m_{2i}})]$ is
again zero due to Eq. (47), (iib) $m_{2i}=m_{2i+1}^{\prime}$ and
$m_{2i}^{\prime}=m_{2i+1}$ for which
$A=\mathrm{Tr}[L_{u}^{\pm}\Phi_{u}(L_{u}^{\mp})]=0$, and (iic)
$m_{2i}=m_{2i+1}$ and $m_{2i}^{\prime}=m_{2i+1}^{\prime}$ for which Eq. (63)
is satisfied. All other cases like $m_{2i}=m_{2i+1}$ and $m_{2i}^{\prime}\neq
m_{2i+1}^{\prime}$ reduce to case (iia) since we deal with a two-level working
substance.

This exhausts all possible cases and overall only when Eq. (63) is satisfied
we get a contribution to the work distribution.

##  Appendix C Reduction in computational costs

In this section, we device a computational scheme to reduce the complexity
from exponential ($d^{8N}$ with $d$ being the dimension of working substance
Hilbert space and $N$ being the number of cycles) to quadratic ($N^{2}$) as
proposed in the main text. Our proof below focuses on the work PDFs given by
Eqs. (18) and (26),

$\displaystyle P^{\mathrm{RM}}(W)$ $\displaystyle=$
$\displaystyle\sum_{\vec{m}_{N},\vec{m}^{\prime}_{N}}D^{\vec{m}_{N},\vec{m}^{\prime}_{N}}e^{-\frac{1}{8\Sigma^{2}}{\sum_{k=1}^{4N}}(e^{m_{k}}_{k}-e^{m_{k}\prime}_{k})^{2}}g_{4N\Sigma^{2}}(W-\mathcal{W}^{\vec{m}_{N},\vec{m}^{\prime}_{N}}),$
$\displaystyle P^{\mathrm{RC}}(W)$ $\displaystyle=$
$\displaystyle\sum_{\vec{m}_{N},\vec{m}^{\prime}_{N}}D^{\vec{m}_{N},\vec{m}^{\prime}_{N}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{W}^{\vec{m}_{N}}-\mathcal{W}^{\vec{m}{\prime}_{N}})^{2}}g_{\Sigma^{2}}(W-\mathcal{W}^{\vec{m}_{N},\vec{m}^{\prime}_{N}}).$
(77)

Above we have reintroduced the subscript $N$ indicating the cycle number which
was suppressed in the main text. The scheme for the heat PDFs can be derived
in exactly the same way. The above equations can be rewritten as,

$P^{\mathrm{X}}(W)=\sum_{\mathcal{W}_{N}}D^{\mathrm{X}}_{\mathcal{W}_{N}}\>g_{\Sigma^{2}_{X}}(W-\mathcal{W}_{N}),$
(78)

where $\mathrm{X}=\mathrm{RM},\mathrm{RC}$ stands for repeated measurements or
repeated contacts. The variances of the superimposed Gaussians depend on the
type of monitoring with $\Sigma^{2}_{RM}=4N\Sigma^{2}$ and
$\Sigma^{2}_{RC}=\Sigma^{2}$. The sum runs over all possible values of
$\mathcal{W}_{N}=\mathcal{W}^{\vec{m}_{N},\vec{m}^{\prime}_{N}}$ and the
coefficients $D^{\mathrm{X}}_{\mathcal{W}_{N}}$ are given by

$D^{\mathrm{X}}_{\mathcal{W}_{N}}=\mathrm{Tr}[\mathcal{D}^{\mathrm{X}}_{\mathcal{W}_{N}}(\rho)].$
(79)

The operators $\mathcal{D}^{\mathrm{X}}_{\mathcal{W}_{N}}(\rho)$ are defined
as

$\displaystyle\mathcal{D}^{\mathrm{RM}}_{\mathcal{W}_{N}}(\rho)$
$\displaystyle=$
$\displaystyle\mkern-36.0mu\sum_{\stackrel{{\scriptstyle\vec{m}_{N},\vec{m}^{\prime}_{N}}}{{\mathcal{W}_{N}=\mathcal{W}^{\vec{m}_{N},\vec{m}^{\prime}_{N}}}}}\mkern-32.0muD^{\vec{m}_{N},\vec{m}^{\prime}_{N}}(\rho)e^{-\frac{1}{8\Sigma^{2}}\sum_{k=1}^{4N}\left(e^{m_{k}}_{k}-e^{m^{\prime}_{k}}_{k}\right)^{2}},$
$\displaystyle\mathcal{D}^{\mathrm{RC}}_{\mathcal{W}_{N}}(\rho)$
$\displaystyle=$
$\displaystyle\mkern-36.0mu\sum_{\stackrel{{\scriptstyle\vec{m}_{N},\vec{m}^{\prime}_{N}}}{{\mathcal{W}_{N}=\mathcal{W}^{\vec{m}_{N},\vec{m}^{\prime}_{N}}}}}\mkern-32.0muD^{\vec{m}_{N},\vec{m}^{\prime}_{N}}(\rho)e^{-\frac{1}{8\Sigma^{2}}\left(e^{m_{1}}_{1}-e^{m^{\prime}_{1}}_{1}\right)^{2}}.$
(80)

In case of repeated contacts, the general expression reduces to the above due
to the dynamical decoupling between populations and coherence [see Append. B.2
Eq. (66)]. The reduction in computational resources occurs because we only
need to compute the coefficients for a fixed $\mathcal{W}_{N}$ that scale as
$N^{2}$ rather than the exponential scaling of $d^{8N}$.

In order to see this reduction in computational resources, let us apply the
map $S^{\mathbf{m},\mathbf{m}^{\prime}}$, defined in Eq. (8), to both sides of
Eq. (80), giving,

$\displaystyle
S^{\mathbf{m},\mathbf{m}^{\prime}}\left(\mathcal{D}^{\mathrm{RC}}_{\mathcal{W}_{N}}(\rho)\right)$
$\displaystyle=$
$\displaystyle\mkern-36.0mu\sum_{\stackrel{{\scriptstyle\vec{m}_{N},\vec{m}^{\prime}_{N}}}{{\mathcal{W}_{N}=\mathcal{W}^{\vec{m}_{N},\vec{m}^{\prime}_{N}}}}}\mkern-32.0muS^{\mathbf{m},\mathbf{m}^{\prime}}\left(D^{\vec{m}_{N},\vec{m}^{\prime}_{N}}(\tilde{\rho})\right).$
(81)

Here we derive the scaling reduction for the case of repeated contacts and the
proof for repeated measurements can be obtained in the same spirit. The
exponential terms appearing in Eq. (80) have been absorbed in the effective
initial density matrix $\tilde{\rho}$ [see Eq. (67)] since they only depend on
the first contact ($m_{1}$ and $m_{1}^{\prime}$), hence, only alters the
initial state. It is important to note here that our proof works only if the
exponentials appearing in Eq. (77) can be expressed as a product of
exponentials implying that such a dramatic reduction is always possible for
repeated measurements and is impossible if the populations and coherence do
not dynamically decouple in the case of repeated contacts. Using Eq. (7) on
the r.h.s. we obtain,

$\displaystyle
S^{\mathbf{m},\mathbf{m}^{\prime}}\left(\mathcal{D}^{\mathrm{RC}}_{\mathcal{W}_{N}}(\rho)\right)$
$\displaystyle=$
$\displaystyle\mkern-36.0mu\sum_{\stackrel{{\scriptstyle\vec{m}_{N},\vec{m}^{\prime}_{N}}}{{\mathcal{W}_{N}=\mathcal{W}^{\vec{m}_{N},\vec{m}^{\prime}_{N}}}}}\mkern-32.0muD^{\vec{m}_{N+1},\vec{m}^{\prime}_{N+1}}(\tilde{\rho}).$
(82)

Recall that $\mathbf{m}^{(\prime)}$ arises from the sequence
$\vec{m}^{(\prime)}_{N+1}=\vec{m}^{(\prime)}\oplus\mathbf{m}^{(\prime)}$,
defined below Eq. (7). Summing both sides of Eq. (82) over
$\mathbf{m}^{(\prime)}$ and $\mathcal{W}_{N}$ (constrained by
$\mathcal{W}_{N+1}=\mathcal{W}^{\vec{m}_{N+1},\vec{m}^{\prime}_{N+1}}=\mathcal{W}_{N}+\mathcal{W}^{\mathbf{m},\mathbf{m}^{\prime}}$)
we obtain,

$\displaystyle\sum_{\stackrel{{\scriptstyle\mathbf{m},\mathbf{m}^{\prime},\mathcal{W}_{N}}}{{\mathcal{W}_{N+1}=\mathcal{W}_{N}+\mathcal{W}^{\mathbf{m},\mathbf{m}^{\prime}}}}}\mkern-42.0muS^{\mathbf{m},\mathbf{m}^{\prime}}\left(\mathcal{D}^{\mathrm{RC}}_{\mathcal{W}_{N}}(\rho)\right)$
$\displaystyle=$
$\displaystyle\mkern-54.0mu\sum_{\stackrel{{\scriptstyle\vec{m}_{N+1},\vec{m}^{\prime}_{N+1}}}{{\mathcal{W}_{N+1}=\mathcal{W}^{\vec{m}_{N+1},\vec{m}^{\prime}_{N+1}}}}}\mkern-54.0muD^{\vec{m}_{N+1},\vec{m}^{\prime}_{N+1}}(\tilde{\rho})$
(83) $\displaystyle=$
$\displaystyle\mathcal{D}^{\mathrm{RC}}_{\mathcal{W}_{N+1}}(\rho).$

Therefore, the operator $\mathcal{D}^{\mathrm{RC}}_{\mathcal{W}_{N+1}}(\rho)$
is obtained recursively from
$\mathcal{D}^{\mathrm{RC}}_{\mathcal{W}_{N}}(\rho)$ in terms of a sum over
$\mathbf{m}$, $\mathbf{m}^{\prime}$ and $\mathcal{W}_{N}$. The required number
of computations does not scale with $N$ (for $\mathbf{m}^{(\prime)}$ we have
$2^{4}$ terms for $d=2$) and the sum over $\mathcal{W}_{N}$ scales as $N^{2}$.
Thus, the computations performed, due to the grouping with fixed work values,
reduces from an exponential scaling of $2^{8N}$ ($d=2$) to a quadratic
$N^{2}$.

##  Appendix D Equivalence between one- and two-pointer work PDFs for a
perfectly thermalizing adiabatic engine

In this section, we show that for a perfectly thermalizing adiabatic engine
diagnosed via repeated contacts the one- and two-pointer work PDFs are equal
independent of the dimensionality $d$ of the Hilbert space of the working
substance. In order to achieve our objective, we begin with the one- and two-
pointer work PDFs as defined in Eqs. (25) and (26), i.e.,

$\displaystyle P^{\mathrm{RC,1P}}(W)$ $\displaystyle=$
$\displaystyle\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{W}^{\vec{m}}-\mathcal{W}^{\vec{m}\prime})^{2}}g_{\Sigma^{2}}(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}}),$
$\displaystyle P^{\mathrm{RC,2P}}(W)$ $\displaystyle=$
$\displaystyle\sum_{\vec{m},\vec{m}^{\prime}}D^{\vec{m},\vec{m}^{\prime}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{W}^{\vec{m}}-\mathcal{W}^{\vec{m}\prime})^{2}}e^{-\frac{1}{8\Sigma^{2}}(\mathcal{Q}^{\vec{m}}-\mathcal{Q}^{\vec{m}\prime})^{2}}g_{\Sigma^{2}}(W-\mathcal{W}^{\vec{m},\vec{m}^{\prime}}).$

Comparing the above equations, it is easy to see that the PDFs are equal when
either $\mathcal{Q}^{\vec{m}}=\mathcal{Q}^{\vec{m}\prime}$ or when
$D^{\vec{m},\vec{m}^{\prime}}=0$.

We restrict ourselves to a single cycle ($N=1$) perfectly thermalizing engine,
such that the elements of
$\vec{m}^{(\prime)}=(m_{1}^{(\prime)},m_{2}^{(\prime)},m_{3}^{(\prime)},m_{4}^{(\prime)})$.
To proceed, we consider all possible relations between the elements of
$\vec{m}$ and $\vec{m}^{\prime}$, i.e.,

  1. 1.

$m_{2}=m_{2}^{\prime}$ and $m_{3}=m_{3}^{\prime}$ (independent of
$m_{1}^{(\prime)}$ and $m_{4}^{(\prime)}$),

  2. 2.

$m_{2}\neq m_{2}^{\prime}$ and $m_{3}=m_{3}^{\prime}$ (independent of
$m_{1}^{(\prime)}$ and $m_{4}^{(\prime)}$), and

  3. 3.

$m_{3}\neq m_{3}^{\prime}$ (independent of $m_{1}^{(\prime)}$,
$m_{2}^{(\prime)}$, and $m_{4}^{(\prime)}$),

and show that they can be classified into one of the two categories
($\mathcal{Q}^{\vec{m}}=\mathcal{Q}^{\vec{m}\prime}$ or
$D^{\vec{m},\vec{m}^{\prime}}=0$).

Case (i) (_Equality due to matching heat outcomes_): For a single cycle
considered herein the heat outcomes
$\mathcal{Q}^{\vec{m}}=e^{m_{3}}_{3}-e^{m_{2}}_{2}$ and
$\mathcal{Q}^{\vec{m}^{\prime}}=e^{m_{3}^{\prime}}_{3}-e^{m_{2}^{\prime}}_{2}$,
using Eq. (23). Since $m_{2}=m_{2}^{\prime}$ and $m_{3}=m_{3}^{\prime}$,
clearly $\mathcal{Q}^{\vec{m}}=\mathcal{Q}^{\vec{m}\prime}$.

Case (ii) (_Equality due to perfect thermalization_): For a single cycle
engine, the explicit expression for the coefficient
$D^{\vec{m},\vec{m}^{\prime}}$ reads,

$\displaystyle
D^{\vec{m},\vec{m}^{\prime}}=\mathrm{Tr}[\mathcal{P}^{m_{4}}_{4}\tilde{U}\mathcal{P}^{m_{3}}_{3}\,\Phi_{h}\left(\mathcal{P}^{m_{2}}_{2}U\mathcal{P}^{m_{1}}_{1}\rho\mathcal{P}^{m_{1}^{\prime}}_{1}U^{\dagger}\mathcal{P}^{m_{2}^{\prime}}_{2}\right)\,\mathcal{P}^{m_{3}^{\prime}}_{3}\tilde{U}^{\dagger}\mathcal{P}^{m_{4}^{\prime}}_{4}],$
(84)

with the operators $\mathcal{P}^{m_{j}}_{j}$ projecting into the eigenstates
of Hamiltonian $H_{j}$, $\rho$ being the density matrix of the working
substance, and $U$ ($\tilde{U}$) being the unitary forward (reversed) time-
evolution operators. Since the engine is perfectly thermalizing,
$\Phi_{h}(\rho)=\rho_{h}\mathrm{Tr}[\rho]$ using Eq. (34), we obtain

$\displaystyle\Phi_{h}\left(\mathcal{P}^{m_{2}}_{2}U\mathcal{P}^{m_{1}}_{1}\rho\mathcal{P}^{m_{1}^{\prime}}_{1}U^{\dagger}\mathcal{P}^{m_{2}^{\prime}}_{2}\right)=\rho_{h}\mathrm{Tr}[\mathcal{P}^{m_{2}}_{2}U\mathcal{P}^{m_{1}}_{1}\rho\mathcal{P}^{m_{1}^{\prime}}_{1}U^{\dagger}\mathcal{P}^{m_{2}^{\prime}}_{2}]=0.$
(85)

The rightmost equality is because we consider $m_{2}\neq m_{2}^{\prime}$ in
this case giving us the trace in the middle term to be zero. Therefore, for
case (ii) $D^{\vec{m},\vec{m}^{\prime}}=0$.

Case (iii) (_Equality due to adiabaticity_): We further sub-divide case (iii)
into two sub-categories: (a) $m_{3}\neq m_{4}$ or $m_{3}^{\prime}\neq
m_{4}^{\prime}$ and (b) $m_{3}=m_{4}$ and $m_{3}^{\prime}=m_{4}^{\prime}$.

For case (iiia), since the evolution during the work strokes is adiabatic,
i.e., if the system starts in an eigenstate of the initial Hamiltonian, it
will end in the corresponding eigenstate of the final Hamiltonian, either

$\displaystyle\mathcal{P}^{m_{4}}_{4}\tilde{U}\mathcal{P}^{m_{3}}_{3}$
$\displaystyle=$ $\displaystyle 0\quad\mathrm{when~{}}m_{3}\neq
m_{4}~{}~{}\mathrm{or}$
$\displaystyle\mathcal{P}^{m_{3}^{\prime}}_{3}\tilde{U}^{\dagger}\mathcal{P}^{m_{4}^{\prime}}_{4}$
$\displaystyle=$ $\displaystyle 0\quad\mathrm{when~{}}m_{3}^{\prime}\neq
m_{4}^{\prime}.$ (86)

In either case the coefficient $D^{\vec{m},\vec{m}^{\prime}}=0$ [see Eq.
(84)].

For the case (iiib), since $m_{3}\neq m_{3}^{\prime}$, $m_{3}=m_{4}$, and
$m_{3}^{\prime}=m_{4}^{\prime}$ implies $m_{4}\neq m_{4}^{\prime}$. Therefore,
again $D^{\vec{m},\vec{m}^{\prime}}=0$ due to the overall trace in Eq. (84).
Thus, overall for case (iii) $D^{\vec{m},\vec{m}^{\prime}}=0$.

Overall, with the three cases above, all possible combinations of $\vec{m}$
and $\vec{m}^{\prime}$ are excluded that would lead to a difference between
the one- and the two-pointer work PDFs. Moreover, our proof does not depend on
the working substance of the engine and hence holds for all thermalizing
adiabatic quantum Otto engines diagnosed via repeated contacts.

## References

  * Müller [2007] I. Müller,  _A History of Thermodynamics_ (Springer Verlag, Berlin, 2007).
  * Blickle and Bechinger [2000] V. Blickle and C. Bechinger, Realization o micrometre-sized stochastic heat engine, Nature Physics 33, 4427 (2000).
  * Scovil and Schulz-DuBois [1959] H. E. D. Scovil and E. O. Schulz-DuBois, Three-level masers as heat engines, Phys. Rev. Lett. 2, 262 (1959).
  * Alicki [1979] R. Alicki, The quantum open system as a model of the heat engine, J. Phys. A 12, L103 (1979).
  * Kosloff [1984] R. Kosloff, A quantum mechanical open system as a model of a heat engine, J. Chem. Phys. 80, 1625 (1984).
  * Geva and Kosloff [1992] E. Geva and R. Kosloff, A quantum mechanical heat engine operating in finite time, J. Chem. Phys. 97, 4398 (1992).
  * Bender  _et al._ [2000] C. M. Bender, D. C. Brody, and B. K. Meister, Quantum mechanical carnot engine, J. Phys. A 33, 4436 (2000).
  * Feldmann and Kosloff [2003] T. Feldmann and R. Kosloff, Quantum four-stroke heat engine: Thermodynamic observables in a model with intrinsic friction, Phys. Rev. E 68, 016101 (2003).
  * Scully  _et al._ [2003] M. O. Scully, M. S. Zubairy, G. S. Agarwal, and H. Walther, Extracting work from a single heat bath via vanishing quantum coherence, Science 299, 862 (2003).
  * Uzdin [2016] R. Uzdin, Coherence-induced reversibility and collective operation of quantum heat machines via coherence recycling, Phys. Rev. App. 6, 024004 (2016).
  * Francica  _et al._ [2019] G. Francica, J. Goold, and F. Plastina, Role of coherence in the nonequilibrium thermodynamics of quantum systems, Phys. Rev. E 99, 042105 (2019).
  * Dann and Kosloff [2020] R. Dann and R. Kosloff, Quantum signatures in the quantum carnot cycle, New J. Phys. 22, 013055 (2020).
  * Latune  _et al._ [2021] C. L. Latune, I. Sinayskiy, and P. F, Roles of quantum coherences in thermal machines (2021).
  * Quan  _et al._ [2005] H. T. Quan, P. Zhang, and C. P. Sun, Quantum heat engine with multilevel quantum systems, Phys. Rev. E 72, 056110 (2005).
  * Gelbwaser-Klimovsky  _et al._ [2018] D. Gelbwaser-Klimovsky, A. Bylinskii, D. Gangloff, R. Islam, A. Aspuru-Guzik, and V. Vuletic, Single-atom heat machines enabled by energy quantization, Phys. Rev. Lett. 120, 170601 (2018).
  * Halpern  _et al._ [2019] N. Y. Halpern, C. D. White, S. Gopalakrishnan, and G. Refael, Quantum engine based on many-body localization, Phys. Rev. B 99, 024203 (2019).
  * Gelbwaser-Klimovsky  _et al._ [2015] D. Gelbwaser-Klimovsky, W. Niedenzu, P. Brumer, and G. Kurizki, Power enhancement of heat engines via correlated thermalization in a three-level “working fluid”, Sci. Rep. 5, 1 (2015).
  * Wang  _et al._ [2009] H. Wang, S. Liu, and J. He, Thermal entanglement in two-atom cavity qed and the entangled quantum otto engine, Phys. Rev. E 79, 041113 (2009).
  * Watanabe  _et al._ [2020] G. Watanabe, B. P. Venkatesh, P. Talkner, M.-J. Hwang, and A. del Campo, Quantum statistical enhancement of the collective performance of multiple bosonic engines, Phys. Rev. Lett. 124, 210603 (2020).
  * Yi  _et al._ [2017] J. Yi, P. Talkner, and Y. W. Kim, Single-temperature quantum engine without feedback control, Phys. Rev. E 96, 022108 (2017).
  * Ding  _et al._ [2018] X. Ding, J. Yi, Y. W. Kim, and P. Talkner, Measurement-driven single temperature engine, Phys. Rev. E 98, 042122 (2018).
  * Campisi  _et al._ [2015] M. Campisi, J. Pekola, and R. Fazio, Nonequilibrium fluctuations in quantum heat engines: theory, example and possible solid state experiments, New J. Phys. 17, 035012 (2015).
  * Roßnagel  _et al._ [2016] J. Roßnagel, S. T. Dawkins, K. N. Tolazzi, O. Abah, E. Lutz, F. Schmidt-Kaler, and K. Singer, A single-atom heat engine, Science 352, 325 (2016).
  * Maslennikov  _et al._ [2019] G. Maslennikov, S. Ding, R. Hablützel, J. Gan, A. Roulet, S. Nimmrichter, J. Dai, V. Scarani, and D. Matsukevich, Quantum absorption refrigerator with trapped ions, Nat. Comm. 10, 1 (2019).
  * Holmes  _et al._ [2020] Z. Holmes, J. Anders, and F. Mintert, Enhanced energy transfer to an optomechanical piston from indistinguishable photons, Phys. Rev. Lett. 124, 210601 (2020).
  * Klatzow  _et al._ [2019] J. Klatzow, J. N. Becker, P. M. Ledingham, C. Weinzetl, K. T. Kaczmarek, D. J. Saunders, J. Nunn, I. A. Walmsley, R. Uzdin, and E. Poem, Experimental demonstration of quantum effects in the operation of microscopic heat engines, Phys. Rev. Lett. 122, 110601 (2019).
  * Peterson  _et al._ [2019] J. P. Peterson, T. B. Batalhao, M. Herrera, A. M. Souza, R. S. Sarthour, I. S. Oliveira, and R. M. Serra, Experimental characterization of a spin quantum heat engine, Phys. Rev. Lett. 123, 240601 (2019).
  * Quan  _et al._ [2006] H. T. Quan, Y. D. Wang, Y. X. Liu, C. P. Sun, and F. Nori, Maxwell’s demon assisted thermodynamic cycle in superconducting quantum circuits, Phys. Rev. Lett. 79, 180402 (2006).
  * Zheng  _et al._ [2016] Y. J. Zheng, P. Hänggi, and D. Poletti, Occurence of discontinuities in the performence of finite-time quantum Otto cycles, Phys. Rev E 94, 012137 (2016).
  * Denzler  _et al._ [2021] T. Denzler, F. F. G. Santos, E. Lutz, and R. M. Serra, Nonequilibriuum fluctuations of a quantum heat engine (2021), arXiv:2104.13427 .
  * Kurchan [2000] J. Kurchan, A quantum fluctuation theorem (2000), arXiv:cond-mat/0007360 .
  * Tasaki [2000] H. Tasaki, Jarzynski relations for quantum systems and some applications (2000), arXiv:cond-mat/0009244 .
  * Talkner  _et al._ [2007] P. Talkner, E. Lutz, and P. Hänggi, Fluctuation theorems: Work is not an observable, Phys. Rev. E 75, 050102 (2007).
  * Kammerlander and Anders [2016] P. Kammerlander and J. Anders, Coherence and measurement in quantum thermodynamics, Sci. Rep. 6, 1 (2016).
  * Vinjanampathy and Anders [2016] S. Vinjanampathy and J. Anders, Quantum thermodynamics, Cont. Phys. 57, 545 (2016).
  * Nieuwenhuizen and Allahverdyan [2002] T. M. Nieuwenhuizen and A. E. Allahverdyan, Statistical thermodynamics of quantum Brownian motion: Construction of perpetuum mobile of the second kind, Phys. Rev. E 66, 036102 (2002).
  * Allahvedyan [2014] A. E. Allahvedyan, Nonequilibrium quantum fluctuations of work, Phys. Rev. A 90, 032137 (2014).
  * Talkner and Hänggi [2016] P. Talkner and P. Hänggi, Aspects of quantum work, Phys. Rev. E 93, 022131 (2016).
  * Perarnau-Llobet  _et al._ [2017] M. Perarnau-Llobet, E. Bäumer, K. V. Hovhannisyan, M. Huber, and A. Acin, No-go theorem for the characterization of work fluctuations in coherent quantum systems, Phys. Rev. Lett. 118, 070601 (2017).
  * Talkner and Hänggi [2020] P. Talkner and P. Hänggi, Colloquium: Statistical mechanics and thermodynamics at strong coupling: Quantum and classical, Rev. Mod. Phys. 92, 041002 (2020).
  * Wiseman and Milburn [2009] H. M. Wiseman and G. J. Milburn,  _Quantum measurement and control_ (Cambridge university press, 2009).
  * Thingna  _et al._ [2012] J. Thingna, J.-S. Wang, and P. Hänggi, Generalized gibbs state with modified redfield solution: Exact agreement up to second order, J. Chem. Phys. 136, 194110 (2012).
  * Rezek and Kosloff [2006] Y. Rezek and R. Kosloff, Irreversible performance of a quantum harmonic heat engine, New J. Phys. 8, 83 (2006).
  * Plastina  _et al._ [2014] F. Plastina, A. Alecce, T. J. Apollaro, G. Falcone, G. Francica, F. Galve, N. L. Gullo, and R. Zambrini, Irreversible work and inner friction in quantum thermodynamic processes, Phys. Rev. Lett. 113, 260601 (2014).
  * Thingna  _et al._ [2019] J. Thingna, M. Esposito, and F. Barra, Landau-zener lindblad equation and work extraction from coherences, Phys. Rev. E 99, 042142 (2019).
  * von Neumann [1932] J. von Neumann,  _Mathematische Grundlagen der Quantenmechanik_ (Springer Verlag, Berlin, 1932) English translation: 1955 Mathematical Foundations of Quantum Mechanics, (Princeton University Press, Princeton) .
  * Thingna and Talkner [2020] J. Thingna and P. Talkner, Quantum measurements of sums, Phys. Rev. A 102, 012213 (2020).
  * Ito  _et al._ [2019] K. Ito, C. Jiang, and G. Watanabe, Universal bounds for fluctations in small heat engines (2019), arXiv:1910.08096 .
  * Verley  _et al._ [2014] G. Verley, M. Esposito, and C. Van den Broeck, The unlikely carnot efficiency, Nature Comm. 5, 4721 (2014).
  * Shevchenko  _et al._ [2010] S. N. Shevchenko, S. Ashhab, and F. Nori, Landau–zener–stückelberg interferometry, Phys. Rep. 492, 1 (2010).
  * Barnes and Das Sarma [2012] E. Barnes and S. Das Sarma, Analytically solvable driven time-dependent two-level quantum systems, Phys. Rev. Lett. 109, 060401 (2012).
  * Thingna  _et al._ [2014] J. Thingna, P. Hänggi, R. Fazio, and M. Campisi, Geometric quantum pumping in the presence of dissipation, Phys. Rev. B 90, 094517 (2014).
  * Bandyopadhyay and Thingna [2021] J. N. Bandyopadhyay and J. Thingna, Lie algebraic floquet engineering protocol for periodically driven quantum systems (2021), arXiv:2103.15923 .
  * Kosloff and Rezek [2017] R. Kosloff and Y. Rezek, The quantum harmonic otto cycle, Entropy 19, 136 (2017).
  * Watanabe  _et al._ [2017] G. Watanabe, B. P. Venkatesh, P. Talkner, and A. del Campo, Quantum performance of thermal machines over many cycles, Phys. Rev. Lett. 118, 050601 (2017).