# Many-body quantum heat engines based on free-fermion systems

Vincenzo Roberto Arezzo  SISSA, Via Bonomea 265, I-34135 Trieste, Italy
Davide Rossini  Dipartimento di Fisica dell’Università di Pisa and INFN, Largo
Pontecorvo 3, I-56127 Pisa, Italy  Giulia Piccitto  Dipartimento di Matematica
e Informatica, Università di Catania, Viale Andrea Doria 6, I-95128, Catania,
Italy

###### Abstract

We study the performances of an imperfect quantum many-body Otto engine based
on free-fermion systems. Starting from the thermodynamic definitions of heat
and work along ideal isothermal, adiabatic, and isochoric transformations, we
generalize these expressions in the case when the hypotheses of ideality are
relaxed (i.e., nonperfect thermalization with the external baths, as well as
nonperfect quantum adiabaticity in the unitary dynamic protocols). These
results are used to evaluate the work and the power delivered by an imperfect
quantum many-body heat engine in a finite time, whose working substance is
constituted by a quantum Ising chain in a transverse field: We discuss the
emerging optimal working points as functions of the various model parameters.

##  I Introduction

Thermodynamics is one of the most fascinating topics in physics. Since the
Industrial Revolution, the necessity of understanding mechanisms underlying
the conversion of the heat in useful work has focused efforts of the
researchers in this direction, leading to the statement of the three laws of
thermodynamics. Much more recently, with the advent of quantum technologies
and the miniaturization of the devices exchanging heat and work down to the
nanoscale, it has become relevant to understand these issues in the novel
scenario of quantum mechanics. Unfortunately, while the theory of classical
thermodynamics is well posed, its generalization to the quantum realm raises
some conceptual puzzles and is still an open research field. Seminal works
along this direction, proposing a generalization of the concepts of heat and
work for quantum systems [1, 2], date back to the eighties.

Only recently, the renewed interest in the possibility of quantum heat engines
has led to a prolific scientific production, starting from the pioneer
proposal of the maser as a first example of quantum engine [3]. A quantum heat
engine is simply defined as an engine constituted by a quantum systems that
can produce work undergoing a suitable thermodynamic cycle. A series of papers
[4, 5, 6, 7] have proposed possible realizations of few-body heat engines,
based both on the Carnot cycle [8, 9] and on the Otto cycle [10, 11, 12, 13].
Some experimental realizations have been obtained by means of different
platforms [14, 15, 16, 17, 18, 19], as also discussed in Ref. [20].

A promising direction to investigate quantum advantage is that of engines
whose working substance is constituted by quantum many-body systems [21, 22].
Although this research line is hindered by technical difficulties, from both
analytical and numerical point of views, some preliminary results have been
proposed, both with gases of interacting atoms [23, 24, 25, 26, 27, 28, 29]
and with interacting quantum spins [30, 31, 32, 33, 34, 35, 36]. It has also
been shown how criticality may have an impact on the enhancement of the engine
performance [37, 26, 34]. Unless some exceptions [38], all results in the
literature have been obtained for ideal engines undergoing perfect
thermodynamic transformations.

In this paper we study a many-body quantum Otto engine with a free-fermion
medium, in which the assumption of perfect thermodynamic transformations can
be relaxed. We first exploit the exact integrability of the model through the
Bogoliubov de-Gennes formalism to derive some general analytic expressions for
the heat exchanged and the work produced by a free-fermion system of arbitrary
size along some paradigmatic thermodynamic processes: the nonperfect isochoric
(static Hamiltonian system, coupled to an external bath) and the nonperfect
adiabatic (isolated quantum system) transformation separately, as well as the
nonperfect adiabatic transformation followed by a nonperfect isochoric one. To
model nonunitary processes, we describe the microscopics of the system-bath
dynamics through a thermal bath of harmonic oscillators at a given
temperature, quadratically coupled to the system, by means of a nonlocal
Lindblad master equation ensuring thermalization at long times.

In the second part of the paper, we apply the above analytic results to the
quantum Otto cycle, with the aim of finding the best parameters which optimize
the engine performances. In particular, we numerically analyze the work output
and the power delivered by the Ising quantum Otto engine when quantum
adiabaticity is violated, thus generalizing the results obtained by two of us
in Ref. [34] for the ideal scenario and at finite size. We observe the
emergence of a nonmonotonic behavior in the power with the various parameters
of our model, thus signalling the existence of an optimal working point.

The paper is structured as follows. In Sec. II, we recap the basic concepts of
quantum thermodynamics and clarify the differences between classical and
quantum realms. In Sec. III, we derive the expressions for the heat and the
work along ideal isotherm, adiabatic, and isochoric transformations, for a
generic free-fermion system that can be coupled to a thermal bath. We then
introduce the quantum Otto cycle (Sec. IV), whose working substance is a free-
fermion system, and analyze the performance of such engine operating under the
assumption of nonperfect thermalization and/or quantum adiabaticity. In Sec.
V, we present our numerical simulations for the case of an Ising Otto engine
when nonadiabatic transformations are taken into account. Finally, in Sec. VI
we draw our conclusions. The appendixes provide further details on the system-
bath coupling scheme we adopt within the Lindblad formalism (App. A), the
ideal Ising Otto engine in the thermodynamic limit (App. B), and the way in
which the Hamiltonian control parameter should be varied to implement a proper
Otto cycle producing useful work (App. C).

##  II Thermodynamics of quantum systems

One of the basic principles governing the exchange of heat and work between
two or more systems is the first law of thermodynamics,

$dE=\delta Q-\delta W,$ (1)

which formally assesses the energy conservation principle: the variation of
internal energy $E$ of a given system can be interpreted as the difference
between the heat exchanged with the external environment and the work
performed by the system itself. In the following, we adopt the convention that
$\delta Q>0$ (positive) refers to the heat absorbed by the system and $\delta
W>0$ (positive) refers to work performed by the system. Notice that, despite
the energy is an observable and thus its differential $dE$ is exact, those of
the heat $Q$ and of the work $W$ are not. For this reason, we indicate their
infinitesimal variation with $\delta Q$ and $\delta W$. While for classical
systems, given a specific thermodynamic transformation, the heat and the work
are well defined quantities, for quantum systems this is not always the case
and additional care should be taken.

To be more quantitative, we start from the Hamiltonian $H=H(\lambda)$ of a
generic quantum system, which is supposed to depend on some macroscopic
parameter $\lambda$ that can be suitably controlled (e.g., the magnetic field
strength or the volume itself). Its expectation value over the actual
(generally mixed) state $\rho(\lambda)$ of the system defines the internal
energy

$E\equiv\braket{H}_{\rho(\lambda)}=\text{Tr}\big{[}H(\lambda)\,\rho(\lambda)\big{]},$
(2)

where ${\rm Tr}[\,\cdot\,]$ denotes the trace operation, so that

$dE=\text{Tr}\big{[}dH(\lambda)\,\rho(\lambda)\big{]}+\text{Tr}\big{[}H(\lambda)\,d\rho(\lambda)\big{]}.$
(3)

The first term in the r.h.s. accounts for a modification of the spectral
structure of the system as a response to a variation of $\lambda$, while the
second term describes a variation of the system state. One could naively
identify the first term as the work done by the system itself, quantifying the
mechanical energy exchanged with the external reservoir, and the second one as
the heat absorbed by the system from the environment [1, 2, 4]. However, this
interpretation can be misleading, as, for example, when considering nonthermal
reservoirs. In what follows, we consider the heat as the variation of energy
due to the interaction with a thermal reservoir and the work as the variation
of energy due to a change in the Hamiltonian parameters.

Once the heat exchange has been introduced, one can write the Clausius
inequality,

$\sum_{i}\beta_{i}Q_{i}\leq 0,$ (4)

where $Q_{i}$ denote the various heats exchanged with a set of surrounding
reservoirs at temperatures $\beta_{i}^{-1}$ (hereafter we will always work in
units of $\hbar=k_{B}=1$). Equation (4) can be combined with Eq. (1) to show
that, for a single working medium operating between two reservoirs, cold (c)
and hot (h), at different temperatures and varying its energy, only four
operation modes are allowed [10]: (i) Heat Engine (production of work by heat
absorption from the hot reservoir, $Q_{c}<0,\,Q_{h}>0,\,W>0$); (ii)
Refrigerator (heat transfer from the cold to the hot bath by energy absorption
of the medium, $Q_{c}>0,\,Q_{h}<0,\,W<0$); (iii) Thermal Accelerator (heat
transfer from the hot to the cold bath by energy absorption,
$Q_{c}<0,\,Q_{h}>0,\,W<0$); (iv) Heater (heating up the two baths by energy
absorption, $Q_{c}<0,\,Q_{h}<0,\,W<0$).

##  III Thermodynamics of free-fermion systems

We now derive expressions for the work performed and the heat exchanged by
free-fermion systems during some paradigmatic and ideal thermodynamic
transformations. Namely, we focus on perfect isotherm, adiabatic (unitary),
and isochoric transformations.

###  III.1 Free-fermion systems

We consider a lattice model of spinless fermions coupled through a generic
quadratic Hamiltonian

$H=\sum_{i,j}D_{i,j}c_{i}^{\dagger}c_{j}+\tfrac{1}{2}\big{(}O_{i,j}c_{i}^{\dagger}c_{j}^{\dagger}+\text{h.c.}\big{)},$
(5)

where $c_{i}^{(\dagger)}$ are anticommuting fermionic annihilation (creation)
operators on the $i$-th site ($i=1,\ldots,N$), while $D$ and $O$ are complex
$N\times N$ matrices satisfying $D=D^{\dagger}$ and $O=-O^{T}$ (so to respect
hermiticity of $H$). By introducing the Nambu spinor
$\Psi=\big{(}c_{1},\dots,c_{N},c_{1}^{\dagger},\dots,c_{N}^{\dagger}\big{)}^{T}$,
Eq. (5) can be written as

$H=\Psi^{\dagger}\mathbb{H}\Psi+\text{Tr}[D],$ (6)

where the $2N\times 2N$ matrix

$\mathbb{H}=\frac{1}{2}\begin{pmatrix}\phantom{-}D&\phantom{-}O\\\
-O^{*}&-D^{*}\end{pmatrix}$ (7)

is the so called Bogoliubov de-Gennes matrix, which can be diagonalized by a
unitary transformation $\mathbb{U}$, such that

$\mathbb{H^{D}}=\mathbb{U}^{\dagger}\mathbb{H}\mathbb{U}=\text{diag}\big{\\{}\epsilon_{1}\ldots\epsilon_{N},-\epsilon_{1}\ldots-\epsilon_{N}\big{\\}}.$
(8)

Here $\epsilon_{k}>0$ denote the positive eigenvalues of $\mathbb{H}$,
$(k=1,\dots,N)$, while $\mathbb{U}$ is a Bogoliubov rotation that can be cast
in the block form

$\mathbb{U}=\begin{pmatrix}\mathbb{u}&\mathbb{v}\\\
\phantom{*}\mathbb{v}^{*}&\phantom{*}\mathbb{u}^{*}\end{pmatrix}.$ (9)

This transformation is constructed in such a way to preserve the particle-
antiparticle symmetry of the model and defines a new set of anticommuting
fermions $b_{k}^{(\dagger)}$ (the Bogoliubov quasiparticles) through the
relation $\Psi=\mathbb{U}\Phi$, with
$\Phi=\big{(}b_{1},\dots,b_{N},b^{\dagger}_{1},\dots,b_{N}^{\dagger}\big{)}^{T}$.
It is always possible to gauge away the term $\text{Tr}[D]$ in Eq. (6) by
shifting the Hamiltonian, thus hereafter we simply neglect it.

Summarizing, in the diagonal basis, the Hamiltonian (5) reads

$H=\sum_{k>0}\omega_{k}(\lambda)\left(b^{\dagger}_{k}b_{k}-\tfrac{1}{2}\right)$
(10)

and corresponds to a free-fermion model, being $\omega_{k}=2\epsilon_{k}$ the
dispersion relation of the (noninteracting) Bogoliubov $b_{k}$-quasiparticles
[39]. The spectrum of this Hamiltonian, which for our purposes we assume to be
non degenerate, should depend on the control parameter $\lambda$, therefore we
have $\omega_{k}\equiv\omega_{k}(\lambda)$.

###  III.2 Ideal isotherm, adiabatic and isochoric transformations

For a quadratic fermionic system (5) prepared in the thermal state

$\rho=Ze^{-\beta H},$ (11)

we can write the internal energy as

$E=\sum_{k}\omega_{k}(\lambda)\left(Z^{-1}\text{Tr}\big{[}b^{\dagger}_{k}b_{k}e^{-\beta
H}\big{]}-\tfrac{1}{2}\right),$ (12)

where $\beta^{-1}$ denotes the temperature and $Z=\text{Tr}\big{[}e^{-\beta
H}\big{]}$ the partition function. For a thermal density matrix as in Eq.
(11), the populations of the $b_{k}$ quasiparticles follow the Fermi-Dirac
distribution function:

$\text{Tr}\big{[}b^{\dagger}_{k}b_{k}\,\rho\big{]}=\big{[}1+e^{-\beta\omega_{k}(\lambda)}\big{]}^{-1}\equiv
f\big{[}\beta,\omega_{k}(\lambda)\big{]},$ (13)

so that Eq. (12) takes the simple form:
$E=\sum_{k}\omega_{k}(\lambda)\,\big{\\{}f\big{[}\beta,\omega_{k}(\lambda)\big{]}-1/2\big{\\}}$.
The partition function can be easily computed to give:

$Z=e^{\frac{\beta}{2}\sum_{k}\omega_{k}(\lambda)}\prod_{k}\big{[}1+e^{-\beta\omega_{k}(\lambda)}\big{]}.$
(14)

The heat exchanged and the work performed by the system can be thus evaluated
straightforwardly in some paradigmatic transformations according to the
following.

####  III.2.1 Ideal isothermal transformation

An ideal isotherm can be obtained by considering a slow variation of the
control parameter $\lambda$, such that the system is assumed to be always in
thermal equilibrium with an environment at a given fixed temperature
$\beta^{-1}$. In that case, the Helmholtz free-energy is given by

$\mathcal{F}=-\frac{1}{\beta}\log
Z=-\frac{\kappa}{2}-\frac{1}{\beta}\sum_{k}\log\big{[}1+e^{-\beta\omega_{k}(\lambda)}\big{]},$
(15)

where $\kappa=\sum_{k}\omega_{k}(\lambda)$.

The work performed by the system along an isothermal variation of the external
parameter $\lambda$, from $\lambda_{i}$ to $\lambda_{f}$, can be evaluated as
the opposite of the variation of the free-energy, that can be expressed in the
integral form

$\displaystyle W=-\Delta\mathcal{F}$
$\displaystyle=-\int_{\lambda_{i}}^{\lambda_{f}}\frac{\partial\mathcal{F}}{\partial\lambda}d\lambda$
(16)
$\displaystyle=-\int_{\lambda_{i}}^{\lambda_{f}}\bigg{\\{}\sum_{k}\omega^{\prime}_{k}(\lambda)\Big{[}f\big{[}\beta,\omega_{k}(\lambda)\big{]}-\tfrac{1}{2}\Big{]}\bigg{\\}}d\lambda,\;$

where the prime denotes the derivative with respect to $\lambda$. Note that,
depending on the derivative of spectrum $\omega_{k}^{\prime}(\lambda)$, the
integrand in Eq. (16) may behave monotonically with the temperature
$\beta^{-1}$: If $\omega_{k}^{\prime}(\lambda)>0$ $\forall k$, then the free-
energy is an increasing function of $\beta$; if
$\omega_{k}^{\prime}(\lambda)<0$ $\forall k$, then the trend is reversed
(details in App. C).

The exchanged heat can be obtained through the first law of thermodynamics
(1), $Q=\Delta E+W$, and using Eq. (12).

####  III.2.2 Ideal adiabatic transformation

By definition, an adiabatic transformation is such that it occurs without heat
exchange ($Q=0$). We thus refer to this wording at any time the quantum system
is isolated (not coupled to an external environment). This has not to be
confused with the quantum adiabaticity condition, i.e., when the unitary
evolution satisfies the requirements of the quantum adiabatic theorem [40]
(thus a transformation can be adiabatic but nonideal, in the sense that it is
not quantum adiabatic — see Sec. IV.3).

Since, in general, an adiabatic transformation is not quasistatic, one cannot
evaluate the work through the Helmholtz free-energy. However, in the ideal
case, the transformation is induced by a variation of $\lambda$, from
$\lambda_{i}$ to $\lambda_{f}$, which is so slow that one can invoke the
quantum adiabatic theorem [41] (the spectrum is nondegenerate) and calculate
the work done by the system by exploiting the first law of thermodynamics:

$W=-\Delta
E=\big{\langle}H(\lambda_{i})\big{\rangle}_{\rho(\lambda_{i})}-\big{\langle}H(\lambda_{f})\big{\rangle}_{\tilde{\rho}(\lambda_{i})}.$
(17)

Here $\rho(\lambda_{i})$ denotes the equilibrium thermal state for a system
described by the Hamiltonian $H(\lambda_{i})$ at temperature $\beta^{-1}$, at
the beginning of the adiabatic transformation, while
$\tilde{\rho}(\lambda_{i})$ is the quantum adiabatically evolved state. Such
latter state is constructed using the projectors on the eigenstates of
$H(\lambda_{f})$, with the same energy-level populations of
$\rho(\lambda_{i})$, therefore:

$\displaystyle W$
$\displaystyle=\sum_{k}\big{[}\omega_{k}(\lambda_{i})-\omega_{k}(\lambda_{f})\big{]}\Big{\\{}f\big{[}\beta,\omega_{k}(\lambda_{i})\big{]}-\tfrac{1}{2}\Big{\\}},$
$\displaystyle=-\int_{\lambda_{i}}^{\lambda_{f}}\bigg{\\{}\sum_{k}\omega^{\prime}_{k}(\lambda)\Big{[}f\big{[}\beta,\omega_{k}(\lambda_{i})\big{]}-\tfrac{1}{2}\Big{]}\bigg{\\}}d\lambda.$
(18)

Note that in general, as in classical thermodynamics, the work done by the
system in an adiabatic transformation [cf. Eq. (18)] is not larger than the
one done in an isothermal transformation [cf. Eq. (16)]. This follows from the
fact that

$\int_{\lambda_{i}}^{\lambda_{f}}{\frac{\omega_{k}^{\prime}(\lambda)}{1+e^{\beta\omega_{k}(\lambda)}}}\,d\lambda\leq\int_{\lambda_{i}}^{\lambda_{f}}{\frac{\omega_{k}^{\prime}(\lambda)}{1+e^{\beta\omega_{k}(\lambda_{i})}}}\,d\lambda,$
(19)

which can be proved by solving the integral and using the fact that
$-\log\left(1+e^{-x}\right)$ is a concave function.

####  III.2.3 Ideal isochoric transformation

An isochoric transformation is such that the system is put in contact with a
thermal reservoir, without variation in the Hamiltonian parameter $\lambda$.
No work is thus performed during this transformation ($W=0$). We model this
situation in such a way that the system is prepared in a thermal state
$\rho_{i}$ at temperature $\beta^{-1}_{1}$ and is put in contact with a bath
at temperature $\beta_{2}^{-1}$, keeping the Hamiltonian $H$ fixed. After a
time $\tau>t^{\text{th}}$, being $t^{\text{th}}$ the thermalization time, the
system will be in the thermal equilibrium state
$\rho_{f}=e^{-\beta_{2}H}/\text{Tr}[e^{-\beta_{2}H}]$. The heat exchanged
during this process can be thus evaluated as the difference of Hamiltonian
expectation values and reads

$Q=\braket{H}_{\rho_{f}}-\braket{H}_{\rho_{i}}\\\
=\sum_{k}\omega_{k}\,\big{[}f(\beta_{2},\omega_{k})-f(\beta_{1},\omega_{k})\big{]}.$
(20)

In the next section, we discuss how to combine the latter two transformations
to implement an Otto cycle and the generalization of these results when the
assumptions of perfect thermalization and quantum adiabaticity are relaxed.

##  IV The quantum Otto cycle

In this section we derive the expression for the work performed by a real
quantum Otto engine, i.e., an engine whose transformations may violate the
assumptions of perfect adiabaticity and/or perfect thermalization. We
consider, as a working substance (medium), a system of coupled fermions,
according to Eq. (5). After a Bogoliubov transformation, such Hamiltonian can
be mapped into the free-quasiparticle model (10), with dispersion relation
$\omega_{k}(\lambda)$. As we shall see in a moment, the only ingredient needed
to assess the performance of such an engine is the quasiparticles distribution
on the evolved system state.

The quantum Otto cycle operates between the two temperatures $\beta_{c}^{-1}$
(cold) and $\beta_{h}^{-1}$ (hot), and consists of four strokes: two
(thermodynamic) adiabatic transformations (i.e., no heat exchange with the
environment) and two isochoric transformations (i.e., no work on the system)
[11]. The system, initially prepared in the thermal state $\rho_{c}$ at
temperature $\beta_{c}^{-1}$, undergoes the following steps:

  1. 1.

$[$Forward adiabatic transformation$]$: The Hamiltonian parameter $\lambda$ is
varied in a time $T$, from $\lambda_{i}$ to $\lambda_{f}$;

  2. 2.

$[$Hot isochoric thermalization$]$: The Hamiltonian $H_{f}\equiv
H(\lambda_{f})$ is fixed and the system is put in contact with the reservoir
at temperature $\beta_{h}^{-1}>\beta_{c}^{-1}$, for a time $\tau$;

  3. 3.

$[$Backward adiabatic transformation$]$: The Hamiltonian parameter $\lambda$
is varied, in a time $T$ from $\lambda_{f}$ to $\lambda_{i}$, thus reversing
point 1;

  4. 4.

$[$Cold isochoric thermalization$]$: The Hamiltonian $H_{i}\equiv
H(\lambda_{i})$ is fixed and the system is put in contact with the reservoir
at temperature $\beta_{c}^{-1}<\beta_{h}^{-1}$ for a time $\tau$.

###  IV.1 Ideal free-fermion cycle

We consider, as a working medium, a quadratic fermionic system weakly coupled
to two thermal reservoirs, as modeled in App. A. By exploiting the results of
Sec. III.2, the work performed and heat exchanged for an ideal Otto cycle
(strictly speaking, perfect thermalization and perfect quantum adiabatic
sweeps are ideally obtained in the limits $T\to\infty$, $\tau\to\infty$) can
be written as:

$\displaystyle Q_{h}^{\text{id}}=$
$\displaystyle\phantom{-}\sum_{k}\omega_{k}(\lambda_{f})\,\big{(}\Delta
f_{k}\big{)}_{hc}\,,$ (21a) $\displaystyle Q_{c}^{\text{id}}=$
$\displaystyle-\sum_{k}\omega_{k}(\lambda_{i})\,\big{(}\Delta
f_{k}\big{)}_{hc}\,,$ (21b) $\displaystyle W^{\text{id}}=$
$\displaystyle\phantom{-}\sum_{k}\big{[}\omega_{k}(\lambda_{f})-\omega_{k}(\lambda_{i})\big{]}\,\big{(}\Delta
f_{k}\big{)}_{hc}\,,\quad$ (21c) where $\big{(}\Delta f_{k}\big{)}_{hc}\equiv
f\big{[}\beta_{h},\omega_{k}(\lambda_{f})\big{]}-f\big{[}\beta_{c},\omega_{k}(\lambda_{i})\big{]}.$
(21d)

Note that, since after a single ideal cycle the system comes back to the same
initial state, we have $\Delta E=0$ and thus
$W^{\text{id}}=Q^{\text{id}}_{h}+Q^{\text{id}}_{c}$ follows from the first
principle of thermodynamics. The efficiency of the cycle can be thus expressed
as:

$\eta^{\text{id}}\equiv\frac{W^{\text{id}}}{Q_{h}^{\text{id}}}=1-\frac{\sum_{k}\omega_{k}(\lambda_{i})\big{(}\Delta
f_{k}\big{)}_{hc}}{\sum_{q}\omega_{q}(\lambda_{f})\big{(}\Delta
f_{q}\big{)}_{hc}}.$ (22)

We can obtain some basic considerations by introducing the ratio

$r_{k}(\lambda_{i},\lambda_{f})\equiv\omega_{k}(\lambda_{i})/\omega_{k}(\lambda_{f}).$
(23)

In fact, the sign of $\left(\Delta f_{k}\right)_{hc}$ depends on the value
assumed by $r_{k}(\lambda_{i},\lambda_{f})$. If
$r_{k}(\lambda_{i},\lambda_{f})>\beta_{h}/\beta_{c}$, then $\left(\Delta
f_{k}\right)_{hc}$ is positive; otherwise it is negative. Since we have
$\omega_{k}(\lambda)>0\hskip 5.69046pt\forall k$, from the sign of
$\left(\Delta f_{k}\right)_{hc}$ we can infer some information on the way the
engine is operating. In particular,

  * •

If $r_{k}(\lambda_{i},\lambda_{f})>\beta_{h}/\beta_{c}\ \forall k$, then
$Q_{h}>0$ and $Q_{c}<0$, thus the engine works either as a thermal accelerator
(if $W<0$) or as a heat engine (if $W>0$);

  * •

If $r_{k}(\lambda_{i},\lambda_{f})>1\ \forall k$, then $Q_{h}>0$, $Q_{c}<0$,
and $W<0$, thus the engine is a thermal accelerator;

  * •

If $\beta_{h}/\beta_{c}<r_{k}(\lambda_{i},\lambda_{f})<1\ \forall k$, then
$Q_{h}>0$, $Q_{c}<0$, and $W>0$, thus the engine is a heat engine;

  * •

If $r_{k}(\lambda_{i},\lambda_{f})<\beta_{h}/\beta_{c}\ \forall k$, then
$Q_{h}<0$ and $Q_{c}>0$, thus the engine is a refrigerator.

These conditions are valid for any system described by a quadratic fermionic
Hamiltonian. It is important to remark that they are sufficient, but not
necessary, conditions [e.g., if
$r_{k}(\lambda_{i},\lambda_{f})>\beta_{h}/\beta_{c}$ only for some $k$, it is
still possible to have an engine operating as a thermal accelerator or as a
heat engine].

###  IV.2 Nonperfect thermalization

We now relax the hypothesis of perfect thermalization, by assuming the system
to be in contact with the two reservoirs for a time $\tau<t^{\text{th}}$,
where $t^{\text{th}}$ is the thermalization time. According to Ref. [42] the
nonperfect thermalized distribution of the quasiparticles reads

$\braket{b^{\dagger}_{k}b_{k}}(\gamma)=f(\beta,\omega_{k})(1-e^{-\mathcal{J}\gamma})+\braket{b^{\dagger}_{k}b_{k}}_{\rho_{0}}e^{-\mathcal{J}\gamma},$
(24)

where we have put $\gamma\equiv 2\tau$ and hereafter we assume
$\mathcal{J}=1$. Note that $\mathcal{J}$ is a parameter which depends on the
bath properties and on the system-bath coupling, as detailed in App. A, and it
obviously affects the time scale of thermalization. The off diagonal
correlators and the anomalous ones $\braket{b_{k}b_{q}}$ remain zero at any
time, provided the initial condition for the Otto cycle is a thermal state
(say, $\rho_{c}$). Equation (24) ensures that convergence to the perfect
thermal behavior is exponential although, strictly speaking, perfect
thermalization cannot be reached at any finite time
($t^{\text{th}}\to\infty$). Thus the ideal case would correspond to the limit
$\tau\to\infty$.

Let us find an expression for the asymptotic distribution of quasiparticles
after $n$ repetitions of the nonperfect thermalized cycle are implemented. We
introduce the diagonal $N\times N$ matrices [43]

$\displaystyle\mathbb{\Theta}_{c(h)}$ $\displaystyle=$
$\displaystyle\text{diag}\big{\\{}f[\beta_{c(h)},\omega_{k}(\lambda_{i(f)})]\big{\\}}_{k=1,\ldots,N}\>,$
(25) $\displaystyle\mathbb{\Gamma}_{c(h)}^{[n]}$ $\displaystyle=$
$\displaystyle\text{diag}\big{\\{}\text{Tr}[b^{\dagger}_{k}b_{k}\,\tilde{\rho}_{c(h)}^{[n]}]\big{\\}}_{k=1,\ldots,N}\>,$
(26)

where $\tilde{\rho}_{c(h)}^{[n]}$ denotes the state at the end of the
nonperfect thermalized stroke with the bath at temperature $\beta_{c(h)}^{-1}$
of the $n$-th cycle repetition. The system is initially prepared in the
thermal state $\rho_{c}$ at temperature $\beta^{-1}_{c}$, i.e.
$\mathbb{\Gamma}^{[0]}_{c}\equiv\mathbb{\Theta}_{c}$. Since the asymptotic
distribution does not depend on the initial condition, this choice does not
affect the generality of the result. Then, at the end of the ideal adiabatic
transformation, the system is put in contact with the reservoir at temperature
$\beta_{h}^{-1}$ and reaches the state $\rho_{h}^{[1]}$, characterized by the
distribution

$\mathbb{\Gamma}_{h}^{[1]}=\mathbb{\Theta}_{h}(1-e^{-\gamma})+\mathbb{\Gamma}^{[0]}_{c}e^{-\gamma}.$
(27a) After the backward ideal adiabatic transformation is implemented, the
system is put in contact with the reservoir at temperature $\beta_{c}^{-1}$
reaching the state $\rho_{c}^{[1]}$, characterized by the distribution
$\mathbb{\Gamma}^{[1]}_{c}=\mathbb{\Theta}_{c}(1-e^{-\gamma})+\mathbb{\Gamma}_{h}^{[1]}e^{-\gamma}.$
(27b)

The two equations (27) can be easily generalized to the $n$-th iteration, to
obtain

$\displaystyle\mathbb{\Gamma}_{c}^{[n]}=\big{(}\mathbb{\Theta}_{c}+e^{-\gamma}\mathbb{\Theta}_{h}\big{)}(1-e^{-\gamma})+\mathbb{\Gamma}_{c}^{[n-1]}e^{-2\gamma},$
(28a)
$\displaystyle\mathbb{\Gamma}_{h}^{[n]}=\big{(}\mathbb{\Theta}_{h}+e^{-\gamma}\mathbb{\Theta}_{c}\big{)}(1-e^{-\gamma})+\mathbb{\Gamma}_{h}^{[n-1]}e^{-2\gamma}.$
(28b)

By imposing the stationary condition that both the $\mathbb{\Gamma}$’s on the
left and right hand side are iteration independent, one gets the stationary
solution

$\displaystyle\mathbb{\Gamma}_{c}^{\infty}=h(\gamma)\big{(}\mathbb{\Theta}_{c}+e^{-\gamma}\mathbb{\Theta}_{h}\big{)},$
(29a)
$\displaystyle\mathbb{\Gamma}_{h}^{\infty}=h(\gamma)\big{(}\mathbb{\Theta}_{h}+e^{-\gamma}\mathbb{\Theta}_{c}\big{)},$
(29b)

where $h(\gamma)=(1+e^{-\gamma})^{-1}$. Being

$\frac{d\big{(}\mathbb{\Gamma}_{i}^{[n]}\big{)}_{kk}}{d\big{(}\mathbb{\Gamma}_{i}^{[n-1]}\big{)}_{kk}}=e^{-2\gamma}<1,$
(30)

with $\big{(}\mathbb{\Gamma}_{i}^{[n]}\big{)}_{kk}$ the $k$-th diagonal
element of the matrix ($i=c,h$), the fixed points in Eq. (29) are stable and
the convergence is exponential in $n$.

Equations (29) suggest that, for the nonperfectly thermalized cycle (n-th),
the heat absorbed and the work performed per cycle are evaluated from Eq. (21)
by substituting the rescaled Fermi functions

$\displaystyle f[\beta_{c},\omega_{k}(\lambda_{i})]$
$\displaystyle\\!\mapsto\\!$ $\displaystyle
h(\gamma)\big{\\{}f[\beta_{c},\omega_{k}(\lambda_{i})]+e^{-\gamma}f[\beta_{h},\omega_{k}(\lambda_{f})]\big{\\}},$
$\displaystyle f[\beta_{h},\omega_{k}(\lambda_{f})]$
$\displaystyle\\!\mapsto\\!$ $\displaystyle
h(\gamma)\big{\\{}f[\beta_{h},\omega_{k}(\lambda_{f})]+e^{-\gamma}f[\beta_{c},\omega_{k}(\lambda_{i})]\big{\\}},$

thus obtaining

$Q_{h(c)}^{\text{n-th}}=g(\gamma)\,Q_{h(c)}^{\text{id}},\qquad
W^{\text{n-th}}=g(\gamma)\,W^{\text{id}},$ (31)

where $g(\gamma)=(1-e^{-\gamma})h(\gamma)=\tanh(\gamma/2)$. Therefore it
follows that, in this case, the efficiency of the heat engine
$\eta^{\text{n-th}}\equiv W^{\text{n-th}}/Q_{h}^{\text{id}}$ remains the same
as for the ideal case [see Eq. (22)].

Figure 1: The value of $p(\delta,\gamma)$ versus the thermalization time
$\gamma=2\tau$, for different values of the quench duration $\delta$ (color
scale in the legend).

If $W>0$ we can define the power of the engine as function of the
thermalization time ($\tau=\gamma/2$) and of the quench time ($T=\delta/2$):

$\mathcal{P}^{\text{n-th}}(\delta,\gamma)\equiv\frac{W^{\text{n-th}}}{\delta+\gamma}=p(\delta,\gamma)\,W^{\text{id}},$
(32)

with

$p(\delta,\gamma)\equiv\frac{\tanh(\gamma/2)}{\delta+\gamma}.$ (33)

In this case we are not concerning about nonadiabatic effects, i.e., Eq. (32)
is valid if we choose $\delta$ such that, according to the features of the
quenches ($\mathbb{H}(\lambda)$, $\lambda_{i}$, $\lambda_{f}$), the quantum
adiabatic theorem approximately holds (see Sec. IV.3). Of course, this value
of $\delta$ is model dependent. In Fig. 1 we show the behavior of
$p(\delta,\gamma)$ versus $\gamma$, for different values of $\delta$ fixed.
Since the work for the ideal case $W^{\text{id}}$ is independent of $\delta$
and $\gamma$, the power is maximum for cycle parameters maximizing
$p(\delta,\gamma).$ Therefore, we can optimize the performance by fixing
$\delta$ (and assuming that it is suitable for adiabatic theorem to hold) and
maximizing with respect to $\gamma$, in order to have some indication on the
best $\gamma$ for any given quench duration. The maximization

$\partial_{\gamma}p(\delta,\gamma)=\frac{1}{\delta+\gamma}\bigg{[}\frac{1}{2\cosh^{2}(\gamma/2)}-p(\delta,\gamma)\bigg{]}=0$
(34)

leads to the stationarity condition

$\frac{e^{\gamma}-e^{-\gamma}}{2}-\gamma=\delta.$ (35)

Let us define $S(\gamma)=\frac{e^{\gamma}-e^{-\gamma}}{2}-\gamma$. For
$\gamma\gg 1$, $S(\gamma)\sim e^{\gamma}/2$. As a consequence,
$\gamma_{\text{max}}\sim\log(2\delta)$, i.e., the maximum scales
logarithmically with the quench duration. In correspondence of the maximum we
have

$p\big{(}\delta,\log(2\delta)\big{)}\sim\frac{1}{\delta}\frac{(2\delta)^{1/2}-(2\delta)^{-1/2}}{(2\delta)^{1/2}+(2\delta)^{-1/2}}.$
(36)

Therefore, for $\delta\to\infty$, we have $p(\delta,\gamma)\sim 1/\delta\to
0$, i.e., no power is produced by the engine. We observe that for $\gamma\ll
1$ we have $S(\gamma)\sim\gamma^{3}/6$, leading to
$\gamma_{\text{max}}\sim(6\delta)^{1/3}$, with $\delta\ll 1$. By substituting
this value, we notice that for fast adiabatic sweeps we have
$p(\delta,\gamma_{\text{max}})\to 0.5$, paying attention that this result
holds for regimes in which even for small $\delta$ quantum adiabaticity is
still valid. In Fig. 2 we show $\gamma_{\text{max}}$ (blue curve) and
$p(\delta,\gamma_{\text{max}})$ (orange curve) vs $\delta$. The dashed lines
are the expected analytical results for large $\delta$.

Figure 2: The value of $\gamma_{\text{max}}$ that maximizes the power output
(blue curve) and the corresponding value of $p(\delta,\gamma_{\text{max}})$
(orange curve) versus $\delta$. Dashed lines denote the asymptotic
large-$\delta$ behaviors analytically evaluated.

###  IV.3 Real adiabatic processes

We now consider a system prepared in the thermal state $\rho_{1}$ with an
Hamiltonian $H_{i}$ characterized by a spectrum $\omega_{k}(\lambda_{i})$. We
implement a transformation $\lambda_{i}\to\lambda_{f}$ in a time $T$ such that
the quantum adiabatic regime is no more valid [41]. As discussed in Sec.
III.2, the work can be evaluated from the first law of thermodynamics.
However, because of the nonperfect quantum adiabaticity of the transformation,
some excitations are generated during the dynamics.

The unitary dynamics is generated by the operator
$U_{ev}(t)=\mathcal{T}\exp\big{[}e^{-i\int_{0}^{t}dsH(s)}\big{]}$, so that

$W=\ \braket{H_{i}}_{\rho_{1}}-\braket{H_{f}}_{\rho_{1}(T)},$ (37)

where $\rho_{1}(T)=U_{\text{ev}}(T)\,\rho_{1}\,U_{\text{ev}}^{\dagger}(T)$
denotes the post-quench state. The first term in the r.h.s. is known, while
the second one can be easily evaluated in the Nambu spinor notation, using Eq.
(6):

$\displaystyle\braket{H_{f}}_{\rho_{1}(T)}=$
$\displaystyle\sum_{i,j}\text{Tr}\big{[}\Psi^{\dagger}_{i}(\mathbb{H}_{f})_{ij}\Psi_{j}\,U_{\text{ev}}(T)\rho_{1}U_{\text{ev}}^{\dagger}(T)\big{]}$
$\displaystyle=$
$\displaystyle\sum_{i,j}\text{Tr}\big{[}\Psi^{H\dagger}_{i}(\mathbb{H}_{f})_{ij}\Psi^{H}_{j}\rho_{1}\big{]},$
(38)

where we adopted the Heisenberg representation and introduced the time-evolved
operators

$\displaystyle\Psi^{H}$ $\displaystyle=U_{\text{ev}}^{\dagger}(T)\Psi
U_{\text{ev}}(T)$
$\displaystyle=\big{(}c_{1}^{H}(T),\dots,c_{N}^{H}(T),c_{1}^{H\dagger}(T),\dots,c_{N}^{H\dagger}(T)\big{)}^{T}.$
(39)

As discussed in Refs. [44, 45], these operators can be expressed in terms of
the $b_{k}$-fermions by implementing the following time-dependent Bogoliubov
transformation

$\Psi^{H}=\mathbb{V}(T)\Phi,$ (40)

where $\mathbb{V}(T)$ is a matrix obeying the following equation

$\partial_{s}\mathbb{V}(s)=-2i\mathbb{H}(s)\mathbb{V}(s),$ (41)

with initial condition $\mathbb{V}(0)=\mathbb{U}(h_{i})$. Substituting this
result in Eq. (38), we get

$\braket{H_{f}}_{\rho_{1}(T)}=\sum_{i,j}\text{Tr}\big{[}\Phi_{i}^{\dagger}\mathbb{\tilde{H}}_{ij}\Phi_{j}\rho_{1}\big{]},$
(42)

being $\mathbb{\tilde{H}}=\mathbb{V}^{\dagger}(T)\mathbb{H}_{f}\mathbb{V}(T)$.
By construction, the state $\rho_{1}$ is diagonal in the $b_{k}$-fermions
(i.e. $\braket{b_{q}^{\dagger}b_{q^{\prime}}}=0$, if $q\neq q^{\prime}$),
therefore we have

$\braket{H_{f}}_{\rho_{1}(T)}=\sum_{k}\tilde{\omega}_{k}\Big{\\{}f\big{[}\beta,\omega_{k}(\lambda_{i})\big{]}-\tfrac{1}{2}\Big{\\}},$
(43)

where $\tilde{\omega}_{k}\equiv\mathbb{\tilde{H}}_{kk}$ (for $1\leq k\leq N$)
and we used the fact that $\mathbb{V}(T)$ is a Bogoliubov transformation and
so it maintains particle-hole symmetry
$(\mathbb{\tilde{H}}_{kk}=-\mathbb{\tilde{H}}_{k+N,k+N})$. Collecting all
these results, the work at the end of the adiabatic stroke reads

$W=\sum_{k}\big{[}\omega_{k}(\lambda_{i})-\tilde{\omega}_{k}(\lambda_{f})\big{]}\Big{\\{}f[\beta,\omega_{k}(\lambda_{i})]-\tfrac{1}{2}\Big{\\}}.$
(44)

From a different perspective, in the Heisenberg picture, we can introduce a
new set of fermions $d_{k}^{(\dagger)}\equiv b_{k}^{H(\dagger)}$, whose Nambu
spinor $\Xi=(d_{1},\dots,d_{N},d_{1}^{\dagger},\dots,d_{N}^{\dagger})^{T}$
relates to the quasiparticles one through the Bogoliubov transformation
$\mathbb{Q}=\mathbb{U}^{\dagger}(\lambda_{f})\mathbb{V}(T)$,
$\Xi=\mathbb{Q}\Phi$. In this way, we can write

$\braket{\Phi^{\dagger}_{k}\Phi_{q}}_{\rho_{1}(T)}=\braket{\Xi_{k}^{\dagger}\Xi_{q}}_{\rho_{1}}=\sum_{i,j=1}^{2N}\mathbb{Q}^{*}_{ki}\braket{\Phi_{i}^{\dagger}\Phi_{j}}_{\rho_{1}}\mathbb{Q}_{qj}.$
(45)

Note that, in the limit of an ideal adiabatic transformation, the unitary
evolution of the eigenvectors of $\mathbb{H}(\lambda_{0})$ [Eq. (41)] is such
that these ones are instantaneous eigenvectors of $\mathbb{H}(t)$, therefore
$\mathbb{V}^{\dagger}(T)\,\mathbb{H}_{f}\,\mathbb{V}(T)=\mathbb{U}^{\dagger}(\lambda_{f})\,\mathbb{H}_{f}\,\mathbb{U}(\lambda_{f})=\mathbb{H^{D}}_{f}$.
So $\mathbb{\tilde{H}}=\mathbb{H^{D}}_{f}$ and Eq. (37) reduces to Eq. (18).

We also point out that, as shown in Ref. [46], imperfect quantum adiabatic
transformations starting from a thermal state excite the system more than the
ideal quantum adiabatic ones. As a consequence, the former are expected to
produce less work than latter case. In fact, we verified such statement
through several numerical simulations, always obtaining
$\tilde{\omega}_{k}<\omega_{k}(\lambda_{f})$, that is $W^{id}>W$.

###  IV.4 Nonperfect thermalization and real adiabatic processes

Let us finally assume nonperfect thermalization and real adiabatic
transformations. The system is prepared in the thermal state at temperature
$\beta_{1}^{-1}$, $\rho_{0}\equiv\rho_{1}$. We can follow a procedure similar
to that of Sec. IV.2, but we have to identify two other intermediate steps in
the dynamics. Moreover, in the case of real adiabatic transformations, the
nondiagonal components of the correlator of the Bogoliubov quasiparticles can
be nonvanishing. To account for this effect, we derive the results using the
Nambu spinors and, in analogy with Eq. (26), define

$\big{(}\mathbb{\Lambda}_{i}^{[n]}\big{)}_{jl}=\braket{\Phi_{j}\Phi^{\dagger}_{l}}_{\rho_{i}^{[n]}},$
(46)

where $\rho_{i}^{[n]}$ represents the state of the system after the $n$th
thermalization with the bath at temperature $\beta_{i}^{-1}$. We also define

$\big{(}\mathbb{\Lambda}_{i,T}^{[n]}\big{)}_{jl}=\braket{\Phi_{j}\Phi^{\dagger}_{l}}_{\rho_{i}^{[n]}(T)},$
(47)

being $\rho_{i}^{[n]}(T)=\mathbb{U}(T)\rho_{i}^{[n]}\mathbb{U}(T)^{\dagger}$
the time evolved density matrix. We notice that

$\displaystyle\braket{\Phi_{j}\Phi^{\dagger}_{l}}_{\rho_{i}^{[n]}(T)}=$
$\displaystyle\text{Tr}\big{[}\Phi_{j}\Phi^{\dagger}_{l}\rho_{1}^{[n]}(T)\big{]}$
(48) $\displaystyle=$
$\displaystyle\text{Tr}\big{[}\Xi_{j}\Xi_{l}^{\dagger}\,\rho_{1}^{[n]}\big{]}$
$\displaystyle=$
$\displaystyle\text{Tr}\sum_{mn}\big{[}\mathbb{Q}_{jm}\Phi_{m}\Phi^{\dagger}_{n}\,\mathbb{Q}^{\dagger}_{nl}\rho_{1}^{[n]}\big{]}$
$\displaystyle=$
$\displaystyle\sum_{mn}\mathbb{Q}_{jm}\text{Tr}\big{[}\Phi_{m}\Phi_{n}^{\dagger}\,\rho_{1}^{[n]}\big{]}\mathbb{Q}^{\dagger}_{nl}$
$\displaystyle=$
$\displaystyle\Big{(}\mathbb{Q}\mathbb{\Lambda}_{i}^{[n]}\mathbb{Q}^{\dagger}\Big{)}_{jl}.$

We introduce the diagonal $2N\times 2N$ matrix

$\mathbb{\Omega}_{c(h)}=\text{diag}\big{\\{}f\left[\beta_{c(h)},-\omega_{k}(\lambda_{i(f)})\right],f\left[\beta_{c(h)},\omega_{k}(\lambda_{i(f)})\right]\big{\\}},$
(49)

with $i=c,h$. We start in the thermal state
$\mathbb{\Lambda}^{[0]}_{1}=\mathbb{\Omega}_{c}$.

After the first adiabatic stroke we have

$\mathbb{\Lambda}^{[1]}_{1,T}=\mathbb{Q}\mathbb{\Lambda}_{1}^{[0]}\mathbb{Q}^{\dagger}$
(50)

Then the system is put in contact with the reservoir at temperature
$\beta_{2}^{-1}$ for a time $\tau$. Differently from Sec. IV.2, in this case
the system at the end of the adiabatic sweep could present some excitations,
so we have to consider also the evolution of nondiagonal elements of the
correlation functions [42]

$\displaystyle\braket{b^{\dagger}_{k}b_{q}}(\tau)$ $\displaystyle=$
$\displaystyle\braket{b^{\dagger}_{k}b_{q}}_{\rho_{0}}e^{i(\omega_{k}-\omega_{q})\tau-\gamma},$
(51a) $\displaystyle\braket{b^{\dagger}_{k}b^{\dagger}_{q}}(\tau)$
$\displaystyle=$
$\displaystyle\braket{b^{\dagger}_{k}b^{\dagger}_{q}}_{\rho_{0}}e^{i(\omega_{k}+\omega_{q})\tau-\gamma},$
(51b) $\displaystyle\braket{b_{k}b_{q}}(\tau)$ $\displaystyle=$
$\displaystyle\braket{b_{k}b_{q}}_{\rho_{0}}e^{-i(\omega_{k}+\omega_{q})\tau-\gamma}.$
(51c)

At the end of this process we have

$\mathbb{\Lambda}^{[1]}_{2}=\mathbb{\Omega}_{h}(1-e^{-\gamma})+\mathbb{\Phi}\mathbb{\Lambda}^{[1]}_{1,T}\mathbb{\Phi}^{\dagger}e^{-\gamma},$
(52)

where

$\mathbb{\Phi}\\!\equiv\\!\text{diag}\big{[}e^{-i\omega_{1}\\!(h_{i})\tau},\ldots,e^{-i\omega_{N}\\!(h_{i})\tau},e^{i\omega_{1}\\!(h_{i})\tau},\ldots,e^{i\omega_{N}\\!(h_{i})\tau}\big{]}.$
(53)

We perform the adiabatic transformation ending with

$\mathbb{\Lambda}^{[1]}_{2,T}=\mathbb{Q^{\prime}}\Lambda_{1}^{[0]}\mathbb{Q^{\prime}}^{\dagger},$
(54)

being $\mathbb{Q^{\prime}}$ the rotation diagonalizing the Hamiltonian after
the backward adiabatic transformation. Finally we connect the system to the
reservoir at temperature $\beta_{1}^{-1}$ obtaining

$\mathbb{\Lambda}^{[1]}_{1}=\mathbb{\Omega}_{c}(1-e^{-\gamma})+\mathbb{\Phi}^{\prime}\mathbb{\Lambda}^{[1]}_{2,T}\mathbb{\Phi}^{\prime\dagger}e^{-\gamma},$
(55)

where $\mathbb{\Phi}^{\prime}$ is the same as $\mathbb{\Phi}$, but with the
final eigenstates $\omega_{k}(h_{f})$. In what follows, to simplify the
notation, we redefine the matrices $\mathbb{Q}$, $\mathbb{Q}^{\prime}$
including $\mathbb{\Phi}$, $\mathbb{\Phi}^{\prime}$ in their definition [i.e.,
$\mathbb{Q}=\mathbb{\Phi}\,\mathbb{U}^{\dagger}(h_{f})\,\mathbb{V}(T)$]. This
does not change the mean value of the Hamiltonian at the end of the adiabatic
sweeps. At the $n$-th iteration we find

$\displaystyle\mathbb{\Lambda}^{[n]}_{1,T}$ $\displaystyle=$
$\displaystyle\mathbb{Q}\mathbb{\Lambda}^{[n-1]}_{1}\mathbb{Q}^{\dagger},$
(56a) $\displaystyle\mathbb{\Lambda}^{[n]}_{2}$ $\displaystyle=$
$\displaystyle\mathbb{\Omega}_{h}(1-e^{-\gamma})+\mathbb{\Lambda}^{[n]}_{1,T}e^{-\gamma},$
(56b) $\displaystyle\mathbb{\Lambda}^{[n]}_{2,T}$ $\displaystyle=$
$\displaystyle\mathbb{Q^{\prime}}\mathbb{\Lambda}^{[n]}_{2}\mathbb{Q^{\prime\dagger}},$
(56c) $\displaystyle\mathbb{\Lambda}^{[n]}_{1}$ $\displaystyle=$
$\displaystyle\mathbb{\Omega}_{c}(1-e^{-\gamma})+\mathbb{\Lambda}^{[n]}_{2,T}e^{-\gamma}.$
(56d)

Substituting we obtain

$\displaystyle\mathbb{\Lambda}^{[n]}_{1,T}$ $\displaystyle=$
$\displaystyle\mathbb{Q}\mathbb{K}_{1}\mathbb{Q}^{\dagger}+e^{-2\gamma}\mathbb{Q}\mathbb{Q^{\prime}}\mathbb{\Lambda}_{1,T}^{[n-1]}\mathbb{Q^{\prime}}^{\dagger}\mathbb{Q}^{\dagger},$
(57a) $\displaystyle\mathbb{\Lambda}^{[n]}_{2}$ $\displaystyle=$
$\displaystyle\mathbb{K}_{2}+e^{-2\gamma}\mathbb{Q}\mathbb{Q^{\prime}}\mathbb{\Lambda}_{2}^{[n-1]}\mathbb{Q^{\prime}}^{\dagger}\mathbb{Q}^{\dagger},$
(57b) $\displaystyle\mathbb{\Lambda}^{[n]}_{2,T}$ $\displaystyle=$
$\displaystyle\mathbb{Q}^{\prime}\mathbb{K}_{2}\mathbb{Q}^{\prime\dagger}+e^{-2\gamma}\mathbb{Q^{\prime}}\mathbb{Q}\mathbb{\Lambda}_{2,T}^{[n-1]}\mathbb{Q}^{\dagger}\mathbb{Q^{\prime}}^{\dagger},$
(57c) $\displaystyle\mathbb{\Lambda}^{[n]}_{1}$ $\displaystyle=$
$\displaystyle\mathbb{K}_{1}+e^{-2\gamma}\mathbb{Q^{\prime}}\mathbb{Q}\mathbb{\Lambda}_{1}^{[n-1]}\mathbb{Q}^{\dagger}\mathbb{Q^{\prime}}^{\dagger},$
(57d)

having introduced

$\displaystyle\mathbb{K}_{1}$ $\displaystyle=$
$\displaystyle(1-e^{-\gamma})(\mathbb{\Omega}_{c}+e^{-\gamma}\mathbb{Q}^{\prime}\mathbb{\Omega}_{h}\mathbb{Q}^{\prime\dagger}),$
(58a) $\displaystyle\mathbb{K}_{2}$ $\displaystyle=$
$\displaystyle(1-e^{-\gamma})(\mathbb{\Omega}_{h}+e^{-\gamma}\mathbb{Q}\mathbb{\Omega}_{c}\mathbb{Q}^{\dagger}).$
(58b)

We can write explicitly the series in the first equation of dynamical systems
Eq. (57), finding

$\displaystyle\mathbb{\Lambda}^{[n+1]}_{1,T}=$
$\displaystyle\sum_{k=0}^{n-1}e^{-2k\gamma}\big{[}(\mathbb{Q}\mathbb{Q^{\prime}})^{k}\mathbb{Q}\mathbb{K}_{1}\mathbb{Q}^{\dagger}(\mathbb{Q^{\prime}}^{\dagger}\mathbb{Q}^{\dagger})^{k}\big{]}$
$\displaystyle+e^{-2n\gamma}\big{[}(\mathbb{Q}\mathbb{Q^{\prime}})^{n}\mathbb{\Lambda}^{[1]}_{i,T}(\mathbb{Q^{\prime}}^{\dagger}\mathbb{Q}^{\dagger})^{n}\big{]}.$
(59)

The asymptotic solution of this equation reads

$\mathbb{\Lambda}^{\infty}_{1,T}=\sum_{k=0}^{\infty}e^{-2k\gamma}\big{[}(\mathbb{Q}\mathbb{Q^{\prime}})^{k}\mathbb{Q}\mathbb{K}_{1}\mathbb{Q}^{\dagger}(\mathbb{Q^{\prime}}^{\dagger}\mathbb{Q}^{\dagger})^{k}\big{]}.$
(60)

Similar expressions can be written for the other matrices.

Before concluding we note that, due to the unitarity of the matrices
$\mathbb{Q},\mathbb{Q^{\prime}}$, and the fact that $\mathbb{\Omega}_{c(h)}$
is a diagonal matrix with real elements in $[0,1]$, we have:

$\Big{\lvert}\big{[}(\mathbb{Q}\mathbb{Q^{\prime}})^{k}\mathbb{Q}\mathbb{K}_{1}\mathbb{Q}^{\dagger}(\mathbb{Q^{\prime}}^{\dagger}\mathbb{Q}^{\dagger})^{k}\big{]}_{ij}\Big{\rvert}<\left(1-e^{-2\gamma}\right).$
(61)

This means that the terms in Eq. (60) decay exponentially to zero with $k$. In
particular, the modulus of each term is smaller than

$\xi(k)\equiv e^{-2k\gamma}\left(1-e^{-2\gamma}\right).$ (62)

So we can safely approximate the asymptotic value of
$\mathbb{\Lambda}_{1}^{\infty}$ with the first $K$ terms of the series. The
convergence of $\mathbb{\Lambda}^{n}_{1,T}$ to
$\mathbb{\Lambda}^{\infty}_{1,T}$ is exponential in $n$,

$\left(\mathbb{\Lambda}^{[n+1]}_{1,T}-\mathbb{\Lambda}^{[n]}_{1,T}\right)_{ij}\propto
e^{-2\gamma(n-1)}.$ (63)

In what follows we present the results for the performance of an Ising quantum
Otto cycle when nonadiabatic transformations are implemented.

##  V The real Ising quantum Otto cycle

In this section we focus on a quantum Otto engine, whose working medium is
constituted by a chain of $N$ spin-$1/2$ systems interacting through the well
known one dimensional transverse-field Ising Hamiltonian:

$H^{\text{spin}}(t)=-J\sum_{i=1}^{N-1}\sigma_{i}^{x}\sigma_{i+1}^{x}-h(t)\sum_{i=1}^{N}\sigma_{i}^{z},$
(64)

where $\sigma^{\alpha}_{i}$ are the Pauli matrices ($\alpha=x,y,z$) acting on
the $i$-th site. To avoid the presence of degenerate energy eigenstates, we
assume open boundary conditions. Hereafter we also set $J=1$ as the energy
scale of the system and set $h>0$ without loss of generality. The model (64)
exhibits a zero-temperature quantum phase transition at $h_{c}=1$, from a
paramagnetic ($h>h_{c}$) to a ferromagnetic ($h<h_{c}$) phase, after
spontaneously breaking the $\mathbb{Z}_{2}$ symmetry that rotates spins by a
$\pi$-angle around the $z$ axis [47, 48].

The Hamiltonian (64) can be mapped into a quadratic fermionic model as the one
in Eq. (5), through a Jordan-Wigner transformation

$\sigma_{i}^{-}=K_{i}c_{i}\,,\qquad\text{with}\quad
K_{i}=\Pi_{j=1}^{i-1}\sigma_{j}^{z}\,,$ (65)

where $\sigma_{i}^{\pm}=(\sigma_{i}^{x}\pm i\sigma_{i}^{y})/2$ are the
raising/lowering operators for the $i$-th spin. The mapped fermionic model is
usually referred to as the Kitaev chain, which maintains the $\mathbb{Z}_{2}$
symmetry of the Ising chain, in the form of fermionic parity: the Ising spin-
spin coupling transforms into nearest-neighbor hopping and $p$-wave pairing
terms, while the transverse field is mapped into a chemical potential term
[49].

Figure 3: Sketch of the Ising quantum Otto cycle consisting of four strokes:
(1) forward adiabatic (increasing of the transverse field), (2) hot isochoric
(thermalization with the hot reservoir), (3) backward adiabatic (decreasing of
the transverse field), (4) cold isochoric (thermalization with the cold
reservoir). Note that, in the thermalization strokes, each spin is assumed to
be coupled to a single bath at the same temperature $\beta_{h/c}$, in such a
way that thermalization to that nominal temperature $\beta_{h/c}$ is
guaranteed [42].

One can engineer an Ising quantum Otto cycle by implementing the four
transformations (strokes) discussed in Sec. IV, with the transverse field
$h_{i}$ taking the role of the Hamiltonian control parameter $\lambda$ [34].
In general, for the engine being useful, the order of the transformations must
be chosen carefully. As detailed in App. C, here we implement the following
procedure, starting from a configuration in which the working substance (the
system) is supposed to be in the thermal state at the (cold) temperature
$\beta_{c}^{-1}$:

  1. 1.

The system is decoupled from the external environment and the transverse field
$h(t)$ is increased linearly in a finite time $T=\delta/2$, such that
$h(t)=h_{i}+(h_{f}-h_{i})t/T$ (with $t\in[0,T]$), where $h_{f}>h_{i}$;

  2. 2.

The transverse field is kept fixed at $h(t)=h_{f}$, while the system is put in
contact with a hot reservoir at $\beta_{h}^{-1}$ for a finite time
$\tau=\gamma/2$ and let equilibrate with it;

  3. 3.

The system is again decoupled from the external environment and the transverse
field $h(t)$ is decreased linearly in a finite time $T=\delta/2$, such that
$h(t)=h_{f}+(h_{i}-h_{f})t/T$ (with $t\in[0,T]$);

  4. 4.

The transverse field is kept fixed at $h(t)=h_{i}$, while the system is put in
contact with a cold reservoir at $\beta_{c}^{-1}$ for a finite time
$\tau=\gamma/2$ and let equilibrate with it.

A sketch of this cycle is shown in Fig. 3, where the thermalization is
achieved by coupling each spin to a single local bath at the same temperature,
within the Lindblad framework (see App. A for details).

In Ref. [34] it has been shown that, for this choice of the quench parameters
with $h_{f}>h_{i}$ (details in App. C), in the ideal case, this engine can
operate both as a heat engine and as a refrigerator, and that the emergence of
criticality in the many-body medium may play an important role in determining
the performance. Below we relax the assumption of ideal thermodynamic
transformations (which would require an infinite operational time and
represents an unrealistic condition corresponding to zero power) and discuss
the cases where the quantum adiabaticity regime is not occurring in the
strokes 1 and 3.

###  V.1 Perfect thermalization

Figure 4: (a) The work performed by an imperfect Ising quantum Otto engine
with $N=50$ spins operating between two reservoirs at temperature
$\beta_{h}^{-1}=1$ and $\beta_{c}^{-1}=0.5$. The quench amplitude is
$h_{f}-h_{i}=0.5$. Different colors refer to various quench velocities. We
observe a more pronounced dependence on $\delta$ for quenches crossing the
critical point (dotted line). (b) The corresponding power output, for the same
parameters of (a). The plot is obtained by assuming a thermalization time
$\gamma=5$. We note that in the ferromagnetic phase the power exhibits a
nonmonotonic behavior in $\delta$.

Following the results of Sec. IV.3, it is straightforward to derive the
expressions for the heat exchanged and the work performed under the assumption
of perfect thermalization strokes (2 and 4), but quantum nonadiabatic
transformations:

$\displaystyle Q^{\text{real}}_{h}=$
$\displaystyle\phantom{-}\sum_{k}\omega_{k}(h_{f})\tilde{f}\big{[}\beta_{h},\omega_{k}(h_{f})\big{]}\\!-\\!\tilde{\omega}_{k}(h_{f})\tilde{f}\big{[}\beta_{c},\omega_{k}(h_{i})\big{]},$
$\displaystyle Q^{\text{real}}_{c}=$
$\displaystyle-\\!\sum_{k}\tilde{\omega}_{k}(h_{i})\tilde{f}\big{[}\beta_{h},\omega_{k}(h_{f})\big{]}\\!-\\!\omega_{k}(h_{i})\tilde{f}\big{[}\beta_{c},\omega_{k}(h_{i})\big{]},$
$\displaystyle W^{\text{real}}=$
$\displaystyle\phantom{-}\sum_{k}\Big{\\{}\left[\omega_{k}(h_{f})-\tilde{\omega}_{k}(h_{i})\right]\tilde{f}\big{[}\beta_{h},\omega_{k}(h_{f})\big{]}$
$\displaystyle\qquad\;-\big{[}\tilde{\omega}_{k}(h_{f})-\omega_{k}(h_{i})\big{]}\tilde{f}\big{[}\beta_{c},\omega_{k}(h_{i})\big{]}\Big{\\}},$
(66)

where $\tilde{f}\big{[}\beta_{h(c)},\omega_{k}(h_{f(i)})\big{]}\equiv
f\big{[}\beta_{h(c)},\omega_{k}(h_{f(i)})\big{]}-1/2$ and
$\tilde{\omega}_{k}^{f(i)}=\tilde{\mathbb{H}}^{f(i)}_{kk}$, as defined in Sec.
IV.3. Note that $W^{\text{real}}=Q^{\text{real}}_{h}+Q^{\text{real}}_{c}$,
since the perfect thermalization strokes reset the system state to the thermal
one. We can thus define the power of the engine as

$\mathcal{P}^{\text{real}}=\frac{W^{\text{real}}}{\delta+\gamma}.$ (67)

In Fig. 4 we show numerical results for the work (a) and the power output (b)
of an engine made of $N=50$ spins and operating between temperatures
$\beta_{h}^{-1}=1$ and $\beta_{c}^{-1}=0.5$, for a fixed quench amplitude
$h_{f}-h_{i}=0.5$. Unless specified otherwise, in the following all
simulations have been performed with these parameters. Data are plotted
against the initial transverse field $h_{i}$, while the various curves refer
to different quench durations $\delta$. To evaluate the power, we fix a
thermalization time of $\gamma=5$ so that, according to Eq. (62),
$\xi(1)\approx 5\times 10^{-5}$. Under the assumption of $\mathcal{J}=1$, this
is sufficient to achieve a nearly perfect thermalization for the working
medium, therefore here the isochoric transformations can be considered ideal,
although occurring in a finite time $\gamma$.

Figure 5: Power (yellow curve) and work (green curve), evaluated at
$h_{i}=0.76$, versus $\delta$ (cf. data in Fig. 4). While the work grows
monotonically with the quench duration, the power exhibits a nonmonotonic
behavior.

An important emerging feature is the sensitivity of the work on the quench
duration $\delta$, for quenches occurring in proximity to the critical point
$h=1$ (marked by the dotted line), such that $h_{i}\lesssim 1$ and
$h_{f}\gtrsim 1$. For the selected range of $\delta\in[10,200]$, this marked
dependence disappears when operating entirely in the paramagnetic phase
($h_{i}>1$), due to considerably larger values of the relevant energy gaps; in
fact, a phenomenology similar to what we observe for $h$ crossing the critical
point can be recovered when considering much faster quenches ($\delta\ll 10$),
becoming comparable with the spectral gap (we checked that, for $h_{i}\approx
1.5$, one should take $\delta\sim 1$ to detect a clear sensitivity in $\delta$
on the scale of the figure).

Focusing on the region $h_{i}\lesssim 1$ and for sufficiently large values of
$\delta\gtrsim 12$, we can identify the value $h_{i}\approx 0.76$ as the
position of the so-called critical peak, defined as the maximum value assumed
by $W^{\text{real}}$ or by $\mathcal{P}^{\text{real}}$ in the region $h_{i}<1$
[34]. Besides that, both the initial field $h_{i}$ and the quench velocity
$\delta$ may modify the operation regime of the engine ($W<0$, for $h_{i}$
smaller than a given threshold which depends on $\delta$): in fact, the energy
wasted in excitations could prevent the engine from producing work, i.e., the
system dissipates so much energy that one would have to perform work on it, in
order to operate the cycle. Moreover, while curves for the work in Fig. 4(a)
are monotonic in the quench duration, those in Fig. 4(b) exhibit different
scalings with $\delta$, suggesting the existence of an optimal working point
which maximizes the power delivered by the engine.

Figure 6: Work performed by the Ising Otto engine versus the thermalization
time for different quench durations. Data are for $N=50$ spins, $h_{i}=0.76$,
$\beta_{h}^{-1}=1$, and $\beta_{c}^{-1}=0.5$.

To further investigate this feature, in Fig. 5 we plot the value of the work
(green curve) and of the power (yellow curve) evaluated at the critical peak
($h_{i}=0.76$) as a function of $\delta$. These data display a clear
nonmonotonicity of the power, contrary to the monotonic behavior of the work.
Different parameters (i.e., quench amplitudes, reservoir temperatures, and
sizes of the medium) do not qualitatively alter this scenario, although the
work and power output at the optimal working point generally depends on the
size of the energy gaps that are relevant to the unitary dynamics [50].

Figure 7: Contour plot of the power per spin for the nonideal quantum Ising
Otto engine, as a function of $\delta$ and $\gamma$, for $h_{i}=0.76$,
$\beta_{h}^{-1}=1$, and $\beta_{c}^{-1}=0.5$. We adopted a logarithmic color
scale, to better appreciate variations in the parameter space. The three
panels refer to different lengths of the quantum Ising model (i.e., of the
mapped Kitaev chain): $N=20$ (a), $N=30$ (b), and $N=50$ (c). The blue dots
indicate, for any $\delta$, the value of $\gamma$ at which the maximum of the
power occurs. The red dot locates the global maximum of the power output. Data
in the bottom-left corners, below the white dots, are not significant, since
they correspond to a parameter range where the system is not operating as a
heat engine.

###  V.2 Nonperfect thermalization

We now relax the hypothesis of perfect thermalization. In this case, in the
asymptotic regime, we find that

$\displaystyle Q_{c}^{\text{real, n-th}}$ $\displaystyle=$
$\displaystyle\text{Tr}\big{[}\,\mathbb{H}_{d}^{f}\left(\mathbb{\Lambda}_{1,T}-\mathbb{\Lambda}_{2}\right)\big{]},$
$\displaystyle Q_{h}^{\text{real, n-th}}$ $\displaystyle=$
$\displaystyle\text{Tr}\big{[}\mathbb{H}_{d}^{i}\left(\mathbb{\Lambda}_{2,T}-\mathbb{\Lambda}_{1}\right)\big{]},$
(68) $\displaystyle W^{\text{real, n-th}}$ $\displaystyle=$ $\displaystyle
Q_{c}+Q_{h},$

where the trace is over the $2N$ degrees of freedom of the Bogoliubov
matrices.

Equations (V.2) for the stationary heat and the work cannot be expressed
analytically. However, we can exploit Eq. (60) to obtain an approximate
expression for the correlation matrix that is suitable for numerical
computation. In particular we can stop at order $K$, knowing we are making an
error $\varepsilon<8N\xi(K+1)$ for work and $\varepsilon<4N\xi(K+1)$ for heat,
where $\xi(K+1)$ is defined in Eq. (62). This approximation breaks down when
the modulus of the work computational derived is of order
$N\xi\left(K+1\right)$, in this case we cannot be sure that we are neglecting
relevant terms, while for $\left|W\right|\gg 8N\xi(K+1)$ the approximation
works.

In Fig. 6 we show the work performed by the engine at the critical peak [51],
as a function of the thermalization time $\gamma$ and for different quench
durations $\delta$. Curves are all monotonic, both in $\delta$ and in
$\gamma$, and exhibit a fast (exponential) convergence to the asymptotic value
obtained for the ideal thermalization (although the thermalization time scale
$\gamma$ obviously depends on the system-bath coupling details — see App. A).
This implies a nonmonotonic behavior of the power with $\delta$ and $\gamma$,
similarly to the case of quantum nonadiabatic transformations and perfect
thermalization discussed in Sec. V.2. Note also that imperfect thermalization
strokes, especially for very fast quenches, can influence the operation modes
as well.

The corresponding power per spin is plotted in Fig. 7 (color scale) as a
function of $\delta$ and $\gamma$, for various sizes of the working medium
(the mapped Kitaev chain): $N=20$ (a), $N=30$ (b), and $N=50$ (c). In each
panel, the blue dots locate the position of the local maximum for fixed
$\delta$, while the red dot indicates the global one. The data in the bottom-
left corner, below the white points, refer to a region of the parameter space
in which the Otto cycle is not operating as a heat engine, therefore it should
not be considered. A comparison between the three panels reveals that the
system size basically affects the absolute value of the power output, which is
maximized for small $N$. This is in agreement with the fact that, at
criticality, the energy gap closes with the system size and thus the work
extraction for large $N$ is affected by fast quenches. The functional form of
the power (and, consequently, the location of the maxima) with $\delta$ and
$\gamma$ seems to be independent of $N$. Finally, the edge of the operation
modes, when increasing the system size, slightly moves towards smaller
$\delta$ and larger $\gamma$.

We notice that, even though some nonmonotonicity in $\gamma$ emerges, this is
not as relevant as for the case discussed in Fig. 5, suggesting that the role
played by the thermalization time in the optimization of the power output is
marginal. This is not surprising when considering that the thermalization
timescale is smaller than the time needed for implementing the adiabatic
sweep. In our analysis we assumed $\mathcal{J}=1$; of course, a different
choice (corresponding to an alternative modeling of the thermal reservoirs and
of their coupling to the working substance, see App. A) would have lead to
different thermalization timescales, thus quantitatively modifying the
analysis of the performances of the engine.

We mention that the analysis performed so far for the heat engine could be
repeated for the refrigerator mode, as well. In the perfect thermalized case,
we expect the coefficient of performance (i.e., the ratio between the heat
absorbed from the cold bath $Q_{c}$ and the work done on the system $|W|$) to
be maximized by ideal quantum adiabatic processes, as for the heat engine. On
the other hand, for $Q_{c}$ the analysis is more complicated: in fact, it does
not depend on the velocity of the quench in which $h$ increases, but only on
the velocity of the other quench. Therefore the heat absorbed in time from the
cold bath should have a monotonic behavior with the duration of the quench in
which $h$ increases (better results for sudden quenches) and a nonmonotonic
behavior with the duration of the quench in which $h$ decreases.

##  VI Conclusions

In this work we analyzed the performances of an imperfect quantum many-body
Otto engine based on a free-fermion medium. After recalling the basic concepts
of quantum thermodynamics and some properties of quadratic fermionic systems,
we briefly discussed the isothermal, the isochoric, and the adiabatic
(unitary) transformation in this context. We derived analytic expressions for
heat and the work exchanged when such transformations are performed under the
hypotheses of perfect thermalization and quantum adiabaticity. Then we
considered isochoric and unitary transformations in which these hypotheses of
ideality are violated (i.e., nonperfect thermalization due to a finite-time
contact with an external thermal bath, or nonperfect quantum adiabaticity due
to a finite speed in the variation of the Hamiltonian control parameter).
Applying our machinery to the case of a nonperfect quantum Otto cycle, we have
extended results of Ref. [34], finding that, while the work performed is
always maximized in the ideal engine, the power output exhibits a nonmonotonic
behavior with the thermalization time $\tau=\gamma/2$ and the quench time
$T=\delta/2$, suggesting the existence of an optimal working point in the
$(\tau,T)$ parameter space.

A possible extension of this research is directed towards a systematic
characterization of quantum heat engines based on the spectral properties of
the Hamiltonian. As an example, it could be interesting to predict
analytically the operational mode and the performances expected, as well as to
derive analytic conditions to identify the optimal working points, once the
Hamiltonian details are provided. One could also study quenches with a time
pattern different than linear [52, 53] (e.g., in such a way that the variation
of the Hamiltonian parameters are faster far from criticality and slower near
to it) or more complicated approaches favoring adiabatic dynamics [54]: this
kind of schemes may help to increase the power of the engine and thus change
its optimal working point.

As a future perspective, it would be also tempting to generalize our results
to working media composed of truly interacting quantum many-body systems: this
may represent an important step to understand the role of interactions in the
work extraction process. Moreover, even though here we mainly focused on the
average values of work and power, additional important information can be
obtained by characterizing the fluctuations of the observed quantities (both
for the ideal and for the real scenario), in order to determine whether these
kind of engines can be practically useful for the work production [55, 56, 57,
58]. Finally, one could investigate alternative strategies to enhance the
engine performance by considering a working medium coupled to nonMarkovian
[59] or even nonthermal baths [60, 5] (e.g., some reservoir suitably
engineered to invert the population of the system).

We hope that our analysis will serve as a guidance to any possible
experimental realizations of many-body quantum engines, that are far from
being ideal. Proposals to experimentally realize the Kitaev chain with arrays
of coupled superconducting quantum dots have been recently put forward [61,
62], indicating a possible strategy to build up such kind of devices.

###### Acknowledgements.

We acknowledge useful discussions with M. Campisi and G. Chiriacò. GP
acknowledges financial support from PNRR MUR Project PE0000023-NQSTI.

##  Appendix A Thermalization of free-fermion system

To perform the two isochoric transformations, we suppose that the working
medium is in contact with an environment made of $N$ identical thermal baths,
described by a continuum of fermionic harmonic oscillators that ensures proper
thermalization in the asymptotic long coupling-time limit. The environment
Hamiltonian reads

$H_{\text{env}}=\sum_{i}\int
dq\,\varepsilon_{i}(q)\,\kappa^{\dagger}_{i}(q)\,\kappa_{i}(q),$ (69)

with $\varepsilon_{i}(q)\geq 0$ and $\kappa_{i}^{(\dagger)}(q)$ being the
annihilation (creation) fermionic operators. We assume the baths to be
independent, therefore the reduced density matrix of the environment can be
written in the factorized form

$\rho_{\text{env}}=\bigotimes_{i}\rho_{\text{{\hbox{i}}-th bath}},$ (70)

where $\rho_{\text{$i$-th bath}}$ is the thermal density matrix describing the
$i$-th bath at temperature $\beta^{-1}$ (that corresponds to $\beta_{c}^{-1}$
and $\beta_{h}^{-1}$ for the thermalization with the cold and hot reservoir,
respectively).

We assume each of the above baths to be coupled to a single site of the chain
through a quadratic, factorizable, Hamiltonian

$H_{\text{int}}=\sum_{i}\int
dq\,g_{i}(q)\,\big{(}c_{i}+c_{i}^{\dagger}\big{)}\,\big{[}\kappa_{i}(q)+\kappa^{\dagger}_{i}(q)\big{]}.$
(71)

Then we introduce the density of states of the baths

$\mathcal{J}_{i}(\omega)\equiv\pi\int
dq\,|g_{i}(q)|^{2}\,\delta\big{[}\omega-\varepsilon_{i}(k)\big{]}\,.$ (72)

Assuming a very large bath bandwidth with respect to the system frequencies,
one has $\mathcal{J}_{i}(\omega)\simeq\mathcal{J}_{i}$. By tracing out the
environmental degrees of freedom and imposing the Born-Markov approximation
for the baths, we derive, in the energy eigenbasis, the microscopic Lindblad
master equation for the reduced density matrix of the system, whose stationary
solution is a thermal state at temperature $T$ [63, 42]. From this equation
and under the hypothesis of no degeneracies in the spectrum, we can
analytically derive the expressions for the Bogoliubov quasiparticles
correlator at the time $\tau$ of the thermalization process:

$\displaystyle\braket{b_{k}^{\dagger}b_{k}}(\tau)=$
$\displaystyle\,f(\beta,\omega_{k})\big{(}1-e^{-2\mathcal{J}\tau}\big{)}+\braket{b_{k}^{\dagger}b_{k}}(0)\,e^{-2\mathcal{J}\tau},$
$\displaystyle\braket{b_{k}^{\dagger}b_{q}}(\tau)=$
$\displaystyle\braket{b^{\dagger}_{k}b_{q}}_{\rho_{0}}e^{i(\omega_{k}-\omega_{q})\tau-2\mathcal{J}t},$
$\displaystyle\braket{b^{\dagger}_{k}b^{\dagger}_{q}}(\tau)=$
$\displaystyle\braket{b^{\dagger}_{k}b^{\dagger}_{q}}_{\rho_{0}}e^{i(\omega_{k}+\omega_{q})\tau-2\mathcal{J}t},$
$\displaystyle\braket{b_{k}b_{q}}(\tau)=$
$\displaystyle\braket{b_{k}b_{q}}_{\rho_{0}}e^{-i(\omega_{k}+\omega_{q})\tau-2\mathcal{J}t}.$
(73)

Note that we assumed the same density of state $\mathcal{J}_{i}=\mathcal{J}$
for all the local baths. The first expression coincides with Eq. (24) with
$\gamma=2\tau$, while the other expressions are the ones reported in Eqs. (51)
(in the latter, $\mathcal{J}=1$).

##  Appendix B Thermodynamic limit of the Ideal Ising Otto engine

In the thermodynamic limit $N\to+\infty$, it is possible to obtain a simple
analytical expression for the eigenvalues of $\mathbb{H}^{\text{spin}}$:

$\omega_{k}(h)=2\sqrt{1+h^{2}-2h\cos(k)},~{}k\in(0,\pi].$ (74)

In this limit, even though the adiabatic theorem is not valid, because of the
critical closure of the gap for $h=1$, it is instructive to study the behavior
of the heat and the work per spin, defined as

$q_{h(c)}=\lim_{N\to+\infty}Q_{h(c)}/N,\hskip
8.5359ptw=\lim_{N\to+\infty}W/N.$ (75)

By substituting $\sum_{k=1}^{N}\mapsto\frac{N}{\pi}\int_{0}^{\pi}dk$ in Eqs.
(21), we find

$\displaystyle q_{h}^{\text{id}}$ $\displaystyle=$
$\displaystyle\phantom{-}\frac{2}{\pi}\int_{0}^{\pi}{dk\hskip
2.84544pt\omega_{k}(h_{f})\,\big{(}\Delta f_{k}\big{)}_{hc}}\,,$ (76a)
$\displaystyle q_{c}^{\text{id}}$ $\displaystyle=$
$\displaystyle-\frac{2}{\pi}\int_{0}^{\pi}{dk\hskip
2.84544pt\omega_{k}(h_{i})\,\big{(}\Delta f_{k}\big{)}_{hc}}\,,$ (76b)
$\displaystyle w^{\text{id}}$ $\displaystyle=$
$\displaystyle\phantom{-}\frac{2}{\pi}\int_{0}^{\pi}{dk\hskip
2.84544pt\big{[}\omega_{k}(h_{f})-\omega_{k}(h_{i})\big{]}\big{(}\Delta
f_{k}\big{)}_{hc}}\,,\qquad$ (76c)

where $w^{\text{id}}=q^{\text{id}}_{h}+q^{\text{id}}_{c}$, while the ratio in
Eq. (23) becomes

$r_{k}\left(h_{i},h_{f}\right)=\sqrt{\frac{1+h^{2}_{i}-2h_{i}\cos(k)}{1+h^{2}_{f}-2h_{f}\cos(k)}}\,.$
(77)

Figure 8: Operating modes of the ideal Ising Otto engine in the thermodynamic
limit ($N\to\infty$) and for fixed $\beta^{-1}_{h}=1$, $h_{f}-h_{i}=0.5$. The
color code indicates the operating mode: orange stands for a thermal
accelerator, green for a heat engine, red for a heater and blue for a
refrigerator. All the four modes are possible, as shown numerically for finite
$N$, in Ref. [34]. The dashed and dashed-dotted lines represent, respectively,
the equations $|1-h_{i}|/|1-h_{f}|=\beta_{h}/\beta_{c}$ and
$(1+h_{i})/(1+h_{f})=\beta_{h}/\beta_{c}$. The intersection point is the
Carnot point. The area in between the white portions of such lines is a
parameter subspace in which the system can only operate as a refrigerator; the
area under the black portions of such lines is where the system can operate
either as a heat engine or a heat accelerator.

Now, taking the derivative of the square of the ratio,

$\frac{\partial r_{k}^{2}(h_{i},h_{f})}{\partial
k}=\frac{2\left(h_{f}h_{i}-1\right)\left(h_{f}-h_{i}\right)\sin(k)}{\big{[}1+h_{f}^{2}-2h_{f}\cos(k)\big{]}^{2}}$
(78)

and using the fact that $h_{i}<h_{f}$ for construction, we obtain that, if
$h_{i}>h_{f}^{-1}$ then $\partial_{k}r_{k}^{2}(h_{i},h_{f})>0$, otherwise
$\partial_{k}r_{k}^{2}(h_{i},h_{f})<0$. So, if $h_{i}>h_{f}^{-1}$, then
$r_{k}(h_{i},h_{f})$ is an increasing function for $k\in(0,\pi]$ and plugging
the dispersion relation in • ‣ IV.1 we have:

  * •

for $r_{k\to
0^{+}}(h_{i},h_{f})\equiv\frac{|1-h_{i}|}{|1-h_{f}|}>\beta_{h}/\beta_{c}$ the
system is a heat engine;

  * •

for $r_{k=\pi}(h_{i},h_{f})\equiv\frac{1+h_{i}}{1+h_{f}}<\beta_{h}/\beta_{c}$
the system is a refrigerator.

On the other hand, if $h_{i}<h_{f}^{-1}$, then $r_{\kappa}(h_{i},h_{f})$ is a
decreasing function for $k\in(0,\pi]$, thus we have:

  * •

for $r_{k=\pi}(h_{i},h_{f})\equiv\frac{1+h_{i}}{1+h_{f}}>\beta_{h}/\beta_{c}$
the system could be a heat engine or a thermal accelerator;

  * •

for $r_{k\to
0^{+}}(h_{i},h_{f})\equiv\frac{|1-h_{i}|}{|1-h_{f}|}<\beta_{h}/\beta_{c}$ the
system is a refrigerator.

The remaining case $h_{i}=h_{f}^{-1}$ is a special point, since in that case
$r_{k}$ is constant ($r_{k}=h_{i}$). This allows us to say that, for
$\beta_{h}/\beta_{c}<h_{i}<1$, the system is a heat engine, while for
$h_{i}<\beta_{h}/\beta_{c}$ it performs as a refrigerator. The point
$h_{i}=h_{f}^{-1}=\beta_{h}/\beta_{c}$ is called a Carnot point, since in that
case no work is produced nor the heat is exchanged, and the efficiency of the
engine saturates the Carnot bound.

Figure 8 shows a contour plot of the modes in which the Ising Otto engine in
the thermodynamic limit can operate, as a function of the initial field
strength $h_{i}$ and of the cold reservoir temperature $\beta^{-1}_{c}$, while
we fix $\beta^{-1}_{h}=1$ and $h_{f}-h_{i}=0.5$ (see caption). This
complements the numerical analysis performed in Ref. [34] for finite values of
$N$, since it is obtained with a fully analytic approach and by directly
accessing the thermodynamic limit.

##  Appendix C Work in the Ising Otto engine

In this section we comment on the order of the adiabatic strokes that are
necessary to obtain an Otto engine that can produce useful work. For classical
engines, in which the work is univocally associated to a variation of
mechanical energy, the situation is rather intuitive: the cold working medium
is adiabatically compressed and then warmed up in order to further increase
the energy that, during the adiabatic expansion, is released, thus producing
work. This order can be quite easily extended to the case of few-body quantum
engines (e.g., a collection of noninteracting qubits, as in Ref. [10]), in
which the adiabatic compression (expansion) corresponds to a compression
(expansion) of the spectrum of the system.

The situation becomes more involved when considering many-body quantum
systems, since is not immediate to find a way to stretch or compress the full
eigenspectrum. We thus call a compression (an expansion) any quench that
increases (decreases) the internal energy of the system. In quantum heat
engines the work can be associated to a variation of some physical quantity
that is not necessarily the mechanical energy. As a consequence, one does not
need to keep the order of the transformations unaltered, and, on the contrary,
it can also be reversed.

Figure 9: (a): Behavior of $M(h)$ in function of the transverse field $h$ for
$\beta^{-1}=0.1$ (blue curve) and $\beta^{-1}=0.5$ (red curve). (b): Inset of
panel (a) for quench from $h_{i}=0.25$ to $h_{f}=0.5$. The arrows represent in
which direction the transverse field must be varied at any temperature to
perform, by implementing sequentially the two transformations, the total work
represented by the shaded area. (c): The same of panel (b) for quench from
$h_{i}=1$ to $h_{f}=1.5$.

As an example, let us consider the Ising quantum engine. In this case, the
work is associated to a change of the magnetization $M(h)$ (see Refs. [64,
65]). Hence, the work performed by the Ising chain during a quasistatic
isothermal transformation at a fixed temperature $\beta^{-1}$ can be evaluated
according to Eq. (16):

$W=\int_{h_{i}}^{h_{f}}dh\,M(h),\qquad\text{with}\
M(h)\equiv-\frac{\partial\mathcal{F}}{\partial h}.$ (79)

In Fig. 9(a) we show the function $M(h)$ for two different temperatures
$\beta^{-1}=0.1$ (blue curve) and $\beta^{-1}=0.5$ (red curve) versus the
transverse field. If $h$ is varied from smaller to larger values, the work
performed during the transformation is the area beneath such curve. If the
transformation is performed backward, the work acquires a minus sign, compared
as before.

We now consider a sequence of four transformations implementing a Carnot
cycle: a forward quasistatic variation of the transverse field while the
system is in thermal equilibrium at $\beta^{-1}_{1}$, a thermalization with
the bath at temperature $\beta^{-1}_{2}$, a backward quasistatic variation of
$h$ while the system is in thermal equilibrium at $\beta^{-1}_{2}$, and a
thermalization with the bath at temperature $\beta_{1}^{-1}$. The total work
of such an engine is the sum of the works performed along the two adiabatic
strokes. Therefore, the temperatures $\beta_{1}^{-1}$ and $\beta_{2}^{-1}$
must be chosen in order to follow the curves clockwise. As an example, if the
Hamiltonian parameter is varied from $h_{i}=0.25$ to $h_{f}=0.5$, as in Fig.
9(b), the only possible choice is $\beta^{-1}_{1}=0.5$ and
$\beta_{2}^{-1}=0.1$. Conversely, considering a quench from $h_{i}=1$ to
$h_{f}=1.5$, as in Fig. 9(c), one should choose $\beta_{1}^{-1}=0.1$ and
$\beta_{2}^{-1}=0.5$. In both figures, the shaded area represents the work
performed by the engine.

For a generic fermionic quadratic system as the one in Eq. (5), the order to
perform the cycle can be argued by looking at the spectrum. If
$\omega^{\prime}_{k}(\lambda)\geq 0$,
$\forall\lambda\in[\lambda_{i},\lambda_{f}]$, $\forall k$, (here $\lambda$
denotes the control parameter) the work output is positive (produced by the
system) for increasing $\lambda$, since, in the interval
$\lambda\in[\lambda_{i},\lambda_{f}]$,

$-\frac{\partial\mathcal{F}}{\partial\lambda}=\sum_{k}\omega^{\prime}_{k}(\lambda)\Big{[}\tfrac{1}{2}-f\big{[}\beta,\omega_{k}(\lambda)\big{]}\Big{]}>0.$
(80)

Beside this, we have that

$-\frac{\partial}{\partial\beta}\frac{\partial\mathcal{F}}{\partial\lambda}=\sum_{k}\omega^{\prime}_{k}(\lambda)\Big{[}\omega_{k}(\lambda)e^{\beta\omega_{k}(\lambda)}f^{2}\big{[}\beta,\omega_{k}(\lambda)\big{]}\Big{]}>0,$
(81)

where we used the fact that $\omega_{k}(\lambda)>0$ $\forall k$. Equation (81)
proves that the work produced is larger at small temperatures. As a
consequence, the order of the cycle should follow that in Fig. 9(c), in which
the forward adiabatic transformation is performed at $\beta_{c}$. If, on the
other hand, $\omega^{\prime}_{k}(\lambda)\leq 0$,
$\forall\lambda\in[\lambda_{i},\lambda_{f}]$, $\forall k$, using the same
arguments as in the latter case, we conclude that, to produce work, the order
of the transformations must be reversed.

Even though, to simplify the discussion, here we have only considered
quasistatic transformations (instead of the adiabatic processes for the Otto
cycle, we adopted the isotherm processes for the Carnot cycle), we can use the
above results to gain some insights for the adiabatic transformation of the
quantum Otto Ising cycle, as well. In particular, since $W^{\text{ad}}\leq
W^{\text{iso}}$, a necessary condition for the Otto cycle to operate as a heat
engine is that the Carnot cycle operating between the same temperatures
($\beta_{c}^{-1}$, $\beta_{h}^{-1}$) produces work. For example, in Fig. 4, to
obtain an operative heat engine for quenches of the Ising chain from
$h_{i}\in[0.5,2]$ to $h_{f}=h_{i}+0.5$, we have chosen as the initial
temperature $\beta_{1}^{-1}\equiv\beta_{c}^{-1}=0.5$.

## References

  * Alicki [1979] R. Alicki, The quantum open system as a model of the heat engine, J. Phys. A: Math. Gen. 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).
  * Scovil and Schulz-DuBois [1959] H. E. Scovil and E. O. Schulz-DuBois, Three-level masers as heat engines, Phys. Rev. Lett. 2, 262 (1959).
  * Quan  _et al._ [2007] H. T. Quan, Y.-x. Liu, C. P. Sun, and F. Nori, Quantum thermodynamic cycles and quantum heat engines, Phys. Rev. E 76, 031105 (2007).
  * Rossnagel  _et al._ [2014] J. Rossnagel, O. Abah, F. Schmidt-Kaler, K. Singer, and E. Lutz, Nanoscale heat engine beyond the carnot limit, Phys. Rev. Lett. 112, 030602 (2014).
  * Hartmann  _et al._ [2020a] A. Hartmann, V. Mukherjee, W. Niedenzu, and W. Lechner, Many-body quantum heat engines with shortcuts to adiabaticity, Phys. Rev. Res. 2, 023145 (2020a).
  * Kosloff and Rezek [2017] R. Kosloff and Y. Rezek, The quantum harmonic Otto cycle, Entropy 19, 136 (2017).
  * Bender  _et al._ [2000] C. M. Bender, D. C. Brody, and B. K. Meister, Quantum mechanical Carnot engine, J. Phys. A: Math. Gen. 33, 4427 (2000).
  * Bender  _et al._ [2002] C. M. Bender, D. C. Brody, and B. K. Meister, Entropy and temperature of a quantum carnot engine, Proc. R. Soc. A: Math. Phys. Eng. Sci. 458, 1519 (2002).
  * Solfanelli  _et al._ [2020] A. Solfanelli, M. Falsetti, and M. Campisi, Nonadiabatic single-qubit quantum Otto engine, Phys. Rev. B 101, 054513 (2020).
  * Henrich  _et al._ [2007] M. J. Henrich, F. Rempp, and G. Mahler, Quantum thermodynamic Otto machines: A spin-system approach, Eur. Phys. J.: Spec. Top. 151, 157 (2007).
  * Uzdin and Kosloff [2014] R. Uzdin and R. Kosloff, Universal features in the efficiency at maximal work of hot quantum Otto engines, Eur. Phys. Lett. 108, 40001 (2014).
  * Leggio and Antezza [2016] B. Leggio and M. Antezza, Otto engine beyond its standard quantum limit, Phys. Rev. E 93, 022122 (2016).
  * Rossnagel  _et al._ [2016] J. Rossnagel, S. Dawkins, K. 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. Commun. 10, 202 (2019).
  * von Lindenfels  _et al._ [2019] D. von Lindenfels, O. Gräb, C. T. Schmiegelow, V. Kaushal, J. Schulz, M. T. Mitchison, J. Goold, F. Schmidt-Kaler, and U. G. Poschinger, Spin heat engine coupled to a harmonic-oscillator flywheel, Phys. Rev. Lett. 123, 080602 (2019).
  * Peterson  _et al._ [2019] J. P. Peterson, T. B. Batalhñ ao, 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).
  * Sheng  _et al._ [2021] J. Sheng, C. Yang, and H. Wu, Realization of a coupled-mode heat engine with cavity-mediated nanoresonators, Sci. Adv. 7, eabl7740 (2021).
  * Bouton  _et al._ [2021] Q. Bouton, J. Nettersheim, S. Burgardt, D. Adam, E. Lutz, and A. Widera, A quantum heat engine driven by atomic collisions, Nat. Commun. 12, 2063 (2021).
  * Myers  _et al._ [2022] N. M. Myers, O. Abah, and S. Deffner, Quantum thermodynamic devices: From theoretical proposals to experimental reality, AVS Quantum Science 4, 027101 (2022).
  * Mukherjee and Divakaran [2021] V. Mukherjee and U. Divakaran, Many-body quantum thermal machines, J. Phys.: Condens. Matter 33, 454001 (2021).
  * Cangemi  _et al._ [2023] L. M. Cangemi, C. Bhadra, and A. Levy,  (2023), arXiv:2302.00726 [quant-ph] .
  * Jaramillo  _et al._ [2016] J. Jaramillo, M. Beau, and A. del Campo, Quantum supremacy of many-particle thermal machines, New J. Phys. 18, 075019 (2016).
  * Bengtsson  _et al._ [2018] J. Bengtsson, M. N. Tengstrand, A. Wacker, P. Samuelsson, M. Ueda, H. Linke, and S. M. Reimann, Quantum Szilard engine with attractively interacting bosons, Phys. Rev. Lett. 120, 100601 (2018).
  * Chen  _et al._ [2019] Y.-Y. Chen, G. Watanabe, Y.-C. Yu, X.-W. Guan, and A. del Campo, An interaction-driven many-particle quantum heat engine and its universal behavior, npj Quantum Inf. 5, 88 (2019).
  * Fogarty and Busch [2020] T. Fogarty and T. Busch, A many-body heat engine at criticality, Quantum Sci. Technol. 6, 015003 (2020).
  * Carollo  _et al._ [2020] F. Carollo, F. M. Gambetta, K. Brandner, J. P. Garrahan, and I. Lesanovsky, Nonequilibrium quantum many-body Rydberg atom engine, Phys. Rev. Lett. 124, 170602 (2020).
  * Boubakour  _et al._ [2023] M. Boubakour, T. Fogarty, and T. Busch, Interaction-enhanced quantum heat engine, Phys. Rev. Res. 5, 013088 (2023).
  * Watson and Kheruntsyan [2023] R. S. Watson and K. V. Kheruntsyan,  (2023), arXiv:2308.05266 [cond-mat] .
  * 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).
  * Hartmann  _et al._ [2020b] A. Hartmann, V. Mukherjee, G. B. Mbeng, W. Niedenzu, and W. Lechner, Multi-spin counter-diabatic driving in many-body quantum Otto refrigerators, Quantum 4, 377 (2020b).
  * Wang [2138] Q. Wang, Performance of quantum heat engines under the influence of long-range interactions, Phys. Rev. E 102, 012138 (012138).
  * Revathy  _et al._ [2020] B. S. Revathy, V. Mukherjee, U. Divakaran, and A. del Campo, Universal finite-time thermodynamics of many-body quantum machines from Kibble-Zurek scaling, Phys. Rev. Res. 2, 043247 (2020).
  * Piccitto  _et al._ [2022] G. Piccitto, M. Campisi, and D. Rossini, The ising critical quantum otto engine, New J. Phys. 24, 103023 (2022).
  * Solfanelli  _et al._ [2023] A. Solfanelli, G. Giachetti, M. Campisi, S. Ruffo, and N. Defenu, Quantum heat engine with long-range advantages, New J. Phys. 25, 033030 (2023).
  * Williamson and Davis [2023] L. A. Williamson and M. J. Davis, Many-body enhancement in a spin-chain quantum heat engine, Phys. Rev. B 109, 024310 (2023).
  * Campisi and Fazio [2016] M. Campisi and R. Fazio, The power of a critical heat engine, Nat. Commun. 7, 11895 (2016).
  * Del Grosso  _et al._ [2022] N. F. Del Grosso, F. C. Lombardo, F. D. Mazzitelli, and P. I. Villar, Quantum Otto cycle in a superconducting cavity in the nonadiabatic regime, Phys. Rev. A 105, 022202 (2022).
  * foo [a] Below, when not specified, we will always assume summation indexes running from $1$ to $N$. (a).
  * Sakurai and Napolitano [2020] J. J. Sakurai and J. Napolitano,  _Modern Quantum Mechanics, 3rd edition_ (Cambridge University Press, 2020).
  * foo [b] We remind that a transformation can be considered approximately adiabatic if the variation of the system is carried on a time $\delta\gg|\lambda_{i}-\lambda_{f}|/\Delta^{2}$, where $\Delta$ quantifies the relevant energy scale of the model. (b).
  * D’Abbruzzo and Rossini [2021] A. D’Abbruzzo and D. Rossini, Self-consistent microscopic derivation of Markovian master equations for open quadratic quantum systems, Phys. Rev. A 103, 052209 (2021).
  * foo [c] These matrices are diagonal, hence it is possible to work with a single vector containing the diagonal elements of $\mathbb{\Theta}_{c(h)}$ and $\mathbb{\Gamma}^{[n]}_{c(h)}$. However, we prefer to adopt the matrix formalism to simplify the generalization of these calculations in the case of a cycle with real adiabatic transformations. (c).
  * Mbeng  _et al._ [2020] G. B. Mbeng, A. Russomanno, and G. E. Santoro,  (2020), arXiv:2009.09208 [quant-ph] .
  * Piccitto  _et al._ [2023] G. Piccitto, A. Russomanno, and D. Rossini, Entanglement dynamics with string measurement operators, SciPost Phys. Core 6, 078 (2023).
  * Allahverdyan  _et al._ [2004] A. E. Allahverdyan, R. Balian, and T. M. Nieuwenhuizen, Maximal work extraction from finite quantum systems, Europhys. Lett. 67, 565 (2004).
  * Sachdev [2011] S. Sachdev,  _Quantum Phase Transitions, 2nd edition_ (Cambridge University Press, 2011).
  * Rossini and Vicari [2021] D. Rossini and E. Vicari, Coherent and dissipative dynamics at quantum phase transitions, Phys. Rep. 936, 1–110 (2021).
  * Kitaev [2001] A. Y. Kitaev, Unpaired Majorana fermions in quantum wires, Phys.-Usp. 44, 131–136 (2001).
  * foo [d] In the paramegnetic region ($h_{i}>1$) and for the same set of parameters used in Fig. 5, one is however able to achieve a significantly larger power output, since smaller values of $\delta$ are sufficient for the engine to deliver a comparable work. (d).
  * foo [e] We checked that the qualitative behavior of the work and power output is quite insensitive to the specific value of $h_{i}<1$, therefore we are not showing data for other initial conditions. (e).
  * Sen  _et al._ [2008] D. Sen, K. Sengupta, and S. Mondal, Defect production in nonlinear quench across a quantum critical point, Phys. Rev. Lett. 101, 016806 (2008).
  * Barankov and Polkovnikov [2008] R. Barankov and A. Polkovnikov, Optimal nonlinear passage through a quantum critical point, Phys. Rev. Lett. 101, 076801 (2008).
  * Guéry-Odelin  _et al._ [2019] D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Shortcuts to adiabaticity: Concepts, methods, and applications, Rev. Mod. Phys. 91, 045001 (2019).
  * G. and Masanes [2016] R. J. G. and L. Masanes, Work extraction from quantum systems with bounded fluctuations in work, Nat. Commun. 7, 13511 (2016).
  * Holubec and Ryabov [2017] V. Holubec and A. Ryabov, Work and power fluctuations in a critical heat engine, Phys. Rev. E 96, 030102(R) (2017).
  * Denzler and Lutz [2020] T. Denzler and E. Lutz, Efficiency fluctuations of a quantum heat engine, Phys. Rev. Res. 2, 032062(R) (2020).
  * Denzler and Lutz [2021] T. Denzler and E. Lutz, Efficiency large deviation function of quantum heat engines, New J. Phys. 23, 075003 (2021).
  * Wiedmann  _et al._ [2020] M. Wiedmann, J. T. Stockburger, and J. Ankerhold, Non-markovian dynamics of a quantum heat engine: out-of-equilibrium operation and thermal coupling control, New J. Phys. 22, 033007 (2020).
  * 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).
  * Sau and Das Sarma [2012] J. D. Sau and S. Das Sarma, Realizing a robust practical Majorana chain in a quantum-dot-superconductor linear array, Nat. Commun. 3, 964 (2012).
  * Dvir and et al. [2023] T. Dvir and et al., Realization of a minimal Kitaev chain in coupled quantum dots, Nature 614, 445 (2023).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Dorner  _et al._ [2012] R. Dorner, J. Goold, C. Cormick, M. Paternostro, and V. Vedral, Emergent thermodynamics in a quenched quantum many-body system, Phys. Rev. Lett. 109, 160601 (2012).
  * Fusco  _et al._ [2014] L. Fusco, S. Pigeon, T. J. Apollaro, A. Xuereb, L. Mazzola, M. Campisi, A. Ferraro, M. Paternostro, and G. De Chiara, Assessing the nonequilibrium thermodynamics in a quenched quantum many-body system via single projective measurements, Phys. Rev. X 4, 031029 (2014).