Work and entropy production in generalised Gibbs ensembles

We say that an observable A equilibrates if, after some relaxation time, its expectation value is for most times the same $\langle A(t)\rangle \simeq \mathrm{Tr}(A{{\rm{\Omega }}}_{\mathrm{TA}})$ as the expectation value of the infinite time average

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{TA}}(\rho ,H):= \mathop{\mathrm{lim}}\limits_{T\to \infty }\displaystyle \frac{1}{T}{\int }_{0}^{T}{{\rm{e}}}^{-{\rm{i}}{Ht}}\,\rho \,{{\rm{e}}}^{{\rm{i}}{Ht}}{\rm{d}}t,\end{eqnarray} \tag{ 1 }$$

of an initial state ρ of a system described by a Hamiltonian H. A simple calculation shows that the time averaged state corresponds to the de-phased state in the Hamiltonian eigenbasis and for this reason is often called diagonal ensemble. More explicitly, given the distinct energies of the Hamiltonian $\{{E}_{k}\}$ and the projectors onto their corresponding eigenspaces P_k, the time averaged state reads

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{TA}}(\rho ,H)=\displaystyle \sum _{k}{P}_{k}\rho {P}_{k}.\end{eqnarray} \tag{ 2 }$$

The time averaged state corresponds to the maximum entropy state given all the conserved quantities [42]. This observation turns the principle of maximum entropy introduced by Jaynes [43, 44] into a consequence of the quantum dynamics. The principle of maximum entropy states that the probability distribution which best represents the current state of knowledge of the system is the one with largest entropy given the conserved quantities of the system; this principle will be crucial to define our equilibration models.

Although relaxation towards the time averaged state has been proven under very general and naturally fulfilled conditions [23–25, 39], in practice, the diagonal ensemble cannot be used as an equilibration model due to its inefficiency. The description of the equilibrium state by the diagonal ensemble requires the specification of as many conserved quantities as the dimension of the Hilbert space, which scales exponentially in the system size. It is therefore in principle not even possible to save all the data in a computer for a large interacting many-body system, let alone compute the infinite time average efficiently.

In practice, the characterisation of the equilibrium state can in many instances be done by specifying only a few quantities, e.g., the temperature and the chemical potential. The most relevant and common such situation is the canonical ensemble or the Gibbs state, for which only the temperature, or equivalently the energy per particle of the initial state ρ, has to be specified,

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{Gibbs}}(\rho ,H)=\displaystyle \frac{{{\rm{e}}}^{-\beta H}}{Z},\end{eqnarray} \tag{ 3 }$$

where ρ is the state of the system before undergoing the equilibration process, $Z=\mathrm{Tr}({{\rm{e}}}^{-\beta H})$ is the partition function and the inverse temperature $\beta \gt 0$ is fixed by imposing that $\mathrm{Tr}(H{{\rm{\Omega }}}_{\mathrm{Gibbs}})=\mathrm{Tr}(H\rho )$.

For generic, non-integrable models, the thermal state is expected to be indistinguishable from the time averaged state under very mild assumptions which relate to conditions on eigenstates of the Hamiltonian [22, 27, 45] and on the energy distribution of the initial state [46, 47]. While dynamical thermalisation in this sense has not yet been rigorously proven, it is highly plausible, and it can be connected to typicality arguments [48, 49]. The generality of these conditions explains why the canonical ensemble is the corner-stone of the standard thermodynamics. Nevertheless, there are known instances of systems that do not thermalise. One central aim of this work is to study how thermodynamic protocols are modified when the Gibbs ensemble is not a good equilibration model and does not satisfactorily describe the equilibrium state of the system.

Examples of systems which do not fully thermalise to Gibbs states are constituted by integrable systems. The infinite-time averaged states are not well described by the Gibbs ensemble because of the existence of (quasi) local integrals of motion, i.e. conserved quantities Q_i, that retain information about the initial state over an infinite amount of time. Instead, there is strong evidence that they can be well-described by the so-called GGE defined as

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GGE}}(\rho ,H,\{{Q}_{i}\})\propto {{\rm{e}}}^{-\beta H+{\displaystyle \sum }_{j=1}^{q}{\lambda }_{j}{Q}_{j}},\end{eqnarray} \tag{ 4 }$$

where the generalised chemical potential ${\lambda }_{j}$ is a Lagrange multiplier associated with the specific conserved quantity Q_j, $j=1,\ldots ,q$, such that its expectation value is the same as the one of the initial state

$$\begin{eqnarray}&&\mathrm{Tr}({{\rm{\Omega }}}_{\mathrm{GGE}}(\rho ,H,\{{Q}_{j}\})\ {Q}_{k})=\mathrm{Tr}(\rho {Q}_{k}),\end{eqnarray} \tag{ 5 }$$

for each $k=1,\ldots ,q$. The GGE can be understood as an interpolation between the diagonal and the canonical ensembles. The diagonal ensemble maximises the von Neumann entropy $S(\rho )=-\mathrm{Tr}(\rho \mathrm{log}\rho )$ given all the conserved quantities (CQ). The Gibbs ensemble maximises the von Neumann entropy considering only the energy as a conserved quantity. The GGE is situated in between. For a given state ρ and a set of operators (conserved quantities) $\{{Q}_{i}\}$, it is natural to define the set of states compatible with the values the conserved quantities

$$\begin{eqnarray}&&{ \mathcal E }(\rho ,\{{Q}_{i}\}):= \{\sigma | \mathrm{Tr}(\rho {Q}_{i})=\mathrm{Tr}(\sigma {Q}_{i})\}.\end{eqnarray} \tag{ 6 }$$

The GGE is the state that maximises the von Neumann entropy within ${ \mathcal E }(\rho ,\{{Q}_{i}\})$. From this perspective, the ensembles introduced so far can be summarised as

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{TA}}:= {\mathrm{argmax}}_{\sigma \in { \mathcal E }(\rho ,\{{\rm{all}}{\rm{CQ}}\})}S(\sigma ),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GGE}}(\rho ,H,\{{Q}_{i}\}):= {\mathrm{argmax}}_{\sigma \in { \mathcal E }(\rho ,\{H,{Q}_{i}\})}S(\sigma ),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{Gibbs}}(\rho ,H):= {\mathrm{argmax}}_{\sigma \in { \mathcal E }(\rho ,\{H\})}S(\sigma ).\end{eqnarray} \tag{ 9 }$$

A relevant question in the construction of GGEs is how the conserved quantities have to be chosen, which is discussed in appendix A. In general, there is a certain degree of ambiguity of what constants of motion to pick in order to arrive at the appropriate GGE. This discussion is not relevant for the general study pursued in this work, however. It is the aim of this work to study the thermodynamical behaviour of the GGE in full generality, hence we will not have to make any precise assumption about the conserved quantities, unless it is explicitly specified.

To illustrate the above considerations, let us consider a quadratic Hamiltonian of fermions in a one dimensional lattice

$$\begin{eqnarray}&&{H}^{(0)}=\displaystyle \sum _{i=1}^{n}{\epsilon }_{i}{a}_{i}^{\dagger }{a}_{i}+g\displaystyle \sum _{i=1}^{n-1}({a}_{i}^{\dagger }{a}_{i+1}+{a}_{i+1}^{\dagger }{a}_{i}),\end{eqnarray} \tag{ 10 }$$

where n is the total number of sites and a_i (${a}_{i}^{\dagger }$) are the creation (annihilation) operators at the i-site which satisfy the fermionic anti-commutation relations

$$\begin{eqnarray}&&\{{a}_{i},{a}_{j}^{\dagger }\}={\delta }_{i,j},\{{a}_{i},{a}_{j}\}=\{{a}_{i}^{\dagger },{a}_{j}^{\dagger }\}=0.\end{eqnarray} \tag{ 11 }$$

We would like to study how an initially out of equilibrium state relaxes to equilibrium and see that the Gibbs ensemble fails to describe the equilibrium state.

The initial state of the system is taken to be in thermal equilibrium, ${\rho }^{(0)}={{\rm{e}}}^{-\beta {H}^{(0)}}/{ \mathcal Z }$. A quench is then performed to a new Hamiltonian ${H}^{(1)}$,

$$\begin{eqnarray}&&{H}^{(0)}\mapsto {H}^{(1)},\end{eqnarray} \tag{ 12 }$$

in which the energy of the first fermion is modified, ${H}^{(1)}={H}^{(0)}+{\rm{\Delta }}{a}_{1}^{\dagger }{a}_{1}$. After the quench, the population of the first fermion evolves in time $t\gt 0$ as

$$\begin{eqnarray}&&{n}_{1}(t)=\mathrm{Tr}({a}_{1}^{\dagger }{a}_{1}\rho (t))\end{eqnarray} \tag{ 13 }$$

with $\rho (t)={{\rm{e}}}^{-{\rm{i}}{H}^{(1)}t}\,\rho (0){{\rm{e}}}^{{\rm{i}}{H}^{(1)}t}$. As the Hamiltonian is quadratic, it is a problem involving free fermions and can be numerically simulated for very long times and system sizes (see appendix G.1).

In figure 1, we plot the time evolution of the occupation of the first site ${n}_{1}(t)$. As expected, we see that after some relaxation time t, ${n}_{1}(t)$ equilibrates to the value predicted by the GGE—which is relatively far from the one given by the Gibbs equilibration model. The situation described in this example, a quench and the characterisation of the equilibrium state, is extensively studied in the literature, see for a recent review [27]. In order to study thermodynamic processes in which many quenches and equilibrations are performed, it will be necessary to promote the suitability of effective descriptions in terms of GGE states for equilibration processes beyond a single quench.

Zoom In Zoom Out Reset image size

Consider a sequence of changes of the Hamiltonian, as defined by a list of $N+1$ Hamiltonians, ${H}^{(m)}$, where $m=0,1,\ldots ,N$ denotes the step in the protocol and ${H}^{(0)}$ is the initial Hamiltonian. These Hamiltonian transformations ${H}^{(m-1)}\mapsto {H}^{(m)}$ are considered to be quenches, in the sense that they are performed sufficiently fast such that the state of the system ρ is unchanged. Let us denote the time at which the quench ${H}^{(m-1)}\mapsto {H}^{(m)}$ is performed by t_m with ${t}_{m}\lt {t}_{m+1}$ for all m. After a quench, the system evolves under the Hamiltonian ${H}^{(m)}$ for a time ${t}_{m+1}-{t}_{m}$ until a new quench ${H}^{(m)}\mapsto {H}^{(m+1)}$ is performed at time ${t}_{m+1}$. This time interval is taken to be much longer than the equilibration time such that the system can be considered to be in equilibrium. The exact state of the system $\rho (t)$ when m quenches have taken place (${t}_{m}\lt t\lt {t}_{m+1}$) is given by,

$$\begin{eqnarray}&&\rho (t)={{\rm{e}}}^{-{\rm{i}}(t-{t}_{m}){H}^{(m)}}\rho ({t}_{m}){{\rm{e}}}^{{\rm{i}}(t-{t}_{m}){H}^{(m)}},\end{eqnarray} \tag{ 14 }$$

where $\rho ({t}_{m})$ is the state of the system at $t={t}_{m}$ when the Hamiltonian ${H}^{(m)}$ starts to dictate the evolution. The state $\rho ({t}_{m})$ is given by the recursive expression

$$\begin{eqnarray}&&\rho ({t}_{k})={{\rm{e}}}^{-{\rm{i}}({t}_{k}-{t}_{k-1}){H}^{(k-1)}}\rho ({t}_{k-1}){{\rm{e}}}^{{\rm{i}}({t}_{k}-{t}_{k-1}){H}^{(k-1)}},\end{eqnarray} \tag{ 15 }$$

with $\rho ({t}_{0})$ the initial state and k = 1,...,m.

Now, our aim is to construct an effective description of the whole evolution of ρ, in such a way that the state after the mth quench and its posterior equilibration, $\rho (t)$, can be described by an appropriate equilibrium state. We denote such equilibrium state that approximates the real state after m quenches, $\rho (t)$, as ${\omega }_{(\cdots )}^{(m)}$ where $(\cdots )$ is the place holder for one of the three models of equilibration: time-average (TA), GGE or Gibbs. The effective description of (14) is then built in a recursive way as follows,

$$\begin{eqnarray}{\omega }_{\mathrm{TA}}^{(m)} & = & {{\rm{\Omega }}}_{\mathrm{TA}}({\omega }_{\mathrm{TA}}^{(m-1)},{H}^{(m)}),\\ {\omega }_{\mathrm{GGE}}^{(m)} & = & {{\rm{\Omega }}}_{\mathrm{GGE}}({\omega }_{\mathrm{GGE}}^{(m-1)},{H}^{(m)},\{{Q}_{i}^{(m)}\}),\\ {\omega }_{\mathrm{Gibbs}}^{(m)} & = & {{\rm{\Omega }}}_{\mathrm{Gibbs}}({\omega }_{\mathrm{Gibbs}}^{(m-1)},{H}^{(m)}).\end{eqnarray} \tag{ 16 }$$

Here, ${\omega }_{(\cdots )}^{(0)}=\rho ({t}_{0})$ is the intial state, before any quench or evolution has taken place. Note that, when constructing the GGE description, the set of conserved quantities $\{{Q}_{i}^{(m)}\}$ changes for every Hamiltonian ${H}^{(m)}$, as well as the Lagrange multipliers $\{{\lambda }_{j}^{(m)}\}{}_{j=1}^{q}$, or simply the inverse temperature ${\beta }^{(m)}$ in the case of equilibration to the Gibbs ensemble.

In order to provide a motivation and interpretation of equation (16), together with the implicit assumptions that come into play, let us illustrate it with a simple example. Suppose a system initially in state $\rho (0)$ and with Hamiltonian ${H}^{(0)}$. At time t₁, we perform a first quench ${H}^{(0)}\mapsto {H}^{(1)}$ and let the system evolve under ${H}^{(1)};$ at time t₂ we perform second quench ${H}^{(1)}\mapsto {H}^{(2)}$ and let the system evolve under ${H}^{(2)}$ until it equilibrates at time t. For both evolutions, we now consider effective descriptions in terms of GGE states. After the evolution under ${H}^{(1)}$ and immediately before performing the second quench, the system is exactly described by $\rho ({t}_{2})$ as given by equation (15). For a set of conserved quantities $\{{Q}_{i}^{(1)}\}$, the corresponding GGE equilibrium state is given by,

$$\begin{eqnarray}&&{\omega }_{\mathrm{GGE}}^{(1)}={{\rm{\Omega }}}_{\mathrm{GGE}}(\rho ({t}_{1}),{H}^{(1)},\{{Q}_{i}^{(1)}\})\simeq \rho ({t}_{2}),\end{eqnarray} \tag{ 17 }$$

where the symbol '≃' means in this context that the average value of relevant observables is well approximated by ${\omega }_{\mathrm{GGE}}^{(1)}$, that is

$$\begin{eqnarray}&&\mathrm{Tr}(A\rho ({t}_{2}))\simeq \mathrm{Tr}(A{\omega }_{\mathrm{GGE}}^{(1)}).\end{eqnarray} \tag{ 18 }$$

Now, when describing the equilibrium state after the second quench, one can simply apply the same recipe. That is, the state $\rho ({t}^{(1)})$ is the initial state when the evolution under ${H}^{(2)}$ starts. Then, assuming that the new conserved quantities $\{{Q}_{i}^{(2)}\}{}_{i}$ are chosen appropriately and applying the same reasoning one obtains an approximation by taking

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GGE}}(\rho ({t}_{2}),{H}^{(2)},\{{Q}_{i}^{(2)}\})\simeq \rho (t),\end{eqnarray} \tag{ 19 }$$

with t longer than the t₂ plus the subsequent equilibration time. Importantly, note that this effective description is not efficient, in the sense that it requires keeping track of the exact state $\rho ({t}_{2})$ to obtain the equilibrium state at time t. If this is extended to N quenches, having to keep track of the exact evolution until the $(N-1)$th quench is as demanding as keeping track of the whole exact evolution over the process. It is here when the effective description (16) becomes handy, as it can be constructed by keeping track of the value of the conserved quantities only. First of all, coming back to the first evolution, note that by applying (16) with m = 1 we recover (17), i.e., the standard result for single quenches. Now, in order to construct the GGE state corresponding to $\rho (t)$, we assume that the conserved quantities $\{{Q}_{i}^{(2)}\}$ are within the set of physically relevant observables A in (18). That is, we assume that

$$\begin{eqnarray}&&\mathrm{Tr}({Q}_{i}^{(2)}\rho ({t}_{2}))\simeq \mathrm{Tr}({Q}_{i}^{(2)}{\omega }_{\mathrm{GGE}}^{(1)})\end{eqnarray} \tag{ 20 }$$

for all i. In this way, in order to obtain the equilibrium GGE ensemble after the second quench, it is not necessary to keep track of the exact state $\rho ({t}_{2})$, but one can simply use ${\omega }_{\mathrm{GGE}}^{(1)}$ instead. Using (20) we then obtain,

$$\begin{eqnarray}{\omega }_{\mathrm{GGE}}^{(2)} & := & {{\rm{\Omega }}}_{\mathrm{GGE}}({\omega }_{\mathrm{GGE}}^{(1)},{H}^{(2)},\{{Q}_{i}^{(2)}\})\\ & \simeq & {{\rm{\Omega }}}_{\mathrm{GGE}}(\rho ({t}_{2}),{H}^{(2)},\{{Q}_{i}^{(2)}\})\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\simeq \,\rho (t).\end{eqnarray} \tag{ 22 }$$

Extending the same reasoning to the case of N quenches and other models of equilibration other than the GGE, we arrive to an effective description of the form (16).

In the rest of this work we will always use the effective description (16) for the full process consisting on a sequence of quenches and equilibrations. We do not claim by this that this model will accurately describe the real dynamics of any system or protocol, and indeed we explicitly leave here as an open question to identify for which Hamiltonians and conserved quantities condition (20) is satisfied for each quench. Nonetheless, and in exactly the same way as equilibration to the Gibbs state is assumed in the usual scenario in thermodynamics, we will assume that equilibration to statistical ensembles of the form (16) occurs over any protocol, so that we can tackle questions about entropy production and work extraction.

To examine the validity of our model, we provide a numerical comparison in section 5.3 of the real exact evolution and the model of equation (16) for different initial states and protocols, given a system of free fermions. We will see for several examples that the model predicts with great accuracy the amount of work that is extracted in a protocol involving a sequence of quenches.

Another interesting issue is to understand how the accuracy of our model decreases when perfect equilibration is not reached, as a consequence of the time of equilibration being too short. In particular, a crucial point is to understand whether the error scales badly with the number of quenches, so that the model (16) becomes progressively worse as N increases. In section 6.3.4, we investigate the accuracy of our model for a free fermionic system for increasing N and a constant time for the whole process. The results suggest that the error becomes independent on the number of quenches that are implemented while it decreases with the total time of the protocol, so that the model (16) remains a good description of free fermions for finite equilibration times and any number of quenches.

Concatenations of quenches and equilibrations constitute a framework to describe thermodynamic processes—see, e.g., [12, 17, 50]. Within this framework, work is associated with the input energy under quenches, whereas heat is associated with the exchange of energy under equilibration processes. At the level of average quantities, the work cost of a single quench, ${H}^{(m-1)}\mapsto {H}^{(m)}$, reads

$$\begin{eqnarray}&&{W}^{(m)}:= \mathrm{Tr}(\rho ({t}_{m})({H}^{(m)}-{H}^{(m-1)})),\end{eqnarray} \tag{ 23 }$$

where $\rho ({t}_{m})$ is given in (14). The main assumption of this study is precisely that the work cost of a quench is very well approximated by the effective description of the equilibrium state, i.e.

$$\begin{eqnarray}&&{W}^{(m)}=\mathrm{Tr}({\omega }_{(\cdots )}^{(m-1)}({H}^{(m)}-{H}^{(m-1)})),\end{eqnarray} \tag{ 24 }$$

where ${\omega }_{(\cdots )}^{(m-1)}$ is its effective description (16). While we focus our attention on average quantities, primarily for simplicity of the analysis, one could also conceive a study of work extraction under GGE for other work quantifiers [10, 51, 52]. As the equilibration processes happen spontaneously and have no work cost, the total work extracted in the entire protocol is simply given by the sum of the steps

$$\begin{eqnarray}&&W:= \displaystyle \sum _{m=1}^{N}{W}^{(m)}.\end{eqnarray} \tag{ 25 }$$

A particularly relevant scenario is the system-bath setting. We call system S to the part of the total system upon which one has control and it is possible to quench its Hamiltonian H_S. The bath B contains the degrees of freedom upon one has no control and it is the responsible for equilibrating the system S. In order for this equilibration to happen, the dimension of the Hilbert-space of S, $\dim ({{ \mathcal H }}_{S})$, is considered to be much smaller than that of the bath,

$$\begin{eqnarray}&&\dim ({{ \mathcal H }}_{S})\ll \dim ({{ \mathcal H }}_{B})\end{eqnarray} \tag{ 26 }$$

and the total Hamiltonian to be of the form⁴ ,

$$\begin{eqnarray}&&{H}^{(m)}={H}_{S}^{(m)}\otimes {{\mathbb{1}}}_{B}+{{\mathbb{1}}}_{S}\otimes {H}_{B}+V,\end{eqnarray} \tag{ 27 }$$

where the interaction V is supported on S and B and couples the two subsystems. Unlike the standard assumptions in thermodynamics, note that we do not assume that the interaction V is weak or that bath size is infinite. Let us be more explicit about what we mean by that.

Usually, within thermodynamics, it is assumed that the system S equilibrates, upon contact with a bath B, according to

$$\begin{eqnarray}&&{\mathrm{Tr}}_{B}({\omega }_{\beta }^{(m)})={{\rm{\Omega }}}_{\beta }({H}_{S}^{(m)}):= \displaystyle \frac{{{\rm{e}}}^{-\beta {H}_{S}^{(m)}}}{Z},\end{eqnarray} \tag{ 28 }$$

where $\beta \gt 0$ is fixed throughout all the protocol. In contrast, in the model that we consider, given by ${\omega }_{\mathrm{Gibbs}}^{(m)}$ in (16), the inverse temperature changes along the protocol and the Gibbs states describe the whole compound SB. Nonetheless, let us note that the model of equilibration ${{\rm{\Omega }}}_{\beta }$ in (28) represents a particular case of our Gibbsian model ${{\rm{\Omega }}}_{\mathrm{Gibbs}}$ in the limit of weak coupling and infinite bath. In the limit of an infinite bath, the total energy of SB in (16) will not be substantially affected by the energy pumped or subtracted in all the quenches ${H}_{{SB}}^{(m)}\mapsto {H}_{{SB}}^{(m+1)}$ and the parameter ${\beta }^{(m)}$ will remain constant throughout the protocol, ${\beta }^{(m)}\approx \beta$ for all m. In the limit of weak coupling V between S and B, then ${{\rm{\Omega }}}_{{\beta }^{(m)}}({H}_{{SB}}^{(m)})\approx {{\rm{\Omega }}}_{{\beta }^{(m)}}({H}_{S}^{(m)})\otimes {{\rm{\Omega }}}_{{\beta }^{(m)}}({H}_{B}^{(m)})$.

In sum, the model of equilibration ${\omega }_{\mathrm{Gibbs}}^{(m)}$ should be regarded as a correction to the usual setup in thermodynamics given by equation (28). This correction incorporates the fact that the bath is of finite size, which introduces a dependence of the inverse temperature ${\beta }^{(m)}$ and also allows for strong couplings between S and B.

We start by analysing the entropy production of a quasi-static process when all conserved quantities are taken into account. In this case the equilibrium state is given by ${\omega }_{\mathrm{TA}}(u)$. Our first result shows that there is no entropy production in a quasi-static process if the trajectory of Hamiltonians $u\mapsto H(u)$ is smooth.

Result 1 Consider a Hamiltonian trajectory $u\mapsto H(u)={\sum }_{k}{E}_{k}(u)| {E}_{k}(u)\rangle \langle {E}_{k}(u)|$, and a quasi-static process along this trajectory which induces a family of time-average states ${\omega }_{\mathrm{TA}}(u)$. Then, if the trajectory is continuous and the eigenvectors $| {E}_{k}(u)\rangle$ are differentiable, there is no entropy production in such a quasi-static process, that is, $S(0)=S(1)=0$.

The proof and discussion can be found in appendix B. Note that this result is independent of the state which is evolving under H(u). In fact, for a given state, there exist quenches that are not quasi-static but preserve its entropy, such as any quench to a Hamiltonian with the same eigenbasis as the state. This is for instance the case of raising and lowering energy levels.

Now, we consider the case of a generic GGE equilibration where not all the conserved quantities are taken into account. In this case, the equilibration model (16) satisfies the relation,

$$\begin{eqnarray}&&\mathrm{Tr}({\omega }_{\mathrm{GGE}}^{(m)}{Q}_{i}^{(m)})=\mathrm{Tr}({\omega }_{\mathrm{GGE}}^{(m-1)}{Q}_{i}^{(m)}),\end{eqnarray} \tag{ 32 }$$

for all i = 1,...,q. Here the $\{{Q}_{i}^{(m)}\}$ correspond to the q conserved quantities of ${H}^{(m)}$, and equation (32) determines the corresponding Lagrange multipliers ${\lambda }_{i}^{(m)}$ in (4). For such equilibrium states, we also identify conditions so that there is no entropy production. More precisely, we find the following:

Result 2 Consider a quasi-static process along a Hamiltonian trajectory $u\to H(u)$ described by a family of equilibrium states ${\omega }_{\mathrm{GGE}}(u)$ . Then, the entropy of ${\omega }_{\mathrm{GGE}}(u)$ is preserved in such a quasi-static process, provided that the Lagrange-multipliers as determined by (32), form in the limit $N\to \infty$ a set of smooth functions $u\mapsto {\lambda }_{j}(u)$ for j = 1,..., q.

This result is shown simply by taking the continuum limit of equation (32) which yields

$$\begin{eqnarray}&&\mathrm{Tr}\left(\displaystyle \frac{{\rm{d}}{\omega }_{\mathrm{GGE}}(u)}{{\rm{d}}u}{Q}_{j}(u)\right)=0,\quad \forall j=1,\ldots ,m\end{eqnarray} \tag{ 33 }$$

which can be in turn used to show that the entropy production vanishes,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}S}{{\rm{d}}u}=\displaystyle \sum _{j=1}^{m}{\lambda }_{j}(u)\mathrm{Tr}\left(\displaystyle \frac{{\rm{d}}{\omega }_{\mathrm{GGE}}(u)}{{\rm{d}}u}{Q}_{j}(u)\right)=0.\end{eqnarray} \tag{ 34 }$$

Hence, we see that, if the conditions of result 2 are satisfied, the entropy of the effective description in terms of GGE states is also preserved in the limit of a quasi-static process.

Let us now discuss heuristically under which conditions the premise that ${\{u\mapsto {\lambda }_{j}(u)\}}_{1}^{q}$ are smooth functions is expected to be fulfilled. This can be well illustrated by the following example:

Example 1 Consider the case of a two dimensional system for which we take q = 1, that is, the only conserved quantity is the Hamiltonian ${Q}_{1}=H$ itself (the Gibbs equilibration model). Consider initially a non-degenerate Hamiltonian $H(0)=E| 1\rangle \langle 1|$ and an arbitrary initial state $\rho (0)$ with an inverse temperature $\beta (0)\gt 0$ and thus the entropy is smaller than $\mathrm{log}(2)$. Now suppose that the final Hamiltonian $H(1)=0$ has degenerate energy levels. We now show that:

• (i)  
  there is a quasi-static trajectory without a smooth behaviour of the Lagrange-multipliers (in this case $\beta (u)$),
• (ii)  
  this results in a positive entropy production, and
• (iii)  
  how this implies that taking only a single conserved quantity—in this case the energy—does not provide a good approximation of the time-averaged state.

To see the above points, take as Hamiltonian path $H(u)=E(1-u)| 1\rangle \langle 1| =H(0)(1-u)$ and an initial Gibbs state with inverse temperature $\beta (0)$. Then the eigenbasis in the entire process does not change. Now note that the condition (32) implies that the energy is preserved in every equilibration. But since we are dealing with a two-dimensional system, as long as H(u) is non-degenerate, the state itself will remain constant $\omega (u)=\rho (u)$ for any $u\in [0,1)$. This requires that the inverse temperature $\beta (u)\to \infty$ as $u\to 1$: Along the path, the inverse temperature needs to fulfil $\beta (u)=\beta (0)/(1-u)$ to keep the state constant. Therefore, it necessarily diverges as $u\to 1$. To show (ii), simply note that when one reaches H(1), the final state is a maximally mixed state with entropy $\mathrm{log}(2)$, which is larger than the one of the initial state by assumption. To show (iii), observe that the time averaged state would remain constant throughout the protocol, thus it differs from the GGE at u = 1. Similar reasoning as for this example holds true for higher dimensional systems, where the ground state degeneracy of H(1) is higher than that of H(0).

The previous example shows that in some cases the premise of result 2 is not fulfilled, however, these pathological cases often imply that the chosen GGE description is not accurate. For example, in the case of encountering a ground state degeneracy, any conserved quantity in the GGE that discerns the ground states would be enough to fix the problem. However, we leave in general open whether one can find smooth trajectories for $u\mapsto {\lambda }_{j}(u)$ for a given set of conserved quantities and trajectory of Hamiltonians—this may well depend on the specifics of the model and on the ambiguity of what constants of motion to pick in the first place [27].

As discussed above in the case of the GGE ensemble, it is in general necessary to ensure that the Lagrange multipliers $u\mapsto {\lambda }_{j}(u)$ follow a smooth trajectory in order to certify that there is no entropy production in a quasi-static process. This requires to compute the Lagrange multipliers following the model of (16) and keeping track of the conserved quantities. We will see now that the situation simplifies substantially for the case of the Gibbs model of equilibration (where the energy is the only conserved quantity).

Result 3 Consider an initial and final Hamiltonian H(0) and H(1) and initial state ${\omega }_{\mathrm{Gibbs}}(0)={{\rm{e}}}^{-\beta (0)H(0)}/Z$ with finite $\beta (0)\gt 0$. There exist a quasi-static trajectory H(u) so that there is no entropy production if and only if there exist ${\beta }^{* }\gt 0$ so that

$$\begin{eqnarray}&&S({\omega }_{\mathrm{Gibbs}}(0))=S({{\rm{e}}}^{-{\beta }^{* }H(1)}/Z)\end{eqnarray} \tag{ 35 }$$

Note that one of the implications is trivial. The final state is ${{\rm{e}}}^{-{\beta }^{(1)}H(1)}/Z$, hence if there exist no $\beta (1)={\beta }^{* }$ so that (35) is fulfilled, then it is clearly impossible to keep the entropy constant. This can happen if H(1) does not admit any Gibbs state with the initial entropy. The non trivial implication of the previous result is that as long as H(1) admits a Gibbs state with the initial entropy, one can always find a quasi-static trajectory that keeps the entropy constant. Indeed, we find that the quasi-static trajectory achieving it does not need to be fine-tuned. We discuss in appendix C, together with the proof of result 3, that any quasi-static process where the degeneracy of the ground state does not increase along the protocol will indeed keep the entropy constant. This condition is expected to be satisfied for trajectories of generic local Hamiltonians, which have non-degenerate ground spaces for typical choices of the Hamiltonian parameters [53].

We now connect entropy production to reversibility of processes. First, let us note that for the GGE equilibration model (similarly for the Gibbs model since it is a particular case of the former), condition (33) is invariant if one reverses the process. More specifically, given H(u) and ${\omega }_{\mathrm{GGE}}(0)$ as initial state, condition (33) determines the trajectory of states (if the premises of result 2 are met) ${\omega }_{\mathrm{GGE}}(u)$, with u from 0 to 1. Now, we can consider the trajectory $H(\tilde{u})$ with initial state ${\omega }_{\mathrm{GGE}}(\tilde{u}=0)$ with $\tilde{u}=1-u$. One can easily verify that

$$\begin{eqnarray}&&\mathrm{Tr}\left(\displaystyle \frac{{\rm{d}}{\omega }_{\mathrm{GGE}}(\tilde{u})}{{\rm{d}}\tilde{u}}{Q}_{j}(\tilde{u})\right)=0,\quad \forall j=1,\ldots ,m.\end{eqnarray} \tag{ 36 }$$

Hence, the equilibrium state for the trajectory $H(\tilde{u})$ is given exactly by ${\omega }_{\mathrm{GGE}}(\tilde{u}=1-u)$ and thus, the protocol is reversible. In other words, we have seen that for the GGE equilibration model reversible protocols correspond to arbitrarily slow protocols where no entropy is produced on the system and bath together, exactly as is the case for phenomenological thermodynamics.

A well-known feature of phenomenological thermodynamics is that reversible transformations are always beneficial in work-extraction protocols, a phenomenon which is referred to as the minimum work principle. We will later see that this principle naturally holds when the model of equilibration is given by Gibbs states, but its range of applicability is considerably reduced when the equilibrium states are described by GGE. Indeed, we will see explicitly that when the equilibration model is given by a GGE ensemble of free fermions, it can well be beneficial to go through a given protocol quickly and thereby producing entropy.

Before we go on to discuss explicit work extraction protocols, let us stress that the entropy in SB, which can only increase or remain constant, is not simply the sum of the entropies of S and B. This happens because we are considering interacting quantum systems that show correlations between S and B. This is true both in the exact and the effective description. Indeed, in general the von Neumann entropy in SB is smaller than or equal to the sum of local entropies

$$\begin{eqnarray}&&S(\omega )\leqslant S({\mathrm{Tr}}_{B}(\omega ))+S({\mathrm{Tr}}_{S}(\omega )),\end{eqnarray} \tag{ 37 }$$

with equality if and only if $\omega ={\mathrm{Tr}}_{B}(\omega )\otimes {\mathrm{Tr}}_{S}(\omega )$, i.e., when S and B are completely uncorrelated. Thus, entropy-production in our set-up does not always mean that entropy is locally produced in the system and the bath. The generation of entropy is not always associated with the generation of correlations, as in [54], but rather to the mixing induced by equilibration processes. The global entropy may, for example, increase due to a decrease of correlations, but entirely without changing the local states of the system.

As a final remark, note that in the so-called isothermal reversible process the entropy of the system S does not remain constant, while the entropy of the whole compound SB does, as we discuss in examples 1 and 2 in appendix C).

Let us consider the same setup as the one laid out in section 3.1, with an initial state ${\omega }_{\mathrm{Gibbs}}^{(0)}$ and a protocol that performs N quenches according to a certain trajectory $u\mapsto H(u)$, where ${H}^{(m)}:= H(m/N)$. Let us stress that here we do not take the limit of $N\to \infty$ and we keep it general by considering finite N. Let us recall from equation (25) that the total work performed is given by the sum of the individual work ${W}^{(m)}$ in the mth step,

$$\begin{eqnarray}W & = & \displaystyle \sum _{m=1}^{N}{W}^{(m)}=\displaystyle \sum _{m=1}^{N}\mathrm{Tr}({\omega }_{\mathrm{Gibbs}}^{(m-1)}({H}^{(m-1)}-{H}^{(m)}))\\ & = & \mathrm{Tr}({\omega }_{\mathrm{Gibbs}}^{(0)}{H}^{(0)})-\mathrm{Tr}({\omega }_{\mathrm{Gibbs}}^{(N)}{H}^{(N)})\\ & & +\displaystyle \sum _{m=2}^{N}\mathrm{Tr}(({\omega }_{\mathrm{Gibbs}}^{(m-1)}-{\omega }_{\mathrm{Gibbs}}^{(m)}){H}^{(m)}),\end{eqnarray} \tag{ 38 }$$

where in equation (38) we have simply reorganised the terms and added and subtracted the quantity $\mathrm{Tr}({\omega }_{\mathrm{Gibbs}}^{(N)}{H}^{(N)})$. We can now use our model of equilibration, as given by equation (16) that we recall here for completeness,

$$\begin{eqnarray}&&{\omega }_{\mathrm{Gibbs}}^{(m)}={{\rm{\Omega }}}_{\mathrm{Gibbs}}({\omega }_{\mathrm{Gibbs}}^{(m-1)},{H}^{(m)})=\displaystyle \frac{{{\rm{e}}}^{-{\beta }^{(m)}{H}^{(m)}}}{{Z}^{(m)}}\end{eqnarray} \tag{ 39 }$$

for all $m\geqslant 1$, where ${Z}^{(m)}=\mathrm{Tr}({{\rm{e}}}^{-{\beta }^{(m)}{H}^{(m)}})$ and ${\beta }^{(m)}\gt 0$ is determined by the conservation of average energy: $\mathrm{Tr}({\omega }_{\mathrm{Gibbs}}^{(m-1)}{H}^{(m)})=\mathrm{Tr}({\omega }_{\mathrm{Gibbs}}^{(m)}{H}^{(m)})$. One can easily check that energy conservation implies that the last sum in (38) vanishes, which implies that

$$\begin{eqnarray}&&W=\mathrm{Tr}({\omega }_{\mathrm{Gibbs}}^{(0)}{H}^{(0)})-\mathrm{Tr}({\omega }_{\mathrm{Gibbs}}^{(N)}{H}^{(N)}),\end{eqnarray} \tag{ 40 }$$

where ${\omega }_{\mathrm{Gibbs}}^{(N)}$ depends on N and the trajectory H(u).

From equation (40) we see that given a fixed final Hamiltonian H(1), the protocol that costs the minimum amount of work (and maximises the extracted work W) is given by the one that leaves the final state with the least average energy. Since the average energy is monotonic with the entropy for Gibbs states of positive temperature, we conclude that the optimal protocol is the one minimising the entropy of the final state ${\omega }_{\mathrm{Gibbs}}^{(N)}$. Furthermore, as the entropy can only increase throughout the protocol (see section 4), a protocol creating no entropy is optimal.

It has to be stressed that this holds true only as long as the final temperature of the Gibbs state is positive, which happens if

$$\begin{eqnarray}&&\mathrm{Tr}({\omega }_{\mathrm{Gibbs}}^{(N)}{H}^{N})\leqslant \displaystyle \frac{1}{d}\mathrm{Tr}({H}^{(N)}),\end{eqnarray} \tag{ 41 }$$

where d is the dimension of the Hilbert space. Note that the right hand side of the equation typically (e.g., for many body systems with short range interactions) grows linearly with the number of particles. Therefore, if the total system is big enough, we expect condition (41) to be satisfied, and thus the minimum work principle to hold.

Taking together the facts that a protocol creating no entropy is optimal and the results of appendix C—which show conditions so that the quasi-static entropy has no entropy production—one can conclude that the minimal work principle is satisfied for any trajectory so that ${d}_{g}(H(0))\geqslant {d}_{g}(H(u))\geqslant {d}_{g}(H(1))$ (see result 5). As mentioned in section 4.3, this condition is satisfied for trajectories of generic Hamiltonians, which have non-degenerate ground spaces.

Let us now comment on the relation between the minimal work principle and the second law of thermodynamics. First, note that if we fix a trajectory H(u) so that $H(0)=H(1)$, then the final state is a Gibbs state ${{\rm{\Omega }}}_{{\beta }_{N}}(H(0))$. The inverse temperature ${\beta }^{N}$ at the end of the protocol certainly depends on the particular trajectory and the number of quenches performed. However, it is clear by the discussion of section 4 that $S({{\rm{\Omega }}}_{{\beta }_{N}}(H(0)))\geqslant S(\omega (0))$. Hence, since the final state is a Gibbs state with respect to H(0) and with more entropy than the initial Gibbs state and the energy is monotonic with the entropy for Gibbs state, the extracted work is negative. Note that this depends crucially on having Gibbs states as equilibrium states and it will not be reproduced by time-average or GGE models of equilibration as we discuss in the next sections.

Also, the minimum work principle can be used to study work-extraction protocols from non-equilibrium states. As an example, consider as initial conditions a pair of state and Hamiltonian ${\rho }^{(0)}$ and ${H}^{(0)}$ respectively. The goal is to extract work from ${\rho }^{(0)}$ by performing a cyclic protocol, where ${H}^{(N)}={H}^{(0)}$. Note that here the initial state is not in a Gibbs state with respect to the initial Hamiltonian, ${H}^{(0)}$. Nevertheless, after the first quench, it does thermalise to ${\omega }_{\mathrm{Gibbs}}^{(1)}={{\rm{\Omega }}}_{\mathrm{Gibbs}}({\rho }^{(0)},{H}^{(1)})$. From that moment onwards, the minimum work principle can be used, implying that it is always optimal to come back to ${H}^{(0)}$ by a protocol that does not create entropy. The only remaining question is in fact to which Hamiltonian the first quench is performed, an issue that is discussed in appendix D.

We now discuss the minimum work principle for protocols of work extraction when the model of equilibration that is used is the time-average ${{\rm{\Omega }}}_{\mathrm{TA}}$. Let us assume a smooth trajectory of Hamiltonians H(u) and some initial equilibrium state ${\omega }_{\mathrm{TA}}(0)$. Since the trajectory of Hamiltonians is smooth we know that the final state in the quasi-static protocol ${\omega }_{\mathrm{TA}}^{{\rm{q}}.{\rm{s}}.}(1)$ has the same entropy as the initial state ${\omega }_{\mathrm{TA}}(0)$, indeed even the same eigenvalues as ${\omega }_{\mathrm{TA}}(0)$ (see appendix B). The question is, whether this also implies that the quasi-static protocol is optimal in terms of the average work-cost. We will show that this is in general only true if this final state in the quasi-static protocol is also a passive state, meaning that it is diagonal in the energy-eigenbasis and the energy-populations decrease with increasing energy:

$$\begin{eqnarray}&&\mathrm{Tr}(H(1){\omega }_{\mathrm{TA}}^{{\rm{q}}.{\rm{s}}.}(1))=\displaystyle \sum _{k}{({\omega }_{\mathrm{TA}}^{{\rm{q}}.{\rm{s}}}(1))}_{k}^{\downarrow }{E}_{k}(1),\end{eqnarray} \tag{ 42 }$$

where ${({\omega }_{\mathrm{TA}}(1))}^{\downarrow }$ is the vector of eigenvalues of ${\omega }_{\mathrm{TA}}(1)$, ordered such that ${({\omega }_{\mathrm{TA}}(1))}_{k}^{\downarrow }\geqslant {({\omega }_{\mathrm{TA}}(1))}_{l}^{\downarrow }$ if ${E}_{k}(1)\leqslant {E}_{l}(1)$.

Result 4 Given an initial state and a smooth trajectory of Hamiltonians, if the final state in the quasi-static realisation of the protocol is passive, then the the quasi-static realisation of the protocol is optimal.

This result follows, because passive states can only increase their average energy under any unitary transformation [56, 57]:

$$\begin{eqnarray}&&\mathrm{Tr}(H\rho )\leqslant \mathrm{Tr}({HU}\rho {U}^{\dagger }),\quad \rho \,\mathrm{passive}\,{\rm{w}}.{\rm{r}}.{\rm{t}}.\,H.\end{eqnarray} \tag{ 43 }$$

In particular the final state of the quasi-static protocol ${\omega }_{\mathrm{TA}}^{{\rm{q}}.{\rm{s}}.}(1)$ is related to the initial state by some unitary transformation U^* since their spectra are identical. To see that the quasi-static realisation of the protocol is optimal in this case, let us now consider any realisation of the protocol with only a finite number of quenches N and let us denote the final state in a protocol with N quenches as ${\omega }_{\mathrm{TA}}^{N}$. Since the time-average equilibration model can be thought of as applying a mixture of unitaries (evolving the system for some random time) in any finite realisation including N quenches, the final state ${\omega }_{\mathrm{TA}}^{N}$ is related to the initial state by:

$$\begin{eqnarray*}&&{\omega }_{\mathrm{TA}}^{N}=\displaystyle \sum _{i}{p}_{i}{U}_{i}{\omega }_{\mathrm{TA}}(0){U}_{i}^{\dagger }=\displaystyle \sum _{i}{p}_{i}({U}_{i}{U}^{* }){\omega }_{\mathrm{TA}}(1){({U}_{i}{U}^{* })}^{\dagger },\end{eqnarray*}$$

where p_i is some probability distributions of unitaries. But since ${\omega }_{\mathrm{TA}}^{{\rm{q}}.{\rm{s}}.}(1)$ is passive, we henceforth have

$$\begin{eqnarray}&&\mathrm{Tr}(H(1){\omega }_{\mathrm{TA}}^{N})\geqslant \mathrm{Tr}(H(1){\omega }_{\mathrm{TA}}(1)),\end{eqnarray} \tag{ 44 }$$

which proves the claim.

The minimum work principle for cyclic unitary processes was studied in [55] where it was shown that the minimum work-principle holds if: (i) the initial state is passive with respect to the initial Hamiltonian H(0) and (ii) the trajectory of Hamiltonians is such that the initial and final Hamiltonians H(0) and H(1), respectively, do not have a level-crossing w.r.t to each other. Here, by an absence of level-crossing we mean that if ${E}_{i}(0)\geqslant {E}_{j}(0)$, then also ${E}_{i}(1)\geqslant {E}_{j}(1)$ (note that the labelling of the energy-basis is fixed since we require the Hamiltonian trajectory to be smooth). It is now easy to see that under the premise that the initial state is passive, the condition that the final state in the quasi-static realisation is also passive is indeed equivalent to the absence of such a level-crossings. Thus, our result naturally generalises that of [55].

Finally, let us note that given two Hamiltonians H(0) and H(1) and an initial equilibrium state, it is always possible to construct a smooth trajectory of Hamiltonians that connects the two Hamiltonians and such that the final state in the quasi-static protocol is passive and has the same spectrum as it had initially. This can be done with the protocol presented in appendix F. However, note that this protocol requires global control over the Hamiltonians. Once we can only control some part of the Hamiltonian, all the available smooth trajectories might lead to a non-passive final state in the quasi-static realisation, so that it can become beneficial to use a protocol with a finite number of quenches which results in entropy-production.

As in the case of the Gibbs equilibration model, one can easily relate the analysis above to discuss the second law of thermodynamics. First, note that the optimal protocol between H(0) and H(1) is such the final state has the same spectrum and it is passive. Hence, if $H(0)=H(1)$ we conclude that one can extract positive work from the initial equilibrium state ${\omega }_{\mathrm{TA}}(0)$ if and only if it is not passive. Of course this fact is well-known if we consider protocols of work extraction that just apply a unitary transformation to the initial state. Here, we are deriving a similar behaviour with families of protocols that are instead quenches and equilibrations to the time-average state.

In summary, we have identified conditions that ensure that the quasi-static realisation of a given protocol is optimal. This condition generalise the ones found in [55]. Also, we have shown that any state can be brought to its passive form—keeping the same spectrum—by applying a quasi-static protocol over a specific trajectory of Hamiltonians. Altogether, this show that quasi-static protocols are as powerful for work extraction as they can conceivably be.

In this section we briefly analyse notions of work extraction in the case of GGE models of equilibration. Although it is difficult to provide general results for the case of the GGE, without having specified the particular form of the conserved quantities, we do include here a general formulation of the problem at hand as an introduction to particular example of free fermions that we study later. In this situation, the equilibrated states are maximum entropy states

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GGE}}(\rho ,{H}^{(m)},\{{Q}_{j}^{(m)}\}):= {\mathrm{argmax}}_{\sigma \in { \mathcal E }(\rho ,\{{H}^{(m)},{Q}_{j}^{(m)}\})}S(\sigma ).\end{eqnarray} \tag{ 45 }$$

for a collection of constants of motion $\{{Q}_{j}^{(m)}\}$ that are relevant at a given step m of the protocol. For a given protocol, the work extracted is again

$$\begin{eqnarray}&&W=\displaystyle \sum _{m=1}^{N}{W}^{(m)}=\displaystyle \sum _{m=1}^{N}\mathrm{Tr}({\omega }_{\mathrm{GGE}}^{(m-1)}({H}^{(m-1)}-{H}^{(m)})),\end{eqnarray} \tag{ 46 }$$

so that in order to compute the extracted work for a given protocol, one has to keep track of the Lagrange multipliers along that protocol. The optimal work extraction is attained as the supremum of this expression over such protocols. In agreement with our considerations for time-averaged states, here we will find that the minimum work principle is in general not satisfied for many-body models that equilibrate to a GGE. Ultimately, this result is linked to the fact that for GGEs there is in general no direct link between entropy and energy, in strong contrast to the case of Gibbs states. We show this statement by considering specific classes of models for which the GGE is relevant, namely the class of physical systems described by free fermions, a most relevant type of systems that are known to be well described by the GGE. In particular we will show an example where a fast protocol outperforms a slow protocol despite the fact that an effective description by Gibbs states would suggest the opposite.

We consider quadratic fermionic Hamiltonians of the form

$$\begin{eqnarray}&&H=\displaystyle \sum _{i,j=1}^{n}{c}_{i,j}{a}_{i}^{\dagger }{a}_{j},\end{eqnarray} \tag{ 47 }$$

where n is the number of different modes and the fermionic operators satisfy the anti-commutation relations $\{{a}_{i},{a}_{i}^{\dagger }\}={\delta }_{i,j}$, $\{{a}_{i},{a}_{j}\}=\{{a}_{i}^{\dagger },{a}_{j}^{\dagger }\}=0$. The Hamiltonian H can be transformed into

$$\begin{eqnarray}&&H=\displaystyle \sum _{k=1}^{n}{\epsilon }_{k}{\eta }_{k}^{\dagger }{\eta }_{k},\end{eqnarray} \tag{ 48 }$$

where ${\eta }_{k}^{(\dagger )}$ is the annihilation (creation) operator corresponding to the kth eigenmode of the Hamiltonian.

It is well known that equilibrium states of Hamiltonians of the form (48) are not well described by Gibbs states, but rather by GGEs, with the conserved quantities being the occupations of the energy modes ${Q}_{k}={\eta }_{k}^{\dagger }{\eta }_{k}$, k = 1,..., n [27]. Notice that the number of conserved quantities used for the construction of the GGE is the number of distinct modes n and, hence, is linear (and not exponential) in the system size.

We define the correlation matrix $\gamma (\rho )$ of a state ρ as the symmetric matrix having entries

$$\begin{eqnarray}&&{\gamma }_{i,j}(\rho )=\mathrm{Tr}({\eta }_{i}^{\dagger }{\eta }_{j}\rho ).\end{eqnarray} \tag{ 49 }$$

If the state ρ is Gaussian, then $\gamma (\rho )$ contains all information about ρ, and its time evolution under Hamiltonians of the type (48) keeps it Gaussian. In other words, the full density matrix ρ can be reconstructed from just knowing the correlation matrix.

The correlation matrix of the GGE ${{\rm{\Omega }}}_{\mathrm{GGE}}(\rho ,H,\{{\eta }_{k}^{\dagger }{\eta }_{k}\})$ is found by maximising the entropy while preserving all ${Q}_{k}={\eta }_{k}^{\dagger }{\eta }_{k}$, which simply reduces to dephasing the correlation matrix defined in (49) to the diagonal (see appendix G.1 for details). This provides a simple method for obtaining $\gamma ({{\rm{\Omega }}}_{\mathrm{GGE}}(\rho ,H,\{{\eta }_{k}^{\dagger }{\eta }_{k}\}))$.

Note that these GGE descriptions are also Gaussian states. Hence, in the following we can always restrict to Gaussian states. Even when the initial state is not Gaussian, all the results are unchanged if the initial state is replaced by a Gaussian state that has the same correlation matrix. Consequently, in the following, the discussion is reduced to the level of correlation matrices instead of the full density matrices. This allows us to perform numerical simulations of the real time-evolution as well as the effective description of large systems, since they have dimension n × n instead of the ${2}^{n}\times {2}^{n}$ needed to describe the full density matrix.

First we consider optimal protocols for work extraction in a scenario where the Hamiltonian can be transformed to any quadratic Hamiltonian of the form (48). The discussion is similar to that of section 5.2, but in the context of GGE equilibrium states. As in the previous sections, the optimal protocol is the one minimising the final energy, $\mathrm{Tr}({\omega }_{\mathrm{GGE}}^{(N)}{H}^{(0)})={\sum }_{k}{n}_{k}^{(N)}{\epsilon }_{k}^{(0)}$, where we assume the process to be cyclic and ${n}_{k}^{(N)}=\mathrm{Tr}({{\eta }_{k}^{(0)}}^{\dagger }{\eta }_{k}^{(0)}{\omega }_{\mathrm{GGE}}^{(N)})$. In appendix G.2, we show that this minimisation yields $\mathrm{Tr}({\omega }_{\mathrm{GGE}}^{(N)}{H}^{(0)})\geqslant \mathrm{Tr}({\omega }_{\mathrm{GGE}}^{* }{H}^{(0)})$, with

$$\begin{eqnarray}&&\mathrm{Tr}({\omega }_{\mathrm{GGE}}^{* }{H}^{(0)})=\displaystyle \sum _{k=1}^{n}{({d}^{(0)})}_{k}^{\downarrow }{({\epsilon }^{(0)})}_{k}^{\uparrow },\end{eqnarray} \tag{ 50 }$$

where ${d}_{k}^{(0)}$ are the eigenvalues of $\gamma ({\rho }^{(0)})$ and the symbols ${}^{\uparrow }$ and ${}^{\downarrow }$ indicate that the lists are ordered in increasing and decreasing order, respectively. An explicit protocol saturating this bound is constructed in appendix G.2. The optimal protocol is found to be ${\text{}}{reversible}$, so that no entropy is generated, and one needs to perform an arbitrarily large amount of quenches to reach optimality.

In the optimal final state ${\omega }_{\mathrm{GGE}}^{* }$, the diagonal elements of the correlation matrix, corresponding to the population of the energy modes, decay as the energy of the modes increases. This form is reminiscent of the passive states previously introduced. However, in general, states of the form ${\omega }_{\mathrm{GGE}}^{* }$ do not need to be passive: While in passive states the occupation probabilities of the energy eigenstates are decreasing with increasing energy, here only the occupations of the different fermionic modes decrease with the energy of the mode. The total energies are however obtained by combinations of different modes. An example of a state that is non-passive, but where the mode-populations are decreasing with increasing mode-energy is provided in appendix G.2.

Regarding the minimum work principle, one can use a similar line of reasoning as in section 5.2. For a fixed process, the minimum work principle is guaranteed to hold true as long as the possible final states—which are realised by implementing the process at different speeds—have the form (50), i.e., their populations decrease with the energy of the modes. If this condition is not satisfied, the minimum work principle does not hold in general.

In this section we compute the work extracted in different scenarios by (i) a numerical simulation of the exact unitary evolution of the system, (ii) using the effective description in terms of Gibbs states, and (iii) in terms of GGE states. As physical system we consider a chain of fermions, taking as an initial Hamiltonian,

$$\begin{eqnarray}&&{H}^{(0)}=\displaystyle \sum _{i=1}^{n}{\epsilon }_{i}{a}_{i}^{\dagger }{a}_{i}+g\displaystyle \sum _{i=1}^{n-1}({a}_{i}^{\dagger }{a}_{i+1}+{a}_{i+1}^{\dagger }{a}_{i}).\end{eqnarray} \tag{ 51 }$$

First we study the optimal protocol for the case unrestricted Hamiltonian case derived in appendix 5.3, and next we consider the case of local changes of the Hamiltonian. In all cases we find a very good agreement between the real dynamics and the GGE effective description.

Besides comparing the effective descriptions with the real dynamics, we also study the applicability of the minimum work principle. We give an explicit example of a process in which producing entropy is beneficial for work extraction, hence showing that the minimum work principle is violated in this example.

6.3.1. Work extraction with unrestricted Hamiltonians and free fermions

Here, we take as the initial state ${\rho }^{(0)}$ a GGE state whose populations $\gamma {({\rho }^{(0)})}_{i,i}\in (0,1)$ in (49) are chosen i.i.d. from a Gaussian distribution. Again, note that this state is Gaussian. We then apply the protocol described in appendix G.2 for maximal work extraction, and compare the results obtained by the exact dynamics and the GGE model of equilibration. The exact dynamics are computed by, after the the i-th quench, letting the system unitarily evolve under the Hamiltonian ${H}^{(i)}$ for a time much longer than the time scale of equilibration. Figure 2 shows the results obtained using both approaches. It shows a very good agreement, as long as the number of fermions is sufficiently large (in the figure n = 100). Yet small discrepancies are observed, which is due to the fact that we implement global quenches, for which the state may not equilibrate. Note that, when performing local quenches and starting with a Gibbs state, as in figure 3, equilibration of local observables is expected (see section 2) and the agreement is indeed excellent.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can also see in figure 2 how work increases as the process becomes slower, becoming maximal in the limit $N\to \infty$, when reversibility is achieved. This is in agreement with our considerations in section 6.2.

6.3.2. Work extraction with restricted free fermionic Hamiltonians with a Gibbs initial state

Let us now assume that the Hamiltonian can only be locally modified, as discussed in section 3. The Hamiltonian (51) is split in three components. ${H}_{S}={\epsilon }_{1}{a}_{1}^{\dagger }{a}_{1}$ (S is a single fermion), $V=g({a}_{1}^{\dagger }{a}_{2}+{a}_{2}^{\dagger }{a}_{1})$ and ${H}_{B}=H-V-{H}_{S}$. Our capability to change the Hamiltonian is thus reduced to a single parameter: the local energy ${\epsilon }_{1}$. Note that the coupling between the S and the B is not assumed to be weak. The initial state takes the form,

$$\begin{eqnarray}&&{\rho }^{(0)}={\rho }_{S}\otimes {{\rm{\Omega }}}_{\beta }({H}_{B})={\rho }_{S}\otimes \displaystyle \frac{{{\rm{e}}}^{-\beta {H}_{B}}}{Z},\end{eqnarray} \tag{ 52 }$$

where ${\rho }_{S}$ is initially out of thermal equilibrium; for example, in figure 3, it is set to a lower temperature than the bath. As discussed before, we do not need to assume that the initial state ${\rho }_{S}$ (and hence ${\rho }^{(0)}$) is Gaussian, but the work extracted will only depend on its correlation matrix and not on the full density matrix, since the energy is a sum of second moments of the fermionic creation and annihilation operators and all the GGE states constructed in the process are Gaussian automatically.

Figure 3 shows the extracted work from ${\rho }_{S}$ as a function of the number of quenches N, which is computed using the real exact unitary evolution, and the effective description in terms of both GGE and Gibbs states. The agreement between the unitary dynamics and the GGE description is excellent, for any value of N and the parameters, but the Gibbs states fail to describe the process. Even if the bath is initially in a Gibbs state, see equation (52), the posterior evolution of SB can not be correctly described by them. Although the description in terms of Gibbs ensembles is quantitatively incorrect, it is fair to say that it describes some qualitative features of the results. In particular, the exact dynamics satisfies the minimum work principle, and so does the effective description with Gibbs states. This follows because condition (41) is satisfied during the process. However, as we show in the next section, condition (41) can fail to predict the applicability of the minimum work principle.

6.3.3. Work extraction with free fermionic restricted Hamiltonians with a GGE initial state

Equilibrium states when dealing with Hamiltonians of the type (48) are well described by GGE states, it is therefore natural to generalise the initial state (52) to

$$\begin{eqnarray}&&{\rho }_{0}={\rho }_{S}\otimes {\omega }_{\mathrm{GGE}}^{(B)},\end{eqnarray} \tag{ 53 }$$

where ${\omega }_{\mathrm{GGE}}^{(B)}$ is a GGE state with respect to the local Hamiltonian of B, ${H}_{B}={\sum }_{k=1}^{n}{\epsilon }_{k}^{(B)}{\eta }_{k}^{(B)\dagger }{\eta }_{k}^{(B)}.$ Let us now pick a very particular initial state given by

$$\begin{eqnarray}\mathrm{Tr}({\omega }_{\mathrm{GGE}}^{(B)}{\eta }_{k}^{(B)\dagger }{\eta }_{k}^{(B)})=\left\{\begin{array}{ll}1 & k\geqslant K,\\ 0 & k\lt K,\end{array}\right.\end{eqnarray} \tag{ 54 }$$

for some $K\lt n$. That is, only the K most energetic modes are populated. No actual thermal state with positive temperature would have such properties due to the population inversion of the fermionic modes. It is important to acknowledge, however, that if we would chose an effective description as a Gibbs state for such initial states, we would nevertheless obtain a positive effective temperature provided that condition (41) is satisfied. This will be the case as long as the number of populated energy-levels K is small enough. Indeed, for any finite K, but large n, the energy-density in the state is much lower than the critical energy-density needed for negative effective temperatures.

The work extracted in a particular protocol with initial state (53) is plotted in figure 4. The results clearly show how the extracted work decreases with the time spent in the process. Therefore, more work is extracted when more entropy is produced, and the minimum work principle does not apply in this situation. In fact, this is to be expected because both the initial and the final state of the protocol are highly non-passive, and thus the conditions described in section 6.2 are not satisfied. However, when using an effective description in terms Gibbs states, we would have predicted that it is always beneficial to use a quasi-static, reversible protocol since condition (41) is satisfied for the case described in figure 4.

Zoom In Zoom Out Reset image size

6.3.4. Work extraction under imperfect equilibrations

To conclude, we study the accuracy of the model (16) in protocols of work extraction where equilibration is not guaranteed, because the system does not evolve for a sufficiently long time after every quench. In order to study this situation, we consider a state initially in thermal equilibrium with a Hamiltonian ${H}^{(0)}$, which is modified by a sequence of quenches until the Hamiltonian ${H}^{(f)}={H}^{(N)}$ is reached. The main difference with respect to our previous calculations is that we let the total time T of the process be fixed, and take for the equilibration time T/N. Clearly, the accuracy of model (16) is expected to get worse and worse as T decreases. It is less obvious, however, how the accuracy depends on the number of quenches N. On the one hand, as we increase N, the time of equilibration of each quench reduces and hence errors might become larger, and can also accumulate. On the other hand, the distance between subsequents Hamiltonians becomes shorter, so that less time is expected to be needed to reach equilibration.

The results are presented in figure 5. They show that the error, as quantified by $| W-{W}_{\mathrm{GGE}}|$ where W is the extracted work from the real exact dynamics, remains constant for large N, but decreases with the total time T. These results suggest that quasi-static processes can be well described by the model of equilibration (16), in the sense that one does not require an unphysically large time to become arbitrarily close to a quasi-static evolution.

Zoom In Zoom Out Reset image size

As a final remark, we would like to stress that, whereas we have focus our attention in work extraction, the model of equilibration (16) is general and can be applied to other physical quantities—see figure 6 for an example, where we compute the population of one fermion. It is to be expected that the model (16) will give better predictions for local quantities, which are expected to equilibrate in many-body systems.

Zoom In Zoom Out Reset image size

First we consider an idealised scenario where one has full control over the global Hamiltonian H. That is, the Hamiltonians ${H}^{(i)}$ at the ith step of the protocol can be chosen to be any Hamiltonian. Given this maximal level of control, we would like to identify the optimal protocol for work extraction.

We have initial conditions described by a pair of state and Hamiltonian ${\rho }^{(0)}$ and ${H}^{(0)}$ respectively. The goal is to extract work by performing a cyclic protocol, where ${H}^{(N)}={H}^{(0)}$. Importantly, we will no longer assume that the initial state is in a Gibbs state with respect to the initial Hamiltonian ${H}^{(0)}$. After the first quench, the state thermalises to ${\omega }_{\mathrm{Gibbs}}^{(1)}={{\rm{\Omega }}}_{\mathrm{Gibbs}}({\rho }^{(0)},{H}^{(1)})$. Hence, from that moment onwards, the minimum work principle can be applied implying that it is optimal to come back to ${H}^{(0)}$ by a protocol that does not create entropy.

The only remaining question concerning the optimal protocol is to which Hamiltonian ${H}^{(1)}$ the first quench is to be performed. This can be straightforwardly answered by expressing the total work, as in (40),

$$\begin{eqnarray}&&W=\mathrm{Tr}(({\rho }^{(0)}-{\omega }_{\mathrm{Gibbs}}^{(N)}){H}^{(0)}),\end{eqnarray} \tag{ D1 }$$

where by equation (39), we see that ${\omega }_{\mathrm{Gibbs}}^{(N)}$ is a Gibbs state with inverse temperature ${\beta }^{(N)}$. Arguing in the same way as in the minimum work principle, we obtain that the optimal protocol is the one which has no entropy production. Note that a protocol creating zero entropy is only possible for initial states ${\rho }^{(0)}$ such that $S({\rho }^{(0)})=S({{\rm{e}}}^{-{\beta }^{(N)}{H}^{(0)}}/Z)$ for some ${\beta }^{(N)}\gt 0$, as discussed in result 3. Here we provide the steps of a protocol that achieves zero entropy production if that condition is met, which is, as discussed above, the protocol that extracts the maximum amount of work. This protocol reads:

• (1)  
  Apply first a quench from ${H}^{(0)}$ to ${H}^{(1)}=k\mathrm{ln}({\rho }^{(0)})$ for any $k\in {{\mathbb{R}}}^{-}$.
• (2)  
  Let the system equilibrate to ${\omega }_{\mathrm{Gibbs}}^{(1)}:= {{\rm{\Omega }}}_{\mathrm{Gibbs}}({\rho }^{(0)},{H}^{(1)})$ given by (39). The condition of average energy conservation implies that ${\beta }^{(1)}=-1/k$, and thus, ${\omega }_{\mathrm{Gibbs}}^{(1)}={\rho }^{(0)}$.
• (3)  
  Apply a quasi-static process (a sequence of infinitesimal quenches and equilibrations) from ${H}^{(1)}$ to ${H}^{(0)}$. Such process keeps the entropy constant $S({\rho }^{(1)})=S({\rho }^{(N)})$, as discussed in section 4.

This protocol resembles the optimal protocol of work extraction for the model of equilibration of equation (28) [17, 70]; however, the first quench is chosen to a different Hamiltonian ${H}^{(1)}$.

We now consider the restricted case where ${H}^{(i)}\in { \mathcal H }$ and ${ \mathcal H }$ is a given set of Hamiltonians. While we will later be interested in the case where restriction are such that we can only change the initial Hamiltonian locally on the subsystem S, so that

$$\begin{eqnarray}&&{{ \mathcal H }}_{\mathrm{Local}}=\{H\ \ | \ \ H={H}^{(0)}+{H}_{S}\otimes {{\mathbb{1}}}_{B}\},\end{eqnarray} \tag{ D2 }$$

we will keep the discussion completely general.

In the same way as in the case of unrestricted Hamiltonians, a maximum amount of work will be extracted if we minimise the final energy, as expressed by equation (40). Since the final state is by assumption a Gibbs state, it is therefore optimal to end up with a Gibbs state with minimal possible positive temperature (every state with negative temperature has higher energy than all states with positive temperature). This is clearly possible if the initial state already has an effective positive temperature with respect to any Hamiltonian in ${ \mathcal H }$. We will assume from now on that this is the case.

Considering steps (1)–(3) of protocol in section D.1, one can easily see that step (1) cannot be applied if $k\mathrm{ln}({\rho }^{(0)})\notin { \mathcal H }$. Instead, we will make a quench ${H}^{(0)}\mapsto {H}^{(1)}$ with

$$\begin{eqnarray}&&{H}^{(1)}={\mathrm{argmin}}_{H\in { \mathcal H }}\ \ S({{\rm{\Omega }}}_{\mathrm{Gibbs}({\rho }^{(0)},H)}),\end{eqnarray} \tag{ D3 }$$

while steps  (2)–(3) are not modified by the restrictions on ${ \mathcal H }$.

We consider Hamiltonians of the type

$$\begin{eqnarray}&&H=\displaystyle \sum _{i,j}{c}_{i,j}{a}_{i}^{\dagger }{a}_{j},\end{eqnarray} \tag{ G1 }$$

where the operators ${a}_{i},{a}_{i}^{\dagger }$ satisfy the fermionic anti-commutation relations,

$$\begin{eqnarray}&&\{{a}_{i},{a}_{j}^{\dagger }\}={\delta }_{i,j},\end{eqnarray} \tag{ G2 }$$

$$\begin{eqnarray}&&\{{a}_{i},{a}_{j}\}=\{{a}_{i}^{\dagger },{a}_{j}^{\dagger }\}=0.\end{eqnarray} \tag{ G3 }$$

Since the matrix c in (G1) is Hermitian, it can be diagonalised by a unitary operator, $c={{ADA}}^{\dagger }$, where ${{AA}}^{\dagger }=1$ and $D=\mathrm{diag}\{{\epsilon }_{1},\ldots ,{\epsilon }_{n}\}$. The Hamiltonian then can be expressed as,

$$\begin{eqnarray}&&H=\displaystyle \sum _{k}{\epsilon }_{k}{\eta }_{k}^{\dagger }{\eta }_{k},\end{eqnarray} \tag{ G4 }$$

with

$$\begin{eqnarray}&&{\eta }_{k}=\displaystyle \sum _{j}{A}_{j,k}^{* }{a}_{j},\end{eqnarray} \tag{ G5 }$$

$$\begin{eqnarray}&&{\eta }_{k}^{\dagger }=\displaystyle \sum _{j}{A}_{j,k}{a}_{j}^{\dagger }.\end{eqnarray} \tag{ G6 }$$

The fermionic operators ${\eta }_{k}^{\dagger },{\eta }_{k}$ are usually referred to as normal modes. The unitarity of A ensures that the transformation preserves the commutation relations,

$$\begin{eqnarray}&&\{{\eta }_{k},{\eta }_{l}^{\dagger }\}=\displaystyle \sum _{i,j}{A}_{k,i}{A}_{l,j}^{* }\{{a}_{i},{a}_{j}^{\dagger }\}={\delta }_{k,l}.\end{eqnarray} \tag{ G7 }$$

where we used (G3).

In the following, we will describe states within the framework of correlation matrices. Define the entries of the correlation matrix $\gamma (\rho )$ corresponding to ρ as

$$\begin{eqnarray}&&{\gamma }_{a}{(\rho )}_{i,j}=\mathrm{Tr}({a}_{i}^{\dagger }{a}_{j}\rho ).\end{eqnarray} \tag{ G8 }$$

Notice that the diagonal elements represent the occupation probabilities, or populations, of each physical fermion. The correlation matrix in the diagonal basis ${\gamma }_{\eta }{(\rho )}_{i,j}=\mathrm{Tr}({\eta }_{i}^{\dagger }{\eta }_{j}\rho )$ is related to ${\gamma }_{a}$ through ${\gamma }_{\eta }={A}^{T}{\gamma }_{a}{A}^{* }$. The diagonal elements of ${\gamma }_{\eta }$, corresponding to the populations of the normal modes, play an important role, and we denote them by p_k,

$$\begin{eqnarray}&&{p}_{k}=\mathrm{Tr}({\eta }_{k}^{\dagger }{\eta }_{k}\rho ).\end{eqnarray} \tag{ G9 }$$

The time evolution of $\gamma (\rho )$ under H, $\rho (t)={{\rm{e}}}^{-{\rm{i}}{Ht}}\rho {{\rm{e}}}^{{\rm{i}}{Ht}}$, can be easily computed in the Heisenberg picture,

$$\begin{eqnarray}&&\dot{{\eta }_{k}}={\rm{i}}[H,{\eta }_{k}]=-{\rm{i}}{E}_{k}{\eta }_{k},\end{eqnarray} \tag{ G10 }$$

$$\begin{eqnarray}&&{\eta }_{k}(t)={{\rm{e}}}^{-{\rm{i}}{E}_{k}t}{\eta }_{k},\end{eqnarray} \tag{ G11 }$$

where we have used $\{{\eta }_{i},{\eta }_{j}^{\dagger }\}={\delta }_{i,j}$ and ${\eta }_{k}^{2}=0$. Therefore, on the one hand, it follows that

$$\begin{eqnarray}&&{\gamma }_{\eta }(\rho (t))={{\rm{e}}}^{{\rm{i}}{tD}}{\gamma }_{\eta }(\rho ){{\rm{e}}}^{-{\rm{i}}{tD}}\end{eqnarray} \tag{ G12 }$$

with $D=\mathrm{diag}\{{E}_{1},\ldots ,{E}_{n}\}$. In the original basis it reads,

$$\begin{eqnarray}&&{\gamma }_{a}(\rho (t))=U{\gamma }_{a}(\rho ){U}^{\dagger }\end{eqnarray} \tag{ G13 }$$

with $U={A}^{* }{{\rm{e}}}^{{\rm{i}}{tD}}{A}^{T}$. On the other, the time averaged state, which is defined as,

$$\begin{eqnarray}&&\langle \rho {\rangle }_{t}=\mathop{\mathrm{lim}}\limits_{T\to \infty }\displaystyle \frac{1}{T}{\int }_{0}^{T}\rho (t),\end{eqnarray} \tag{ G14 }$$

is represented simply by

$$\begin{eqnarray}&&{\gamma }_{\eta }(\langle \rho {\rangle }_{t})=\langle {\gamma }_{\eta }(\rho ){\rangle }_{t}={\rm{\Gamma }}[{\gamma }_{\eta }(\rho (t))],\end{eqnarray} \tag{ G15 }$$

where Γ corresponds to a de-phasing operation. In fact, this correlation matrix is the same as the one of the GGE where the conserved quantities are the normal modes ${\eta }_{k}^{\dagger }{\eta }_{k}$, i.e.,

$$\begin{eqnarray}&&\gamma ({{\rm{\Omega }}}_{\mathrm{GGE}}(\rho ,H,\{{\eta }_{k}^{\dagger }{\eta }_{k}\}))=\gamma ({\langle \rho \rangle }_{t}).\end{eqnarray} \tag{ G16 }$$

Note however, that this does not imply that the full quantum state of the GGE is the same as the time-averaged state.

A particularly important class of fermionic states is given by Gaussian states. In this situation of fixed particle number, such Gaussian states are completely determined from the correlation matrix. In particular, eigenstates and thermal states of free fermionic Hamiltonians are Gaussian states, but clearly also the GGEs given above, as they are obtained by maximizing the entropy given the expectation values of the operators ${\eta }_{k}^{\dagger }{\eta }_{k}$.

If a state ρ is Gaussian, the entropy of ρ can be calculated as

$$\begin{eqnarray}&&S(\rho )=\displaystyle \sum _{k}{\rm{H}}({d}_{k}),\end{eqnarray} \tag{ G17 }$$

where d_k are the eigenvalues of $\gamma (\rho )$, and ${\rm{H}}(p)=-p\mathrm{ln}p-(1-p)\mathrm{ln}(1-p)$. This fact allows us to study entropy-production numerically for large systems.

Here we find a bound for work extraction protocols, which, as discussed in the main text, is equivalent to finding a lower bound on the final energy, $\mathrm{Tr}({\omega }_{\mathrm{GGE}}^{(N)}{H}^{(0)})$, with ${H}^{(0)}={\sum }_{k}{\epsilon }_{k}^{(0)}{\eta }_{k}^{\dagger }{\eta }_{k}^{(0)}$. From our considerations in section G.1, it follows that under the joint operation of a quench,

$$\begin{eqnarray}&&{H}^{(i)}=\displaystyle \sum _{i,j}{c}_{i,j}^{(i)}{a}_{i}^{\dagger }{a}_{j}\mapsto {H}^{(i+1)}=\displaystyle \sum _{i,j}{c}_{i,j}^{(i+1)}{a}_{i}^{\dagger }{a}_{j},\end{eqnarray} \tag{ G18 }$$

followed by an equilibration process, our effective description in terms of GGE states takes the form

$$\begin{eqnarray}&&{\gamma }_{a}({\omega }_{\mathrm{GGE}}^{(i+1)})={A}_{(i+1)}^{* }{\rm{\Gamma }}[{A}_{(i+1)}^{T}{\gamma }_{a}({\omega }_{\mathrm{GGE}}^{(i)}){A}_{(i+1)}^{* }]{A}_{(i+1)}^{T},\end{eqnarray} \tag{ G19 }$$

where Γ is a de-phasing operation, and ${c}^{(i+1)}={A}_{i+1}{{DA}}_{i+1}^{\dagger }$, with D a diagonal matrix. Let $\{{d}_{k}^{(i+1)}\}$ and $\{{d}_{k}^{(i)}\}$ be the eigenvalues of ${\gamma }_{a}({\omega }_{\mathrm{GGE}}^{(i+1)})$ and ${\gamma }_{a}({\omega }_{\mathrm{GGE}}^{(i)})$, respectively. Under (G19), they are related through a doubly stochastic matrix,

$$\begin{eqnarray}&&{d}_{k}^{(i+1)}=\sum {C}_{k,l}{d}_{l}^{(i)}\end{eqnarray} \tag{ G20 }$$

with ${\sum }_{k}{C}_{k,l}={\sum }_{l}{C}_{k,l}=1$. Therefore, the eigenvalues of the final state ${\gamma }_{a}({\omega }_{\mathrm{GGE}}^{(N)})$ can also be expressed as a stochastic combination of the eigenvalues of ${\gamma }_{a}({\rho }^{(0)})$, $\{{d}_{k}^{(0)}\}$. It now follows from basic notions of the theory of majorisation that,

$$\begin{eqnarray}&&\mathrm{Tr}({\omega }_{\mathrm{GGE}}^{(N)}{H}^{(0)})\geqslant \displaystyle \sum _{k=1}^{n}{({d}_{k}^{(0)})}^{\downarrow }{({\epsilon }_{k}^{(0)})}^{\uparrow }=\mathrm{Tr}({\omega }_{\mathrm{GGE}}^{* }{H}^{(0)})\end{eqnarray} \tag{ G21 }$$

where ${}^{\uparrow }$ and ${}^{\downarrow }$ reflect lists ordered in increasing (decreasing) order. This provides the bound (50).

We now construct an explicit protocol that achieves this bound in the limit $N\to \infty$, where N is the number of quenches performed. Let ${\gamma }_{a}({\rho }^{(0)})$ be the correlation matrix of ${\rho }^{(0)}$ as in (G8). First, find some U that diagonalises $\gamma ({\rho }^{(0)})$,

$$\begin{eqnarray}&&U\gamma ({\rho }^{(0)}){U}^{\dagger }=D,\end{eqnarray} \tag{ G22 }$$

and make a quench to

$$\begin{eqnarray}&&{H}^{(1)}={U}^{T}{H}^{(0)}{U}^{* }.\end{eqnarray} \tag{ G23 }$$

Since $\gamma ({\rho }^{(0)})$ is diagonal in the new basis, it follows that ${\omega }_{\mathrm{GGE}}^{(1)}={\rho }^{(0)}$, i.e., the state is not changed during the equilibration process. Now, slowly rotate back to the original Hamiltonian, by performing $N/2$ quenches (followed by equilibration processes) until ${H}^{(0)}$ is reached. At the end the state, ${\rho }^{(N/2)}$ is (approximately) diagonal with respect to the original Hamiltonian, ${H}^{(0)}$. Next, find some V that order the populations of $\gamma ({\rho }^{(N/2)})$, so that $V\gamma ({\rho }^{(N/2)}){V}^{\dagger }$ satisfies (50). As before, perform a quench to

$$\begin{eqnarray}&&{H}^{(1)}={V}^{T}{H}^{(0)}{V}^{* },\end{eqnarray} \tag{ G24 }$$

and slowly come back to the original Hamiltonian by performing $N/2$ quenches. This process give rise to the desired final state ${\omega }_{\mathrm{GGE}}^{* }$ in the limit of infinitesimally slow transformations, i.e., in the limit $N\to \infty$. The optimal protocol is therefore reversible, and it agrees with our intuition that slow processes are better for work extraction.

Importantly, note that these results for the free fermions are completely analogue to the case of time average equilibrium state, as detailed in section 5.1. Indeed, the optimal final state resembles a passive state, which is the optimal final state for work extracting protocols using time-averaged states. However, it should be stressed that the GGE equilibration model considered for free fermions does not coincide in general with the time averaged state. Indeed, this difference can be highlighted by looking at the final state obtained for the time average model in comparison with the final state of the GGE equilibration for free fermions. In the former, one ends up with a passive state. This implies, for n fermions, ${2}^{n}$ energy populations decrease with the energy. On the other hand, for the GGE model of equilibration considered here, the final state ${\omega }_{\mathrm{GGE}}^{* }$ is such only the n populations of the normal modes need to be in decreasing order. These two states are in general not the same.

For example, consider a three-fermion system with Hamiltonian

$$\begin{eqnarray}&&H={\epsilon }_{1}{\eta }_{1}^{\dagger }{\eta }_{1}+{\epsilon }_{2}{\eta }_{2}^{\dagger }{\eta }_{2}+{\epsilon }_{3}{\eta }_{3}^{\dagger }{\eta }_{3}\end{eqnarray} \tag{ G25 }$$

and a state ρ with $\mathrm{Tr}({\eta }_{i}^{\dagger }{\eta }_{i}\rho )={p}_{i}$ with $i=1,2,3$. The quantum state ρ and H can be written as

$$\begin{eqnarray}H & = & \mathrm{diag}\{0,{\epsilon }_{1},{\epsilon }_{2},{\epsilon }_{3},{\epsilon }_{1}+{\epsilon }_{2},{\epsilon }_{2}+{\epsilon }_{3},{\epsilon }_{1}+{\epsilon }_{3},{\epsilon }_{1}+{\epsilon }_{2}+{\epsilon }_{3}\}\\ \rho & = & \mathrm{diag}\{(1-{p}_{1})(1-{p}_{2})(1-{p}_{3}),{p}_{1}(1-{p}_{2})(1-{p}_{3}),\\ & & {p}_{2}(1-{p}_{1})(1-{p}_{3}),{p}_{3}(1-{p}_{1})(1-{p}_{2}),{p}_{2}{p}_{1}(1-{p}_{3}),\\ & & {p}_{2}{p}_{3}(1-{p}_{1}),{p}_{1}{p}_{3}(1-{p}_{2}),{p}_{1}{p}_{2}{p}_{3}\}.\end{eqnarray} \tag{ G26 }$$

If we now choose ${\epsilon }_{1}=1$, ${\epsilon }_{2}=2$, ${\epsilon }_{3}=2.5;$ and p₁ = 0.4, p₂ = 0.3, and ${p}_{3}=0.1;$ we obtain that ρ is not passive but has the form of ${\omega }_{\mathrm{GGE}}^{* }$.

The origin of the difference is the set of operations in which every state is defined. Passive states arise as optimal states for work extraction protocols if any unitary operation can be performed to the system, or, equivalently, every cyclic process in which the system remains thermally isolated. On the other hand, states ${\omega }_{\mathrm{GGE}}^{* }$ become optimal when the set of operations corresponds to (arbitrary) quenches to quadratic Hamiltonians, which is in general more constraint that the set of unitary operations. Within this constrained set of operations, they become optimal.