A many-body heat engine at criticality

For the performance of the engine the filling fraction $N/M$ plays an important role. At $T_{C}=0$ and for an incommensurate filling, $N\neq M$, the particles are delocalized in the lattice and can move within the box. However, for a commensurate filling, $N=M$, a pinning transition occurs for any infinitesimal lattice strength, whereby each particle becomes more strongly localized at an individual lattice site, which significantly restricts its motion [54, 55]. The behaviour of this insulating phase is then determined by the energy gap in the single particle spectrum (see Fig. 1(d)) which has a size of approximately $V_{0}/2$ for shallow lattices and $2\sqrt{V_{0}}$ for deep lattices [38]. The differences in the accessible single particle excitation spectrum for $N/M$ therefore lead to different behaviours when running the engine and in Fig. 1(e,f) one can clearly see that peak performance is achieved at unit filling. At this point the particles in the cold adiabat fill the lowest energy band and as $V_{0}$ is increased the energy gap is widened. Thermal excitations induced by the hot bath then allow particles to jump the gap and therefore more energy can be extracted as the lattice depth is decreased along the hot adiabat. As the high performance regime is the interesting one for heat engines, we will focus on the case of unit filling in the following.

The advantage of exploiting the critical point in a many-body quantum heat engine (mQHE) with $N$ particles can be quantified by comparison with an ensemble of $N$ single-particle quantum heat engines (sQHE), see Fig. 1(b). Each sQHE obeys the Hamiltonian given by Eq. (1) with a box length of $L=\pi/k_{0}$, so that exactly one lattice well is present, $V_{l}=V_{0}\cos(k_{0}x+\pi/2)$. The Otto cycle is then carried out using the same lattice height and bath temperatures, however in the mQHE the final state is strongly influenced by the presence of the interparticle interactions and the periodicity of the optical lattice. To quantify the difference between the mQHE and the sQHE, we calculate the ratio of their respective efficiencies and powers

such that $\eta^{*}(N)>1$ and $P^{*}(N)>1$ indicate that the many-body state gives a performance boost [27].

In Fig. 2(a) we show the efficiency ratio for an adiabatic cycle as a function of the lattice depths $V_{i}$ and $V_{f}$. For $V_{f}>V_{i}$ the cycle produces positive work and therefore acts as an engine. One can see that large many-body cooperative effects can be achieved in the regime where both lattices are weak, and where therefore the particles in the mQHE are still partially overlapping. This results in a non-trivial, non-flat single-particle energy spectrum (see Fig. 1(d)) and therefore in enhanced efficiency and power output over the sQHE. When both lattices are deep, $V_{0}\gtrsim 30$, the particles are highly localized in individual lattice sites and the single particle energy spectrum becomes degenerate forming flat bands. In this limit all many-body cooperative effects are lost and the mQHE becomes equivalent to the sQHE. Since for weak initial lattice depths, $V_{i}\lesssim 10$, the mQHE shows enhanced performance for a range of values of $V_{f}$, we will focus on this region of the parameter space in what follows, specifically considering the limiting case of initially having free particles ($V_{i}\rightarrow 0$).

Indeed, for mQHE cycles which operate at low reservoir temperatures, $\theta=\frac{E_{R}}{k_{B}T_{H}}\sqrt{V_{f}}>1$, and which ramp to deep lattices, $V_{f}\gg 1$, it is possible to find an approximate expression for the many-body performance boost,

where $\Delta=2\sqrt{V_{f}}-1$ is the energy gap (see Supplemental Material for details). From this one can immediately see that at unit filling the mQHE will always outperform the sQHE once $V_{f}>4$. Furthermore, increasing the number of particles and reaching the state of double filling, $N=2M$, where the two lowest states of each lattice site are occupied, does not lead to improved performance. In the limits $\theta,M\rightarrow\infty$ the efficiency ratio can be written as $\eta^{*}(2M)\equiv\frac{\eta(2M)}{\eta(1)}\rightarrow\frac{1-4(\Delta-1)^{-1}}{1-3\Delta^{-1}}$, showing that for $V_{f}>1$ the efficiency of the sQHE is always larger than the mQHE at double filling, which is due to the anharmonicity of the individual lattice sites leading to reduced gaps between higher lying energy states.

In Fig. 2(b) we show the numerically obtained values of the ratios $\eta^{*}_{AD}(N)$ and $P^{*}_{AD}(N)$ as a function of $V_{f}$ in comparison to the approximation in Eq. (4). One can see that the exact ratios peak at lower values of $V_{f}$, which is due to the fact that the particles in the many-body state are still partially delocalized and therefore many-body cooperativity is stronger. For deeper lattices, $V_{f}\gtrsim 50$, both ratios head towards one, as stronger localisation makes the lattice sites become effectively independent. The decay of the many-body advantage is well described by the approximation in Eq. (4) (solid yellow line in Fig. 2(b)) and given by a $1/\Delta$ dependence. While the decay is universal, the position and height of the maximum are depending on the other parameters of the system, in particular $T_{H}$ (see inset of Fig. 2(b)). In general, a significant many-body advantage exists by operating the mQHE in weak lattices and at low temperatures when the commensurate system remains close to the quantum critical point. At higher temperatures the existing thermal energies diminish the importance of energy gap and the quantum criticality is washed out. To demonstrate this we show in Fig. 2(c) the efficiency at maximum power (optimised over the lattice depth $V_{f}$, see inset in Fig. 2(c)) for two different temperatures of the cold adiabat: one deep in the quantum regime $T_{C}=0.01$, which ensures that the system is in its ground state and therefore close to the quantum critical point; the other at a slightly higher temperature $T_{C}=0.1$, where the effect of the quantum criticality is reduced. When $T_{H}$ is small, the mQHE with the lower $T_{C}$ can be seen to be more efficient and close to the Curzon-Ahlborn efficiency, $\eta_{CA}=1-\sqrt{T_{C}/T_{H}}$, which is a good indicator of the performance of the Otto cycle [56, 29, 46]. Furthermore, it is worth noticing that at higher $T_{C}$ the mQHE is outperformed by the sQHE as the thermal energy leads to less localisation within the box potential and the energy gap is washed out. Therefore, this quantum critical mQHE only shows enhanced performance at low temperatures when the system is close to the critical point.

In deep lattices the power and efficiency ratios are equivalent for any number of lattice sites at unit filling ($N=M$) and they are exactly described by Eq. (4) (see Fig. 2(d)). As the many-body advantage is proportional to $(1-1/N)$, one can see a rapid increase in both quantities for increasing particle number until $N\sim 10$, after which it asymptotically approaches $1+\left(\Delta-3\right)^{-1}$ in the thermodynamic limit. In more shallow lattices the efficiency and power ratios asymptotically reach different, but overall larger values, while the dependence on $N$ remains consistent with the behaviour observed for deep lattices (see inset of Fig. 2(d)). Indeed, one does not need to create large many-body states to see a marked improvement in engine performance, rather only a few dozen particles are sufficient for observing the effects of many-body cooperativity in this system.

While all the results above are obtained for a reversible cycle with undergoes adiabatic dynamics, this results in negligible power output due to the long timescales for each cycle, therefore necessitating fast engine cycles for finite power-output. However, fast driving through a critical point will inevitably result in non-adiabatic dynamics and irreversible work being produced, with the latter being defined as the difference between the average work of the non-adiabatic ($NA$) and adiabatic driving, $\langle W_{irr}\rangle=\langle W\rangle_{NA}-\langle W\rangle_{AD}$. This ultimately leads to reduced performance of the QHE [57, 58, 59]. To explore the effect of irreversible work on the engine performance we numerically calculate the unitary dynamics of the single particle states during the compression and expansion strokes, $\psi_{n}(x,t)=e^{-\frac{i}{\hbar}\int_{0}^{t}H(t^{\prime})dt^{\prime}}\psi_{n}(x,0)$, describing the insertion and removal of the optical lattice over a finite time $t_{f}$ (in units of $2\pi/E_{R}$). We parametrise the lattice strength as $V_{\lambda}(t)=\lambda(t)V_{f}\cos^{2}(k_{0}x)$, with $\lambda(t)=t^{3}/t_{f}^{3}[1+3(1-t/t_{f})+6(1-t/t_{f})^{2}]$, which is sufficient to explore the dynamical properties of a finite time engine stroke, but is not necessarily optimal for the system [60, 61, 62, 63]). As our focus is on the non-adiabatic dynamics initiated by the lattice ramp we will neglect the dynamics during the coupling to the different heat baths and assume that the thermalization times $t_{2}$ and $t_{4}$ are much shorter than the times for the work strokes $t_{1}$ and $t_{3}$ [42, 64]. Taking $t_{1}=t_{3}\equiv t_{f}$, the total time for the cycle is $\tau\approx 2t_{f}$.

For a finite time cycle at commensurate filling the efficiency only slowly approaches the adiabatic efficiency (see Fig. 3(a)) due to the large amount of irreversible work created when driving the system at the pinning transition (see Fig. 3(b-c)). In comparison, incommensurate fillings produce significantly less irreversible work as excitations are far from the energy gap and therefore the adiabatic efficiency can be reached for significantly shorter ramp times. Also note that more irreversibility is created during the raising of the barrier compared to the lowering, as the opening of the energy gap adds to the nonequilibrium excitations.

Even with the advantage gained from the energy gap, the resulting irreversible dynamics on short timescales set a limit on the performance of the engine cycle. In fact, this problem does not just appear in dynamics about a critical point, but is present in any non-adiabatic driving of quantum heat engines. To improve engine performance on finite timescales different shortcut to adiabaticity (STA) approaches have been suggested [42, 65, 28, 66, 43, 46, 67]. However, while STAs have been successfully developed for non-interacting and mean-field systems, designing them for strongly interacting many-body systems poses new challenges when scale invariance can not be exploited [50], and is especially difficult due to the orthogonality catastrophe in larger systems [48, 68]. We therefore employ a variational approach which can find the optimal driving amplitude $V^{\text{STA}}_{n}(t)$ for each of the single particle functions $\psi_{n}(x,t)$ which are used in the Slater determinant to construct the many-body state [51]. However, while this in principle can optimise the dynamics of each $\psi_{n}(x,t)$ individually, in practice a single lattice ramp must act on the entire many-body state and the chosen $V^{\text{STA}}_{n}(t)$ may create unwanted excitations in different $\psi_{m}(x,t)$, for $m\neq n$.

We therefore consider two different approximate STAs to optimize the many-body dynamics. First, we choose the average of the STA pulses for all states up to the energy gap, $\langle V^{\text{STA}}(t)\rangle=\sum_{n=1}^{M}V^{\text{STA}}_{n}(t)/M$, and numerically time evolve the single particle states with this finite time ramp. While this shows an improvement over the non-optimized ramp $V_{\lambda}(t)$ for all timescales (see Fig. 3(d-e)), it is only marginal as the optimization is averaged over the whole system. We therefore also consider the ramp $V^{\text{STA}}_{M-1}(t)$, which specifically optimizes the most irreversible single particle state, $\psi_{M-1}(x)$, which sits just below the gap and possesses the most excess energy after the $V_{\lambda}$ ramp (see insets in Fig. 3(d-e)). This STA results in a larger efficiency gain as excitations of this state are mostly suppressed, and the adiabatic limit is quickly reached when the lattice is weak. However, this STA becomes ineffective for fast cycles, as large modulations in the approximate STA ramp can induce excitations in the rest of the system, which is a limitation of using these approximate techniques to design STAs for many-body states.

In following we will describe how the efficiency and power ratios, $\eta_{AD}^{*}(N)$ and $P_{AD}^{*}(N)$, can be approximated in the deep lattice ($V_{f}\gg 1$) and low temperature ($\theta=\frac{E_{R}}{k_{B}T_{H}}\sqrt{V_{f}}>1$) regime as given in Eq. 4 in the main text.

For this we first look at the compression stroke for a system with $N$ particles at zero temperature ($T_{C}=0$). At the beginning of the compression stroke the particles fill a box of length $L$ which has single particle energies given by $\epsilon^{B}_{n}=(n+1)^{2}$. The groundstate energy of the $N$-particle system is therefore simply

At the end of the compression stroke the lattice has a depth $V_{f}$ and at unit filling the lowest band is fully populated with one particle per lattice site. We can therefore treat this as $N$ isolated single particles confined to individual wells whose potential can be approximated by

Treating the quartic term as a perturbation then gives the single particle energy spectrum in one deep lattice site as

where $n=\{0,1,2,\dots\}$. Here the first term is the harmonic oscillator energy, $\hbar\omega=2\sqrt{VE_{R}}$ with the rescaled $V_{f}=VE_{R}$, and the second term is the correction due to the anharmonicity of the lattice site. The energy gap is the difference between two lowest single particle states, $\Delta=\epsilon^{l}_{1}-\epsilon^{l}_{0}=2\sqrt{V_{f}}-1$ and the total energy of $N$ particles in the lowest energy band is

Next we calculate the energy of the many-body state during the expansion stroke at temperature $T_{H}$. At the start of this stroke the lattice depth is fixed at $V_{f}$ and the partition function is well approximated by that of the harmonic oscillator

which is justified when $\sqrt{V_{f}}\gg 1/4$. Alternatively, one may use the partition function as calculated from the energies given in Eq. (7), resulting in

which will yield qualitatively similar results to the low temperature harmonic oscillator approximation, however leads to more complex expressions. The total energy of the $N$ particle system is therefore $N$-times the single particle energy, $\langle H_{T_{H}}(V_{f})\rangle=\frac{N}{\mathcal{Z}}\sum_{m=0}^{\infty}\epsilon^{l}_{m}e^{-(2m+1)\theta}$. This gives

where the first term describes the thermal state of a harmonic oscillator and the second term is the correction due to the anharmonicity of the lattice site.

Finally, at the end of the expansion stroke the lattice is removed and we must describe the thermal state of $N$ particles in the box potential. As the thermal statistics of the particles are still described by the lattice band structure we need to evaluate the total energy as a distribution over the each band. Since each band can contain $N$ particles which have single particle energies described by $(n+1)^{2}$, the total energy of the $m^{th}$ band is given by

Therefore, the total energy for $N$ particles is given by $\langle H_{T_{H}}(0)\rangle=\frac{1}{\mathcal{Z}}\sum_{m=0}^{\infty}\tilde{\epsilon}^{B}_{m}e^{-(2m+1)\theta}$ which gives

The average work done during the compression and expansion strokes can then be calculated from $\langle W_{C}\rangle=\langle H_{0}(V_{f})\rangle-\langle H_{0}(0)\rangle$ and $\langle W_{H}\rangle=\langle H_{T_{H}}(0)\rangle-\langle H_{T_{H}}(V_{f})\rangle$ respectively, while the heat exchange with the hot bath is $\langle Q_{H}\rangle=\langle H_{T_{H}}(V_{f})\rangle-\langle H_{0}(V_{f})\rangle$. With these quantities the efficiency and power of the adiabatic mQHE can be calculated as in the main text.

Similarly the energies for the sQHE cycle can be straightforwardly calculated, with the energy of a single particle in the lattice potential simply given by $\langle H_{0}^{1}(V_{f})\rangle=\langle H_{0}(V_{f})\rangle/N$ and $\langle H_{T_{H}}^{1}(V_{f})\rangle=\langle H_{T_{H}}(V_{f})\rangle/N$, while in the box potential they are

As above, the efficiency and power of the sQHE can then be calculated.

Finally, to compare the performance of the mQHE to the sQHE we calculate the ratios of the efficiency $\eta^{*}=\eta(N)/\eta(1)$ and power $P^{*}=P(N)/(NP(1))$. Using the above approximations for an engine operating in the deep lattice and low temperature regime we find that we can write both ratios as

To design an STA for the time-dependent ramp acting on the whole system we use a variational approach [66, 50, 51]. This method relies on minimizing the Lagrangian and finding the optimal $V_{n}(t)$ ramp for each single particle state individually. The success of this approach depends strongly on the choice of the ansatz for the time evolution of the corresponding single particle state. For our work we choose a simple ansatz of the form of a superposition between the the initial and the target state [50]

where $\psi_{n}^{I}(x,t)$ is the initial $n^{th}$ single particle state and $\psi_{n}^{F}(x,t)$ is the corresponding target state, with $b(t)$ being a dynamical phase. We use the time dependent parameter $\varepsilon(t)$ to switch the single particle state from its initial to the target state. To be able to obey the boundary conditions $\varepsilon(0)=0$, $\varepsilon(t_{f})=1$, and $\dot{\varepsilon}(0)=\dot{\varepsilon}(t_{f})=\ddot{\varepsilon}(0)=\ddot{\varepsilon}(t_{f})=0$, we choose a functional form of the switching parameter of $\varepsilon(t)=\sum_{j=0}^{5}a_{j}t^{j}$.

The resulting optimised lattice ramp is then given by

where the variables are given by

Here $\xi^{2}$ is the width of the single particle state, $\alpha$ is the contribution from the lattice, $\beta$ the kinetic energy and $b$ is the chirp.