Work extraction in an isolated quantum lattice system: Grand canonical and generalized Gibbs ensemble predictions

TL;DR

This paper investigates work extraction in isolated quantum lattice systems using quenches and equilibration to GGE and GE, identifying conditions for maximal work extraction without entropy increase.

Contribution

It introduces a class of quenches that preserve entropy and maximize work extraction in both noninteracting and weakly interacting quantum systems.

Findings

Quenches that do not produce entropy enable maximal work extraction.

Equilibration to GGE or GE can be achieved without entropy increase under certain quenches.

Maximal work extraction is possible in both integrable and weakly interacting systems.

Abstract

We study work extraction (defined as the difference between the initial and the final energy) in noninteracting and (effectively) weakly interacting isolated fermionic quantum lattice systems in one dimension, which undergo a sequence of quenches and equilibration. The systems are divided in two parts, which we identify as the subsystem of interest and the bath. We extract work by quenching the on-site potentials in the subsystem, letting the entire system equilibrate, and returning to the initial parameters in the subsystem using a quasi-static process (the bath is never acted upon). We select initial states that are direct products of thermal states of the subsystem and the bath, and consider equilibration to the generalized Gibbs ensemble (GGE, noninteracting case) and to the Gibbs ensemble (GE, weakly interacting case). We identify the class of quenches that, in the thermodynamic limit, results in GE and GGE entropies after the quench that are identical to the one in the initial state (quenches that do not produce entropy). Those quenches guarantee maximal work extraction when thermalization occurs. We show that the same remains true in the presence of integrable dynamics that results in equilibration to the GGE.

Equations

$${\hat{\mathcal{H}}}_s(V_s)$$

$${\hat{\mathcal{H}}}_b(V_b)$$

$$\hat{\mathcal{H}}(V_s,V_b)={\hat{\mathcal{H}}}_s(V_s)+{\hat{\mathcal{H}}}_b(V_b)-({\hat{c}}_{L_s}^{†}{\hat{c}}_{L_s+1}+\text{H.c.}).$$

$$W=\text{mbox}Tr\left[({\hat{\rho }}^I-{\hat{\rho }}^F)[{\hat{\mathcal{H}}}_s(V_s^I)+{\hat{\mathcal{H}}}_b(V_b)]\right],$$

$${\hat{\rho }}^I={\hat{\rho }}_s^I\otimes {\hat{\rho }}_b^I.$$

$${\hat{\rho }}_s^I$$

$${\hat{\rho }}_b^I$$

$$Z^{I}_{s(b)}=\prod_{i=1}^{L_{s(b)}}\left[1+e^{-\beta^{I}_{s(b)}\big{(}\epsilon_{s(b)}^{i}-\mu^{I}_{s(b)}\big{)}}\right],$$

$$E_{s(b)}^I=\text{mbox}Tr[{\hat{\rho }}_{s(b)}^I{\hat{\mathcal{H}}}_{s(b)}^I]\text{and}{\mathcal{N}}_{s(b)}^I=\text{mbox}Tr[{\hat{\rho }}_{s(b)}^I{\hat{\mathcal{N}}}_{s(b)}],$$

$${\hat{\mathcal{H}}}^I\equiv {\hat{\mathcal{H}}}_s(0)+{\hat{\mathcal{H}}}_b(0),$$

$$S^I=S_s^I+S_b^I,$$

$$S_{s(b)}^I=\ln Z_{s(b)}^I+{\beta }_{s(b)}^I[E_{s(b)}^I-{\mu }_{s(b)}^I{\mathcal{N}}_{s(b)}^I].$$

$$E^I=E_s^I+E_b^I\text{and}{\mathcal{N}}^I\equiv {\mathcal{N}}_s^I+{\mathcal{N}}_b^I,$$

$$E_Q=E^I+V{n}_s^I{L}_s\text{and}{E}_Q^{\prime }=E^I+V{n}_I{L}_s,$$

$$W_{\text{max}}^{\text{GE}}=\text{mbox}Tr\left[({\hat{\rho }}^I-{\hat{\rho }}_{\text{max}}^{\text{GE}})[{\hat{\mathcal{H}}}_s(0)+{\hat{\mathcal{H}}}_b(0)]\right],$$

$$n_s^I$$

$$n_b^I$$

$$E_s^I$$

$$E_b^I$$

$${\hat{\rho }}_Q^{\text{GE}}=\exp (-{\beta }_Q[\hat{\mathcal{H}}(V,0)-{\mu }_Q\hat{\mathcal{N}}])/Z_Q,$$

$$E_s^{\text{GE}}$$

$$E_b^{\text{GE}}$$

$$E_s^{\text{GE}}$$

$$E_s^{\text{GE}}+E_b^{\text{GE}}=E_s^I+E_b^I+V{n}_s^I{L}_s$$

$$S^I$$

$$S_Q^{\text{GE}}$$

$${\rho }_{ij}^{\text{GE}}={\delta }_{ij}-{\left[I+e^{-\beta (H-\mu )}\right]}_{ji}^{-1},$$

$$E_{F,1}^{\text{GE}}$$

$$E_{F,2}^{\text{GE}}$$

$$E_{F,N}^{\text{GE}}$$

$$n_{s,q_m}^{\text{GE}}\simeq n_{s,Q}^{\text{GE}}+\frac{(n_{s,F}^{\text{GE}}-n_{s,Q}^{\text{GE}})m}{N},$$

$$W^{\text{GE}}(N)=E^I-E_{F,N}^{\text{GE}}$$

$$\simeq V{L}_s\left(\frac{n_{s,F}^{\text{GE}}+n_{s,Q}^{\text{GE}}}{2}-n_s^I\right)-\frac{V{L}_s}{2N}\left(n_{s,F}^{\text{GE}}-n_{s,Q}^{\text{GE}}\right).$$

$$S^{\text{GE}}=-\sum _{\alpha =1}^L[I_{\alpha }^{\text{GE}}\ln I_{\alpha }^{\text{GE}}+(1-I_{\alpha }^{\text{GE}})\ln (1-I_{\alpha }^{\text{GE}})],$$

Full text

I Introduction

Work is a familiar concept in the context of classical thermodynamics, which deals with systems with a large number of particles. A goal of classical thermodynamics is to identify protocols that provide an efficient way of converting heat into work. In recent years, there has been a lot of interest in developing a thermodynamic framework to deal with small systems that can be far from equilibrium jarzynski2011equalities ; seifert2012stochastic . In particular, following seminal work by Jarzynski jarzynski1997nonequilibrium , fluctuation theorems have become an area of intense activity. Fluctuation theorems have been generalized to understand the stochastic fluctuations of work done on non-equilibrium quantum systems crooks1999entropy ; kurchan2000quantum ; tasaki2000jarzynski ; campisi2011colloquium . Another recent development comes from the so-called single-shot information theory, which was initially postulated to study finite-size effects in quantum cryptography. Single-shot information theory has become a useful tool in understanding work extraction in the context of quantum thermodynamics aaberg2013truly ; horodecki2013fundamental .

In parallel, extraordinary advances in experiments with ultracold atomic gases cazalilla_11 ; bloch_dalibard_08 have motivated much research on the far-from-equilibrium dynamics and the description after equilibration of isolated many-body quantum systems d2016quantum ; eisert_friesdorf_15 ; polkovnikov_sengupta_11 . Of particular interest has been the dynamics following the so-called quantum quenches calabrese_cardy_06 , in which the system is initially in a stationary state of some time-independent Hamiltonian and that Hamiltonian is suddenly changed into a new one that is also time-independent. If the Hamiltonian after the quench is quantum chaotic, i.e., if its distribution of many-body energy level spacings is of the Wigner-Dyson type, one expects thermalization to occur santos2010onset ; d2016quantum . Namely, one expects that after equilibration observables are described by traditional statistical mechanics d2016quantum ; rigol2008thermalization ; rigol_14 . This can be understood to be a consequence of eigenstate thermalization deutsch1991quantum ; srednicki1994chaos ; rigol2008thermalization ; rigol2012alternatives . On the other hand, if the Hamiltonian after the quench is integrable, one expects generalized thermalization to occur. Namely, one expects that after equilibration observables are described by a generalized Gibbs ensemble (GGE), which takes into account the presence of an extensive set of nontrivial conserved quantities rigol2007relaxation ; calabrese_essler_11 ; ilievski_denardis_15 (see Refs. vidmar_rigol_16 ; essler_fagotti_16 ; cazalilla_chung_16 ; caux_16 for recent reviews). This can be understood to be a consequence of generalized eigenstate thermalization cassidy2011generalized ; vidmar_rigol_16 ; caux_essler_13 .

Here, we explore work extraction in isolated (integrable) noninteracting and (effectively) weakly interacting fermionic quantum lattice systems in one dimension as described by a quadratic Hamiltonian. The isolated systems are divided in two parts, which we identify as the subsystem of interest and the bath. The specific protocol we consider is motivated by a possible straightforward implementation in experiments, and consists of: (i) quenches of on-site potentials in the subsystem; (ii) equilibration to the GGE (noninteracting case) or the grand canonical ensemble (GE, weakly interacting case); and (iii) return to the initial values of the on-site potentials in the subsystem by means of quasi-static process with equilibration to the GGE and the GE (the bath is never acted upon). The work extracted is computed as the difference between the initial and the final energy of the entire system. Since the average energy of a thermal state (with non-negative temperature) can only increase due to unitary operations (quantum quenches), because of passivity pusz1978passive ; lenard1978thermodynamical , the initial states cannot be thermal equilibrium states of the entire system (we would like to be able to extract work). Instead, they are selected to be direct products of grand canonical states of the subsystem and the bath at the same temperature but at different chemical potentials (different site occupations).

In a recent study, Perarnau-Llobet et al. perarnau2016work discussed upper bounds for the work that can be extracted in processes involving equilibration to the GE or the GGE in isolated quantum systems. In the context of our protocol, we identify a class of quenches that do not produce entropy when equilibration occurs to the GE or the GGE, which automatically ensures maximal work extraction for equilibration to the GE. We show that those quenches saturate the bound for work extraction under equilibration to the GGE.

The paper is organized as follows. In Sec. II, we introduce the model, quench protocol, and the initial states considered. We discuss the results for equilibration to the GE in Sec. III, and for equilibration to the GGE in Sec. IV. In Sec. V, we present a comparison between the results obtained for the GE and for the GGE. A summary of our results is presented in Sec. VI.

II Model, quench protocol, and initial states

We study noninteracting (and, effectively, weakly interacting) fermions in one-dimensional (1D) lattices with open boundary conditions. The system is divided in two parts, which we identify as the subsystem of interest and the bath. They are described by the tight-binding quadratic Hamiltonians $\hat{\mathcal{H}}_{s}$ and $\hat{\mathcal{H}}_{b}$, respectively

[TABLE]

where $\hat{c}^{{\dagger}}_{i}$ ($\hat{c}_{i}$) is the fermionic creation (annihilation) operator at site $i$, $\hat{n}_{i}=\hat{c}^{{\dagger}}_{i}\hat{c}_{i}$ is the number operator, and $L_{s}$ ($L$) is the size of the subsystem (entire system). We have set the hopping amplitudes in the subsystem and bath to unity, and $V_{s}$ ($V_{b}$) is the on-site potential of the subsystem (bath). Dynamics are studied under the total Hamiltonian

[TABLE]

We prepare the initial state to be a direct product of GE density matrices of the subsystem and bath with $\hat{\mathcal{H}}_{s}(V^{I}_{s})$ and $\hat{\mathcal{H}}_{b}(V_{b})$, respectively (see Sec. II.1). [This can be done by (weakly) connecting the subsystem and the bath to two reservoirs (see Fig. 1).] Next, we connect the subsystem and the bath, and at the same time quench the on-site potentials of the subsystem from $V^{I}_{s}$ to $V^{F}_{s}$ (see Fig. 1). The entire system is then allowed to equilibrate to the GE and the GGE under $\hat{\mathcal{H}}(V^{F}_{s},V_{b})$. After equilibration, which is ensured by taking the density matrix of the entire system to be the appropriate GE or GGE density matrix, we apply a large number $N$ of weak quenches followed by equilibration to the GE or the GGE. In each of those $N$ weak quenches, the on-site potential in the subsystem is changed by $(V^{I}_{s}-V^{F}_{s})/N$, so that at the end we have the subsystem at the initial value $V^{I}_{s}$. In a final ($N+1$) quench, we turn off the hopping between the subsystem and the bath and let the system equilibrate to the GE or the GGE. Note that the latter is a local quench, i.e., it does not produce extensive changes in thermodynamic quantities. This completes our cyclic process (see Fig. 1). (One can prepare again the initial direct product of GE density matrices by connecting the subsystem and the bath to the two reservoirs.)

Relaxation to the GGE is assumed in order to describe what happens under the dynamics dictated by the quadratic Hamiltonian $\hat{\mathcal{H}}(V_{s},V_{b})$. It is actually straightforward to prove that the infinite-time average of the entire one-body density matrix of a noninteracting system (from which all observables in our noninteracting system can be computed) is, in the absence of degeneracies in the single-particle spectrum, identical to that of the GGE he2013single ; ziraldo_santoro_13 . However, the density matrix of the entire system, at any time after a quench, is never that of the GGE (the dynamics is unitary). This is also true for the one-body density matrix, because the fermions are noninteracting wright_rigol_14 . Hence, one may be wary about replacing the density matrix “after equilibration” by that of the GGE. We report numerical results that support the appropriateness of this procedure (see also Ref. perarnau2016work ), as the exact time evolution in which one waits random times after equilibration and the GGE replacement produce nearly identical results for the work extraction.

Relaxation to the GE is assumed in order to describe what happens under the dynamics dictated by the quadratic Hamiltonian $\hat{\mathcal{H}}(V_{s},V_{b})$ plus very weak integrability-breaking interactions. We imagine the weak integrability-breaking interactions as allowing the systems to thermalize at long times rigol_14 ; bertini_essler_15 ; brandino_caux_15 ; rigol2016fundamental ; bertini_essler_16 , while being weak enough not to significantly change the thermal expectation value of macroscopic observables (such as the energy, which is needed to compute the work extracted) from the result in the noninteracting limit. The fact that one can replace the density matrix of the time-evolving state of a quantum chaotic system after equilibration, without needing to wait random times as in the integrable case, by that of the GE has been discussed in Refs. ji_fine_2011 ; ikeda_sukumichi_15 ; d2016quantum .

The work extracted due to our cyclic process, $W$, is defined as

[TABLE]

which is the difference between the energy in the initial and final states alicki2013entanglement ; d2016quantum , where $\hat{\rho}^{I}$ ($\hat{\rho}^{F}$) is the density matrix of the initial (final) state. For all calculations reported here, we take $V^{I}_{s}=V_{b}=0$ and $V^{F}_{s}=V>0$.

II.1 Initial States

We consider initial states that are product states of the subsystem and the bath, i.e., whose density matrix can be written as

[TABLE]

We take the density matrices of the subsystem and the bath to be grand canonical,

[TABLE]

respectively. $\hat{\mathcal{N}}_{s}$ ($\hat{\mathcal{N}}_{b}$) is the total number of particles operator of the subsystem (bath), $\beta^{I}_{s}$ ($\beta^{I}_{b}$) and $\mu^{I}_{s}$ ($\mu^{I}_{b}$) are the inverse temperature and chemical potential of the subsystem (bath), respectively, and $Z^{I}_{s}$ ($Z^{I}_{b}$) is the grand-canonical partition of the subsystem (bath) he2012initial :

[TABLE]

where $\epsilon_{s(b)}^{i}$ are the single-particle eigenenergies of the subsystem (bath). The total energy and particle number of the subsystem (bath) are

[TABLE]

respectively. The occupation per site in the subsystem is then $n^{I}_{s}={\mathcal{N}}^{I}_{s}/L_{s}$.

If the chemical potential and the inverse temperature of the subsystem and the bath were chosen to be same, then $\hat{\rho}^{I}$ is the thermal equilibrium state of

[TABLE]

which is the initial and final Hamiltonian of our cyclic process. Since thermal equilibrium states are passive, one would not be able to extract work in a cyclic process starting from such an initial state d2016quantum . Hence, one needs to choose $\mu^{I}_{s}\neq\mu^{I}_{b}$ and/or $\beta^{I}_{s}\neq\beta^{I}_{b}$. Next we show that, in order to be able to extract work with our quench protocol, we need $\mu^{I}_{s}<\mu^{I}_{b}$.

Since the initial state is a product of two thermal states, its entropy can be written as

[TABLE]

where $S^{I}_{s(b)}$ is the GE entropy of the subsystem (bath) he2012initial :

[TABLE]

Given the total initial energy and the total initial number of particles

[TABLE]

respectively, one can construct the density matrix $\hat{\rho}^{\text{GE}}_{I}$ of a thermal state of the entire system that matches the initial total energy and number of particles. Namely, $\hat{\rho}^{\text{GE}}_{I}=\exp[-\beta_{I}^{\prime}(\hat{\mathcal{H}}^{I}-\mu_{I}^{\prime}\hat{\mathcal{N}})]/\mbox{Tr}\,[\exp(-\beta_{I}^{\prime}[\hat{\mathcal{H}}^{I}-\mu_{I}^{\prime}\hat{\mathcal{N}}])]$, such that $E^{I}=\mbox{Tr}\,[\hat{\rho}^{\text{GE}}_{I}\hat{\mathcal{H}}^{I}]$ and ${\mathcal{N}}^{I}=\mbox{Tr}\,[\hat{\rho}^{\text{GE}}_{I}\hat{\mathcal{N}}]$. Since $V^{I}_{s}=V_{b}=0$, that state has a uniform occupation of the sites, $n_{I}={\mathcal{N}}^{I}/L$. We call the entropy corresponding to such a thermal state $S^{\text{GE}}_{I}$.

Now we can consider the strong quench of the on-site potential ($V_{s}^{I}=0\rightarrow V_{s}^{F}=V$) starting from $\hat{\rho}^{I}$ and $\hat{\rho}^{\text{GE}}_{I}$. They result in energies after the quench

[TABLE]

respectively. ($n^{I}_{s}$ and $n_{I}$ are the site occupations in the subsystem for $\hat{\rho}^{I}$ and $\hat{\rho}^{\text{GE}}_{I}$, respectively.) For $V>0$ and $n^{I}_{s}\geq n_{I}$, $E_{Q}\geq E^{\prime}_{Q}$. Since the entropy in the GE is a monotonic function of the energy, we immediately realize that equilibration to the GE results in entropies $S^{\text{GE}}_{Q}$ and $S^{{}^{\prime}\text{GE}}_{Q}$, corresponding to $E_{Q}$ and $E^{\prime}_{Q}$, respectively, which satisfy $S^{\text{GE}}_{Q}\geq S^{{}^{\prime}\text{GE}}_{Q}$. We also know that, as a result of the quench, $S^{{}^{\prime}\text{GE}}_{Q}\geq S^{\text{GE}}_{I}$, from which it follows that $S^{\text{GE}}_{Q}\geq S^{\text{GE}}_{I}$.

Finally, let us consider the quasi-static process that brings the system back to the initial Hamiltonian. In the limit $N\rightarrow\infty$, when equilibration to the grand-canonical ensemble is assumed in every weak quench, the entropy of the thermal state at the end of the quasi-static process $S^{\text{GE}}_{F}$ is $S^{\text{GE}}_{F}\simeq S^{\text{GE}}_{Q}$ aizenman1981third . Additionally, we see that, for $V>0$ and $n^{I}_{s}\geq n_{I}$, $S^{\text{GE}}_{F}\geq S^{\text{GE}}_{I}$. Again, since the entropy in the GE is a monotonic function of the energy, we conclude that the final energy of the system after the cyclic process is larger than the initial energy. As a result, no work can be extracted [$W$ in Eq. (3) is negative]. Hence, in order to be able to extract work, we need $n^{I}_{s}<n_{I}$. This can be ensured by choosing $\mu^{I}_{s}<\mu^{I}_{b}$. Note that our analysis and conclusion are completely independent of the values of $\beta^{I}_{s}$ and $\beta^{I}_{b}$. In what follows, we take $\beta^{I}_{s}=\beta^{I}_{b}=\beta^{I}$ for all our calculations. This choice plays an important role in devising the protocol that maximizes the work extracted.

III Grand Canonical Ensemble (GE)

For isolated integrable quantum systems, such as those described by quadratic Hamiltonians, the expectation values of observables after equilibration following a quench are not described by traditional ensembles of statistical mechanics. (This is true even if the initial state before the quench is a thermal equilibrium state of a quantum chaotic Hamiltonian rigol2016fundamental .) The reason for this lack of thermalization is the presence of an extensive set of nontrivial conserved quantities. However, very weak integrability-breaking interactions are expected to ensure that the system thermalizes at long times rigol_14 ; bertini_essler_15 ; brandino_caux_15 ; rigol2016fundamental ; bertini_essler_16 , even if they do not significantly change the energy of the system from their noninteracting values.

With this in mind (see also Sec. II), in this section we replace the density matrix of the entire system after a quench with that of a GE whose energy and total number of particles match those of the noninteracting system after the quench. That GE density matrix is assumed to describe observables of interest here, such as site occupations, after equilibration in the presence of weak integrability-breaking interactions. What happens in the absence of interactions is the subject of the next section. Studying quadratic Hamiltonians allows us to gain analytic insights about the specific quench protocol that saturates theoretical bounds in the thermodynamic limit. It also allows us to study numerically finite systems that can be small or as large as desired.

III.1 Quasi-static process in the thermodynamic limit

Given our definition of work extracted [see Eq. (3)], maximal work is associated with minimal energy after completing the cyclic process. Assuming the system thermalizes, this corresponds to the case in which the system has minimal GE entropy at the end of the cyclic process. Therefore, our goal in order to maximize work is to design a cyclic process that keeps the entropy constant after every quench (the entropy cannot decrease). If this is achieved, we saturate the upper bound for the work that can be extracted in a cyclic process perarnau2016work

[TABLE]

where $\hat{\rho}^{\text{GE}}_{\text{max}}$ is the density matrix of a GE that has the same entropy as the initial state, $S^{\text{GE}}_{\text{max}}=S^{I}$. (It also must have the same number of particles, but this is enforced no matter the protocol implemented.) In the limit $N\to\infty$ (quasi-static protocol), and for large systems (we note that equilibration times increase with increasing system size), the entropy at the end of the cyclic process equals that after the strong quench, i.e., $S^{\text{GE}}_{F}\simeq S^{\text{GE}}_{Q}$. Hence, all we need to do is to find a protocol by means of which the GE entropy after the strong quench is that of the initial state, $S^{\text{GE}}_{Q}=S^{I}$.

In the thermodynamic limit, when $L_{s}\to\infty$ and $L_{b}\to\infty$ for $\eta=L_{s}/L$ finite, the initial site occupations of the subsystem and the bath can be obtained as

[TABLE]

where $g(\epsilon)=(\pi\sqrt{4t^{2}-\epsilon^{2}})^{-1}$ is the density of states. Similarly, the initial energies of the subsystem and the bath read

[TABLE]

After the quench $V_{s}^{I}=0\rightarrow V_{s}^{F}=V$, the GE density matrix of the entire system is

[TABLE]

where $Z_{Q}=\mbox{Tr}\,[\exp(-\beta_{Q}[\hat{\mathcal{H}}(V,0)-\mu_{Q}\hat{\mathcal{N}}])]$. $\beta_{Q}$ and $\mu_{Q}$ are computed such that the GE energy and number of particles match the results after the quench, namely, $E_{Q}$ [see Eq. (12)] and ${\mathcal{N}}^{I}$ [see Eq. (11)], respectively.

Within the local density approximation, the GE energy of the subsystem and the bath after the quench can be obtained as

[TABLE]

respectively. We can rewrite the energy of the subsystem as

[TABLE]

Neglecting the $O(1)$ contribution of the hopping between the subsystem and the bath, the total energy

[TABLE]

[see Eqs. (11) and (12)], where $E^{I}_{s}$ and $E^{I}_{b}$ can be computed using Eq. (15), and $n^{I}_{s}$ can be computed using Eq. (14). A trivial solution to Eq. (19) is obtained for $E^{\text{GE}}_{s}=E^{I}_{s}+Vn^{I}_{s}L_{s}$ and $E^{\text{GE}}_{b}=E^{I}_{b}$, which require $\beta_{Q}=\beta^{I}$, $\mu_{Q}=\mu^{I}_{b}$, as well as $V=\mu^{I}_{b}-\mu^{I}_{s}$.

This solution trivially satisfies that the total number of fermions after the quench remains the same as before the quench, and that the entropy of the initial state

[TABLE]

and of the GE describing the system after the quench

[TABLE]

where $\Pi(x)=-x\ln x-(1-x)\ln(1-x)$, are the same. This is possible because the initial state of our quench, which is not a thermal equilibrium state of the initial Hamiltonian, is very close to a thermal equilibrium state of the Hamiltonian after the quench. Such quenches can be implemented in a wide range of settings, including interacting systems.

A straightforward example in the context of the quadratic Hamiltonian (8) is the case in which initially the subsystem and the bath have the same chemical potential $\mu^{I}_{b}=\mu^{I}_{s}=0$ but different inverse temperature ($\beta^{I}_{s}$ and $\beta^{I}_{b}$, respectively). In this case, one can extract work by quenching the hopping amplitudes in the subsystem (no quench of the on-site potentials). Maximal work can be extracted for a strong quench of the hopping amplitude in the subsystem from $t\rightarrow t_{Q}$ with $t_{Q}/t=\beta^{I}_{s}/\beta^{I}_{b}$, and then returning to $t$ using a quasi-static process ($N$ weak quenches in the subsystem).

III.2 Work extraction and entropy differences vs the number of quenches

Next, we would like to gain an understanding of what happens in finite systems and for a finite number of quenches. For this, we use numerical calculations. Since the Hamiltonian of interest here is quadratic, all observables in thermal equilibrium can be computed from the one-body density matrix, $\rho^{\text{GE}}_{ij}=\mbox{Tr}\,[\hat{\rho}^{\text{GE}}\hat{c}^{{\dagger}}_{i}\hat{c}_{j}]$, which can be obtained as

[TABLE]

where $I$ is the identity matrix, and $H$ is the matrix representing our Hamiltonian in the single-particle basis $\hat{\mathcal{H}}=\sum_{ij}\hat{c}^{{\dagger}}_{i}H_{ij}\hat{c}_{j}$ muramatsu_99 ; assaad_02 ; rigol2005finite .

Given the total energy $E_{Q}$ [see Eq. (12)], and the GE site occupancy in the subsystem $n^{\text{GE}}_{s,Q}$, after the quench $V_{s}^{I}=0\rightarrow V_{s}^{F}=V$, it is straightforward to see that the energy of the entire system after the $N$ weak quenches in which thermalization occurs, $E^{\text{GE}}_{F,N}$, depends on $N$ [in this analysis we ignore the final ($N+1$) local quench in which the subsystem and the bath are disconnected]. Denoting the site occupation in the subsystem after thermalization following the $m^{\text{th}}$ weak quench as $n^{\text{GE}}_{s,q_{m}}$, one can write

[TABLE]

Since the on-site potential of the subsystem is reduced the same amount after each weak quench, the chemical potential in the GE after the $m^{\text{th}}$ weak quench ($\mu_{q_{m}}$) exhibits a linear decrease with $m$ (the temperature decreases slightly after every weak quench): $\mu_{q_{m}}\simeq\mu_{Q}+[(\mu_{F}-\mu_{Q})m]/N$ [see Fig. 2(a)], where $\mu_{Q}$ and $\mu_{F}$ are the chemical potentials after the strong quench and at the end of the cyclic process, respectively. Given $\mu_{q_{m}}$, the site occupations in the subsystem $n^{\text{GE}}_{s,q_{m}}$ can be obtained computing an integral such as the one in Eq. (14). The relation between $\mu_{q_{m}}$ and $n^{\text{GE}}_{s,q_{m}}$ is, in general, not a linear one. However, when $|\mu_{Q}|$ and $|\mu_{F}|$ are much smaller than the bandwidth, the relation is linear and

[TABLE]

where $n^{\text{GE}}_{s,F}$ is the site occupation in the subsystem at the end of the cyclic process. This relation works remarkably well when $V\lesssim 2$ for $\mu^{I}_{b}=-\mu^{I}_{s}$ [see Fig. 2(b)]. We choose $\mu^{I}_{b}=-\mu^{I}_{s}$ so that the system is at half-filling for $L_{s}=L_{b}$.

Using Eq. (24), it is straightforward to compute the work extracted when ignoring the final local quench in which the subsystem and the bath are disconnected

[TABLE]

As shown in Fig. 3(a), this expression is in excellent agreement with the exact numerical results for $L=1000$, $V=1.0$ and $1.2$, and $\mu^{I}_{b}=-\mu^{I}_{s}=0.5$, when the final local quench in which the subsystem and the bath are disconnected is taken into account. The top curve ($V=1.0$) shows results when the strong quench fulfills the condition for maximal work extraction. The vertical dashed-dotted lines depict the predictions of Eq. (25) for $N\rightarrow\infty$.

The entropy in the GE can be written as he2012initial

[TABLE]

where $I_{\alpha}=(\exp[\beta(\epsilon_{\alpha}-\mu)]+1)^{-1}$ is the occupation of the single-particle states in the GE, and $\epsilon_{\alpha}$ are the single-particle eigenenergies. It follows from Eq. (26) that the derivative of the GE entropy at the end of the cyclic process $S^{\text{GE}}_{F}$, with respect to the total number of weak quenches $N$, is $dS^{\text{GE}}_{F}/dN=\beta_{F}\sum_{\alpha}(\epsilon_{\alpha}-\mu_{F})\,dI_{\alpha}/dN$, where the final inverse temperature $\beta_{F}$ and chemical potential $\mu_{F}$ depend on $N$. Since $\sum_{\alpha}I_{\alpha}$ equals the total number of particles, which is conserved during our cyclic process, $\sum_{\alpha}\mu_{F}\,dI_{\alpha}/dN=0$. Hence, as expected, $dS^{\text{GE}}_{F}/dN=\beta_{F}\sum_{\alpha}\epsilon_{\alpha}\,dI_{\alpha}/dN=\beta_{F}\,dE^{\text{GE}}_{F,N}/dN$.

Using Eq. (25), we see that

[TABLE]

For $N$ sufficiently large, $\beta_{F}$ and $n^{\text{GE}}_{s,F}$ are independent of $N$ [see the inset in Fig. 3(a) for the behavior of $\beta_{F}$ vs $N$], and $S^{\text{GE}}_{F}\propto 1/N$.

In Fig. 3(b), we show the difference between the GE entropy per site at the end of the cyclic process $S_{F}^{\text{GE}}/L$ and the GE entropy per site after the strong quench $S_{Q}^{\text{GE}}/L$. Figure 3(b) shows that, when the final local quench is ignored, the entropy difference vanishes with $1/N$ as predicted. When the last local quench is taken into account, the entropy difference can be seen to saturate with increasing $N$. As expected, with increasing system size, the effect of the local quench becomes negligible and the difference approaches the prediction in Eq. (27).

III.3 Work Extraction and entropy differences vs the quench parameter

Here we study the effect of changing the strong quench strength (set by the value of $V$) in the work extracted in the limit $N\to\infty$, as well as on the GE entropy after the strong quench, for finite system sizes. The energy at the end of the cyclic process for $N\to\infty$ is computed as follows: (i) We determine the entropy of the GE that has the same energy and number of particles as our system after the strong quench. (ii) We determine the GE that has the same entropy and number of particles determined in (i) but for the Hamiltonian $\hat{\mathcal{H}}(0,0)$ [see Eq. (2)], i.e., the Hamiltonian of the system after the $N\to\infty$ weak quenches. (iii) We compute the energy of the GE in (ii) after the local quench in which the subsystem and the bath are disconnected. The difference between the initial energy and the energy determined in (iii) is the work extracted in the limit $N\to\infty$.

In Fig. 4(a), we plot results for the work extracted per site vs $V/(\mu_{b}^{I}-\mu_{s}^{I})$ for three values of $\mu_{b}^{I}=-\mu_{s}^{I}$, for $L=1000$. The dashed lines are the maximal work bounds predicted by Eq. (13). The numerical results show that, as advanced for $V=\mu_{b}^{I}-\mu_{s}^{I}$, the work extracted for each value of $\mu_{b}^{I}=-\mu_{s}^{I}$ nearly saturates the maximal work bound. The inset in Fig. 4(a) shows that the small difference between the numerical result and the bound, for $\mu_{b}^{I}=-\mu_{s}^{I}=0.5$, vanishes as $1/L$ with increasing system size. Figure 4(a) also shows that increasing the difference between $\mu_{b}^{I}$ and $\mu_{s}^{I}$ increases the maximal work one can extract when the system thermalizes.

Figure 4(b) shows the difference between the GE entropy after the strong quench and the initial entropy, per site, vs $V/(\mu_{b}^{I}-\mu_{s}^{I})$. Results are shown for the same three values of $\mu_{b}^{I}=-\mu_{s}^{I}$ and system size as in Fig. 4(a). The entropy difference per site can be seen to be minimal when $V=\mu_{b}^{I}-\mu_{s}^{I}$, and the inset shows that it vanishes as $1/L$ with increasing system size.

IV Generalized Gibbs Ensemble

In this section, we study what happens when the systems are truly noninteracting. While this might appear to be a theoretical exercise of no relevance to experimental systems (or microscopic quantum devices), one should bear in mind that if interactions are very weak there exists the possibility that the dynamics for experimentally relevant time scales (operation times) is well described by a noninteracting Hamiltonian. The same can be said about systems that are interacting but close to some integrable point cazalilla_11 . For experimentally relevant time scales, their dynamics can be described by an integrable Hamiltonian and they do not thermalize, even if at very long times (not accessible in experiments) one expects that integrability-breaking effects result in thermalization. Beautiful experiments with ultracold atoms in 1D geometries have shown such a lack of thermalization kinoshita_wenger_06 ; gring_kuhnert_12 ; langen_erne_15 , while others have demonstrated that thermalization does occur in (nearly) isolated quantum systems if they are not close to integrable regimes trotzky_chen_12 ; kaufman_tai_16 ; clos_porras_16 .

As mentioned before, the breakdown of thermalization in integrable systems is due to the existence of an extensive number of nontrivial conserved quantities. In the noninteracting system of interest here, the conserved quantities ${\hat{I}_{j}}$ are the occupations of the single-particle eigenstates of the Hamiltonian after the quench (they are conserved because the particles do not interact with each other). There are as many of those as lattice sites, i.e., there is an extensive number of them. In integrable systems in general, and in our noninteracting system in particular, observables after equilibration are expected to be described by the GGE rigol2007relaxation (see Refs. vidmar_rigol_16 ; essler_fagotti_16 ; cazalilla_chung_16 ; caux_16 for recent reviews).

The GGE density matrix rigol2007relaxation , which is obtained maximizing the entropy under the constraints imposed by the conserved quantities (following Jaynes jaynes_57a ) and has been justified microscopically in terms of generalized eigenstate thermalization cassidy2011generalized ; vidmar_rigol_16 ; caux_essler_13 , can be written as

[TABLE]

where $Z_{\text{GGE}}=\mbox{Tr}\,[\exp(-\sum_{\alpha}\lambda_{\alpha}\hat{I}_{\alpha})]$ is the partition function of the GGE. The Lagrange multipliers $\lambda_{\alpha}$ are determined by the condition $\mbox{Tr}\,[\hat{\rho}^{\text{GGE}}\hat{I}_{\alpha}]=I^{I}_{\alpha}\equiv\mbox{Tr}\,[\hat{\rho}^{I}\hat{I}_{\alpha}]$, in which $\hat{\rho}^{I}$ is the density matrix of the initial (nonstationary) state. In the fermionic system of interest here rigol2007relaxation

[TABLE]

and the GGE entropy is

[TABLE]

Unlike for systems that thermalize and hence can be described by the GE, there is no simple way to determine the occupation of the single-particle eigenstates of the Hamiltonian after a quench. In addition, contrary to the GE, those occupations are not a monotonic function of the single-particle eigenenergies. As a result, within the GGE, the entropy is not necessarily a monotonic function of the energy. Therefore, the analytical arguments used in the context of the GE are not valid within the GGE.

In what follows, we report and discuss numerical results for the cyclic protocol in Fig. 1 when we replace the exact density matrix of the system after equilibration by the GGE density matrix. As shown in Fig. 5, within our protocol, numerical results for work extraction using exact dynamics (waiting random times after equilibration) and the GGE density matrix are in excellent agreement.

IV.1 Work Extraction

Given our initial state, which is a product of GE density matrices, the first question we address is the effect that the number of weak quenches $N$ has on the work extracted in the cyclic process. As mentioned before, within the GGE, the entropy is not necessarily a monotonic function of the energy. This means that there is no a priori reason to expect that a quasi-static return to the initial Hamiltonian, following the strong quench, allows us to extract the most work. Actually, for a specific non-passive initial state, in Ref. perarnau2016work a quasi-static process was shown not to be optimal for extracting work in a noninteracting fermionic system.

In Fig. 5, we show the work extracted per site within the GGE description, $W^{\text{GGE}}/L$, for the same two cyclic processes as in Fig. 2. Figure 5 shows that $W^{\text{GGE}}/L$ increases with $N$, and that, as $N\rightarrow\infty$, the work extracted for $V=1.0$ is greater than for $V=1.2$, as when thermalization occurs. The inset in Fig. 5 shows that the difference between the GGE entropy per site at the end of the cyclic process and after the strong quench, $(S^{\text{GGE}}_{F}-S^{\text{GGE}}_{Q})/L$, decreases with increasing $N$. Also, for any given $N$, the entropy difference is smaller for $V=1.0$ than for $V=1.2$.

Next, we study how changing the strength of the strong quench changes the work extracted within the GGE description for a large, but finite, number $N$ of weak quenches. We report results for $N=1000$. (Unlike for the GE, to determine the work extracted in the limit $N\to\infty$ within the GGE, one needs to do numerical calculations for finite $N$ and extrapolate the results to $N\to\infty$.) Figure 6 shows that, similarly to the results obtained for the GE, maximal work is extracted when $V=\mu^{I}_{b}-\mu^{I}_{s}$. This is understandable in terms of entropy production (or the lack thereof) as, for $V=\mu^{I}_{b}-\mu^{I}_{s}$, the strong quench in our protocol does not produce entropy within the GGE description. This follows from the inequalities $S^{\text{GE}}_{Q}\geq S^{\text{GGE}}_{Q}\geq S^{I}$. Since $S^{\text{GE}}_{Q}\simeq S^{I}$ for $V=\mu^{I}_{b}-\mu^{I}_{s}$, it then follows that $S^{\text{GGE}}_{Q}\simeq S^{I}$. It also follows that $S^{\text{GE}}_{Q}\simeq S^{\text{GGE}}_{Q}$, which, together with the fact that the energy of the GE and the GGE must match after the strong quench, hints that the GGE density matrix after the strong quench is very close to that of the GE (as both occupations of single-particle eigenstates and their ordering must match). Hence, the GGE density matrix after the strong quench is (almost) passive.

Given a density matrix $\hat{\rho}^{I}$, the work that can be extracted within a GGE description of equilibration during a cyclic process has an upper bound perarnau2016work

[TABLE]

where $\epsilon_{\alpha}$ are the single-particle energy eigenvalues of $\hat{\mathcal{H}}^{I}$ (in ascending order), and $I^{\prime}_{\alpha}$ are the occupations of single-particle eigenstates in $\hat{\rho}^{I}$ reordered in descending order so that $\sum_{\alpha}\epsilon_{\alpha}I^{\prime}_{\alpha}$ is the minimal energy given that set of occupations. Since the occupations of single-particle energy eigenstates are the same in the initial state and in the hypothetical final state with minimal energy, both states have identical entropies.

Figure 6 shows that the maximal work extracted in our cyclic process is very close to that bound. As mentioned before, for $V=\mu^{I}_{b}-\mu^{I}_{s}$ and $N\rightarrow\infty$, our protocol ensures that $S_{F}^{\text{GGE}}\simeq S^{I}$, i.e., the initial and final sets of occupations of the single-particle eigenstates of $\hat{\mathcal{H}}^{I}$ are expected to be the same. In addition, since the density matrix of the GGE after the strong quench is (almost) passive, $N\rightarrow\infty$ ensures that the density matrix of the GGE at the end of the cyclic process is (almost) passive, i.e., our final state is (almost) the hypothetical final state in Eq. (31).

We have studied what happens with the small differences seen in Fig. 6 between the numerical calculations using the GGE and the predictions of Eq. (31), $\Delta W^{\text{GGE}}/L=(W^{\text{GGE}}_{\text{max}}-W^{\text{GGE}})/L$ at $V=\mu^{I}_{b}-\mu^{I}_{s}$, when one changes the system size $L$ and the total number of small quenches $N$. In Fig. 7, we show results for $V=\mu^{I}_{b}-\mu^{I}_{s}=1.0$. As expected, $\Delta W^{\text{GGE}}/L$ decreases with increasing $N$ and $L$. With increasing $N$, the results approach a power law decay $\propto 1/L$. This confirms our expectation that as $N\rightarrow\infty$ and then $L\rightarrow\infty$, our protocol for $V=\mu^{I}_{b}-\mu^{I}_{s}$ saturates the bound in Eq. (31).

V Comparison between GE and GGE descriptions

In this section, we discuss the differences between the GE and the GGE descriptions of our cyclic process. We should stress that, after a quench starting from the same initial state, the energy and number of particles within the GE and GGE descriptions are identical. What are different are the density matrices describing the system. This is what leads to different results after subsequent quenches and, ultimately, to different work extracted after completing cyclic processes.

As we discussed in Secs. III and IV, the protocols devised to extract maximal work within the GE and GGE descriptions do not increase the entropy of the system. Since in the GE the energy is a monotonically increasing function of the entropy, and for any given energy the entropy is maximal, whenever the GGE has the same entropy as the GE it must have a higher energy. As a result, $W^{\text{GE}}_{\text{max}}\geq W^{\text{GGE}}_{\text{max}}$. This can be seen if one compares the results in Fig. 4(a) and in Fig. 6(a). A recent work has proposed a protocol to extract $W^{\text{GE}}_{\text{max}}-W^{\text{GGE}}_{\text{max}}$ in a noninteracting setting verstraelen2017unitary .

A quantity that exhibits a qualitatively different behavior in the GE and GGE with increasing the quench strength is the average site occupation in the subsystem $n_{s,Q}$ (and, consequently, in the bath). In Fig. 8, we plot $n_{s,Q}$ as a function of $V$ for three values of $\mu^{I}_{b}=-\mu^{I}_{s}$ (i.e., at half-filling). While one can see that in the GE, for $V\leq 4$, $n_{s,Q}$ decreases smoothly as $V$ increases, $n_{s,Q}$ exhibits a nonmonotonic behavior in the GGE with a sort of kink at $V=4$. For $V\geq 4$, $n_{s,Q}$ does not change in the GGE with increasing $V$. This is because for $V\geq 4$ all fermions that are in the subsystem (bath) before the quench remain in the subsystem (bath) after the quench, which is the result of the subsystem and the bath having a local on-site potential difference that is larger than the band-width, i.e., because of energy conservation the fermions cannot hop between the subsystem to the bath in the absence of interactions.

In order to be more quantitative in the comparison of the GE and the GGE, we calculate the relative differences in the site ($n^{\text{GE/GGE}}_{i,Q}$) and single-particle energy eigenstate ($I^{\text{GE/GGE}}_{\alpha,Q}$) occupations after the strong quench:

[TABLE]

Results for those two quantities are reported in Fig. 9. One can see there that $\Delta n$ [Fig. 9(a)] and $\Delta I$ [Fig. 9(b)] have deep minima (note the logarithmic scale) at $V=\mu^{I}_{b}-\mu^{I}_{s}$, which correspond to the quenches for which maximal work is extracted within both ensemble descriptions. As discussed in Sec. IV.1, for such strong quenches the GE and GGE density matrices are expected to be very close to each other. The insets in Fig. 9 show that, as expected, $\Delta n$ and $\Delta I$ for $V=\mu^{I}_{b}-\mu^{I}_{s}$ vanish with increasing system size (almost linearly with $1/L$).

In all numerical results reported so far, we considered the case in which the subsystem and the bath have the same size, $L_{s}=L_{b}$. In Fig. 10, we report results obtained when changing the ratio between the size of the subsystem and the size of the entire system, $\eta=L_{s}/L$. We focus on the protocol for which maximal work can be extracted, $V=\mu^{I}_{b}-\mu^{I}_{s}$. Figure 10 shows that, both for the GE and GGE descriptions, the work extracted per site in the entire system is maximal when $L_{s}=L_{b}$ ($\eta=1/2$). On the other hand, the inset in Fig. 10 shows that the work extracted per site in the subsystem is a monotonically decreasing function of $\eta$. Again, this is true both for the GE and GGE. Depending on whether one wants to extract the most work or the most work per site of the subsystem, one needs to select the subsystem size to be equal to that of the bath or much smaller than that of the bath, respectively.

The GE results in Fig. 10 can be understood in the context of the theoretical framework discussed in Sec. III for the thermodynamic limit. The strong quench with $V=\mu^{I}_{b}-\mu^{I}_{s}$ does not produce entropy, and is used to extract maximal work. For this quench, the chemical potential and temperature of the entire system in the GE that describes the thermalized observables are $\mu_{Q}=\mu^{I}_{b}$ and $\beta_{Q}=\beta^{I}$, respectively. The site occupations in the subsystem remain unchanged from their initial values. Hence, substituting $n^{\text{GE}}_{s,Q}=n^{I}_{s}$, $V=\mu^{I}_{b}-\mu^{I}_{s}$, and $N\to\infty$ in Eq. (25), one obtains $W^{\text{GE}}\simeq L_{s}(\mu^{I}_{b}-\mu^{I}_{s})(n^{s}_{F}-n^{I}_{s})/2$. The average site occupation in the final state is $n^{\text{GE}}_{s,F}=\eta\,n^{I}_{s}+(1-\eta)n^{I}_{b}$, where $\eta=L_{s}/L$, so the (maximal) work extracted in our cyclic process is

[TABLE]

The dashed lines in the main panel in Fig. 10 and its inset depict the results from Eq. (33), and can be seen to be in excellent agreement with the numerical results for $N\rightarrow\infty$ and $L=1000$.

In Fig. 11, we plot $W/[L(\mu^{I}_{b}-\mu^{I}_{s})]$ vs $n^{I}_{b}-n^{I}_{s}$ for two values of $\eta$ and for two initial temperatures $T^{I}=(\beta^{I})^{-1}$. Results within the GE description are shown in Fig. 11(a), while results within the GGE description are shown in Fig. 11(b). The results in Fig. 11(a) are in excellent agreement with the predictions of Eq. (33) whenever the difference in site occupations in the subsystem and the bath, as well as the initial temperature, are not too large. They allow one to also gain a qualitative understanding of what happens within the GGE description because, as seen in Fig. 11(b), the GGE results are qualitatively similar to those obtained within the GE. We note that, both in the GE and GGE, $W/[L(\mu^{I}_{b}-\mu^{I}_{s})]$ for a given value of $n^{I}_{b}-n^{I}_{s}$ decreases with increasing $T^{I}$.

VI Summary

We studied work extraction within a fermionic quadratic model in an isolated (unitarily evolving) 1D lattice system, and considered equilibration both to the GE, to describe what happens in a weakly interacting quantum chaotic system, as well as to the GGE, to describe what happens in the noninteracting limit. We considered initial states that are products of thermal states of the subsystem and the bath (in which we divided the isolated system). We devised a cyclic protocol that begins by connecting the subsystem and the bath, and quenching the local on-site potentials in the subsystem (we called that quench the “strong” quench in our protocol). After equilibration, we applied a quasi-static process in which the local on-site potentials in the subsystem are brought to the values before the strong quench by means of $N$ weak quenches, after which the subsystem and the bath are disconnected (a local quench).

We calculated the work extracted when changing the strength $V$ of the strong quench, and the number $N$ of weak quenches, both within the GE and the GGE. We found the value of $V$ for which maximal work can be extracted both in the GE and the GGE, and discussed why for that value of $V$ our cyclic protocols saturate the theoretical bounds. We studied the effect of changing the ratio between the sizes of the subsystem and the bath, but focused in the case in which they are the same. This might be of interest to understand microscopic devises, and is different from the subsystem-much-larger-than-bath approach in traditional thermodynamics. Within our cyclic protocol, and in the parameter regime studied, no qualitative differences were found in the work extracted when considering equilibration to the GE and the GGE. Exploration of initial states for which systems described by the GGE can be used to extract work without involving quasi-static processes perarnau2016work is a topic of future research.

Acknowledgements.

This work was supported by the Army Research Office Grant No. W911NF1410540. The computations were done at the Institute for CyberScience at Penn State.