Quantum supremacy of many-particle thermal machines

We consider the working susbtance to be described by a quantum many-body Hamiltonian

$$\begin{eqnarray}&&H=\displaystyle \sum _{i=1}^{N}\left[-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}+\displaystyle \frac{1}{2}m\omega {(t)}^{2}{{\bf{r}}}_{i}^{2}\right]+\displaystyle \sum _{i\lt j}V({{\bf{r}}}_{i}-{{\bf{r}}}_{j}),\end{eqnarray} \tag{ 1 }$$

where $N$ is the total number of particles and $\omega (t)$ is the trap frequency of the isotropic harmonic confinement. We assume $V({\bf{r}}/b)={b}^{2}V({\bf{r}})$ so that the unitary dynamics generated by (1) during the expansion and compression strokes exhibits scale-invariance [22]. The time-evolution of an eigenstate Ψ at t = 0 is [23, 24]

$$\begin{eqnarray}&&{\rm{\Psi }}({{\bf{r}}}_{1},\ldots ,{{\bf{r}}}_{N},t)={ \mathcal N }{{\rm{e}}}^{{\rm{i}}\phi }{\rm{\Psi }}\left(\displaystyle \frac{{{\bf{r}}}_{1}}{b},\ldots ,\displaystyle \frac{{{\bf{r}}}_{N}}{b},t=0\right),\end{eqnarray} \tag{ 2 }$$

where the time-dependent phase reads $\phi ={\sum }_{i}m\dot{b}{{\bf{r}}}_{i}^{2}/(2{\hslash }b)$ and we have neglected a global dynamical phase which plays no role in our analysis. The normalization constant is ${ \mathcal N }={b}^{-{Nd}/2}$ for a system in d spatial dimensions. The scaling factor $b=b(t)\gt 0$ is the solution of the Ermakov differential equation $\ddot{b}+\omega {(t)}^{2}b={\omega }_{0}^{2}/{b}^{3},$ where we choose the constant ${\omega }_{0}=\omega (0)$ to be the trap frequency at the beginning of the stroke. As shown in appendix A.1, following a modulation of the trap frequency the mean nonadiabatic energy reads $\langle H(t)\rangle ={Q}^{* }(t){\langle H(t)\rangle }_{{\rm{ad}}},$ where

$$\begin{eqnarray}&&{Q}^{* }(t)={b}_{{\rm{ad}}}^{2}\left[\displaystyle \frac{1}{2b{(t)}^{2}}+\displaystyle \frac{\omega {(t)}^{2}b{(t)}^{2}}{2{\omega }_{0}^{2}}+\displaystyle \frac{\dot{b}{(t)}^{2}}{2{\omega }_{0}^{2}}\right],\end{eqnarray} \tag{ 3 }$$

is the nonadiabatic factor that accounts for the amount of energy excitations over the adiabatic dynamics, and ${b}_{{\rm{ad}}}={[{\omega }_{0}/\omega (t)]}^{1/2}$ is the scaling factor in the adiabatic limit, obtained from the Ermakov equation by setting $\ddot{b}(t)\approx 0$. In this limit, ${Q}^{* }(t)$ reduces to unity and the mean energy is given by ${\langle H(t)\rangle }_{{\rm{ad}}}=\langle H(0)\rangle /{b}_{{\rm{ad}}}^{2}$. From the stationarity of the initial state it follows that $b(0)=1$ and $\dot{b}(0)=0$. Remarkably, ${Q}^{* }(t)$ coincides with the nonadiabatic factor introduced by Husimi for the single-particle time-dependent harmonic oscillator [25], as we show in appendix A.2. As a result, at the end of the strokes decoupled from the heat reservoirs, ${\langle H\rangle }_{B}={Q}_{{AB}}^{* }{b}_{{\rm{ad}}}^{-2}{\langle H\rangle }_{A}$ and ${\langle H\rangle }_{D}={Q}_{{CD}}^{* }{b}_{{\rm{ad}}}^{2}{\langle H\rangle }_{C},$ where ${b}_{{\rm{ad}}}={({\omega }_{1}/{\omega }_{2})}^{1/2}$ and ${Q}_{{AB}({CD})}^{* }$ is the nonadiabatic factor for the compression (expansion) stroke. It follows that the total input work per cycle $\langle W\rangle =\langle {W}_{1}\rangle +\langle {W}_{3}\rangle$ reads

$$\begin{eqnarray}&&\langle W\rangle =\left({Q}_{{AB}}^{* }\displaystyle \frac{{\omega }_{2}}{{\omega }_{1}}-1\right){\langle H\rangle }_{A}+\left({Q}_{{CD}}^{* }\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}}-1\right){\langle H\rangle }_{C}.\end{eqnarray} \tag{ 4 }$$

From a practical point of view, the output power of the engine is more interesting as it accounts for the work per cycle $P\equiv -\langle W\rangle /({\tau }_{1}+{\tau }_{2}+{\tau }_{3}+{\tau }_{4}),$ where ${\tau }_{i},i=1,\ldots ,4$ corresponds to the duration of each stroke of the Otto cycle. Similarly, the efficiency of the many-particle QHE run in finite-time is found to be

$$\begin{eqnarray}&&\eta =1-\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}}\left(\displaystyle \frac{{Q}_{{CD}}^{* }{\langle H\rangle }_{C}-\tfrac{{\omega }_{2}}{{\omega }_{1}}{\langle H\rangle }_{A}}{{\langle H\rangle }_{C}-{Q}_{{AB}}^{* }\tfrac{{\omega }_{2}}{{\omega }_{1}}{\langle H\rangle }_{A}}\right).\end{eqnarray} \tag{ 5 }$$

The maximum efficiency is achieved under slow driving in the adiabatic limit when the QHE operates at vanishing output power P as a result of the requirement for a long cycle time $\tau ={\sum }_{i=1}^{4}{\tau }_{i}$. In this limit (${Q}_{{AB}({CD})}^{* }\to 1$ or equivalently $| \dot{\omega }/{\omega }^{2}| \to 0$, see [26]) equation (5) reduces to the Otto efficiency ${\eta }_{O}=1-{\omega }_{1}/{\omega }_{2}$ which is shared as an upper bound by the single-particle quantum Otto cycle. By contrast, realistic engines operating in a finite time achieve a finite output power at the cost of introducing nonadiabatic energy excitations that represent quantum friction. The engine efficiency is then reduced as ${Q}_{{AB}({CD})}^{* }\geqslant 1$,

$$\begin{eqnarray}&&\eta \leqslant {\eta }_{{\rm{nad}},O}\equiv 1-{Q}_{{CD}}^{* }\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}}.\end{eqnarray} \tag{ 6 }$$

This bound is independent of the nature of the working medium (number of particles and inter-particle interactions) and is tighter than the Otto efficiency bound $\eta \leqslant {\eta }_{O}$. Equations (4), (5), and (6) are the first results of this article giving explicit formulae for the total output work and efficiency, and showing the fundamental bound on the efficiency of a nonadiabatic QHE with the wide family of systems (1) as a working medium [23, 24]. We have recently discussed the optimization of a many-particle QHE using shortcut to adiabaticity [27]. In what follows we focus on enhancing the performance of the QHE without external control and derive the conditions (number of particles, interaction strength, trap frequencies, as well as the temperatures of the hot and cold reservoirs) for quantum supremacy.

To explore the interplay between the nonadiabatic dynamics and many-body effects, we consider a working medium consisting of $N$ bosons confined in an harmonic trap and interacting through an inverse-square pairwise potential [28, 29]

$$\begin{eqnarray}&&H=\displaystyle \sum _{i=1}^{N}\left[-\displaystyle \frac{{{\hslash }}^{2}}{2m}\displaystyle \frac{{\partial }^{2}}{\partial {z}_{i}^{2}}+\displaystyle \frac{1}{2}m\omega {(t)}^{2}{z}_{i}^{2}\right]+\displaystyle \frac{{{\hslash }}^{2}}{m}\displaystyle \sum _{i\lt j}\displaystyle \frac{\lambda (\lambda -1)}{{({z}_{i}-{z}_{j})}^{2}},\end{eqnarray} \tag{ 7 }$$

where $\lambda \geqslant 0$ is the coupling strength of the inter-particle interaction. This instance of (1) is the Calogero–Sutherland gas and describes an ideal Bose gas for $\lambda =0$ as well as hard-core bosons (Tonks–Girardeau gas) for $\lambda =1$ [30, 31]. The thermodynamics of the latter is identical to that of polarized fermions. For other values of λ, the Hamiltonian (7) describes hard-core bosons with inverse-square interactions and can be reinterpreted as an ideal gas of particles obeying generalized-exclusion statistics [32–34]. While the eigenstates are a paradigm of strong correlations with Bijl–Jastrow form, the spectrum is purely linear [28, 35], $E(\{{n}_{j}\})=({\hslash }\omega N/2)[1+\lambda (N-1)]+{\sum }_{j=1}^{\infty }j{\hslash }\omega {n}_{j},$ where n_j is the occupation number in the jth mode of energy $j{\hslash }\omega$, satisfying the normalization ${\sum }_{j=1}^{\infty }{n}_{j}=N$. To determine the output power and efficiency of the QHE, we first note that the partition function is given by ${Z}_{N}^{(\lambda )}={{\rm{e}}}^{-\beta \lambda \tfrac{{\hslash }\omega }{2}N(N-1)}{Z}_{N}^{(0)}$ where the partition function for an ideal Bose gas ${Z}_{N}^{(0)}={{\rm{e}}}^{-\beta {\hslash }\omega /2}{\prod }_{k=1}^{N}{(1-{{\rm{e}}}^{-\beta k{\hslash }\omega })}^{-1}$ can be computed via recurrence relations [36], see appendix B.1 for further details. Knowledge of the partition function allows one to compute the equilibrium mean thermal energy via the identity $\langle H\rangle =-{\partial }_{\beta }$ ${\rm{ln}}{Z}_{N}^{(\lambda )}$, i.e., $\langle H\rangle$ $=\tfrac{{\hslash }\omega }{2}N[1+\lambda (N-1)]$ $+{\hslash }\omega {\sum }_{k=1}^{N}k{({{\rm{e}}}^{\beta k{\hslash }\omega }-1)}^{-1}.$ In the adiabatic limit, ${Q}^{* }(t)=1$, both the output power per cycle and the efficiency become independent of the interaction strength λ, see appendix B.2. By contrast, for an arbitrary nonadiabatic driving protocol with ${Q}^{* }(t)\gt 1$, they both decrease monotonically as a function of λ. It follows that the output power is optimal for an ideal Bose gas. As a paradigm for nonadiabatic effects, we shall consider in the following a sudden-quench driving of the trap frequency between ${\omega }_{1}$ and ${\omega }_{2}$. In this case, the nonadiabatic factor takes the form ${Q}_{{AB}({CD})}^{* }=({\omega }_{1}^{2}+{\omega }_{2}^{2})/(2{\omega }_{1}{\omega }_{2})$ that is symmetric with respect to the exchange ${\omega }_{1}\leftrightarrow {\omega }_{2}$. As a result, the efficiency (5) of a realistic QHE with a short time per cycle has as fundamental upper limit (6), ${\eta }_{\mathrm{sq}}\leqslant {\eta }_{{\rm{nad}},O}=(1-{({\omega }_{1}/{\omega }_{2})}^{2})/2$, i.e., not higher than 50%. Despite this limit, we show next that nonadiabatic quantum effects can enhance the efficiency of a multi-particle QHE over the single-particle counterpart.

We determine the conditions for optimal output power ${P}_{{\rm{sq}}}^{(N,\lambda )}$ by optimizing the ratio ${\omega }_{1}/{\omega }_{2}$ for fixed ${\omega }_{1}$, temperatures ${\beta }_{{\rm{c}}},{\beta }_{{\rm{h}}}$, number of particles $N$, and interaction strength λ. We denote the efficiency at optimal output power by ${\eta }_{{\rm{sq}}}^{(N,\lambda )}$. To characterize many-particle effects, we introduce the following ratios

$$\begin{eqnarray}&&{r}_{{\rm{sq}}}^{(N,\lambda )}\equiv \displaystyle \frac{{P}_{{\rm{sq}}}^{(N,\lambda )}}{N\times {P}_{{\rm{sq}}}^{(1,\lambda )}},\ {\rho }_{{\rm{sq}}}^{(N,\lambda )}\equiv \displaystyle \frac{{\eta }_{{\rm{sq}}}^{(N,\lambda )}}{{\eta }_{{\rm{sq}}}^{(1,\lambda )}},\end{eqnarray} \tag{ 8 }$$

where the first (second) ratio compares the optimal output power (efficiency at optimal output power) of a $N$-particle QHE with that of $N$ single-particle QHEs. Their dependence on the interaction strength is displayed in figure 2, where it is shown that they can be greater than 1 for $\lambda \lt 1/2$ and optimal for $\lambda =0$. In our following analysis, we will show that for some temperature regimes, the single-particle QHE exhibits classical performance up to first order corrections, while the many-particle QHE can retain a leading quantum term of the same magnitude as the single-particle classical one. In this regime, the two ratios defined in equation (8) compare the quantum many-particle performance with that of an ensemble of classical single-particle QHEs. We consider that a QHE exhibits quantum supremacy when both ratios given in equation (8) are greater than one, see figure 3. We have verified that the conditions for quantum supremacy do not substantially depend on the definition of (8). In particular, a similar definition in which the efficiency and power of both single- and many-particle QHEs are evaluated at the same frequency ratio yields essentially the same conditions for quantum supremacy, see appendix B.3. To analyze the optimal output power we resort to the numerically exact optimization shown in figures 2 and 3 where the ratios in equation (8) are plotted for different values of $N$ and λ. The performance of a many-particle QHE is shown to depend strongly on the regime of temperature, particle number, and interaction strength.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the high-temperature limit of the hot reservoir ${\sigma }_{{\rm{h}}}\equiv N{\hslash }{\beta }_{{\rm{h}}}{\omega }_{2}\ll 1$ and of the cold reservoir ${\sigma }_{{\rm{c}}}\equiv N{\hslash }{\beta }_{{\rm{c}}}{\omega }_{1}\ll 1$, we find that the optimal frequency is given by ${({\omega }_{1}/{\omega }_{2})}_{\mathrm{sq}}\approx {a}^{1/4}\left[1+\tfrac{1}{8}{g}_{N}(\lambda -\tfrac{1}{2}){\sigma }_{{\rm{c}}}\right]$ where ${g}_{N}=(N-1)/N$ and $a\equiv {\beta }_{{\rm{h}}}/{\beta }_{{\rm{c}}}$ denotes the ratio of temperatures (see appendix D for details), whence it follows that

$$\begin{eqnarray}&&{\eta }_{\mathrm{sq}}^{(N,\lambda )}=\displaystyle \frac{1-\sqrt{a}}{2+\sqrt{a}}+\displaystyle \frac{(3-2{a}^{\tfrac{3}{4}})\sqrt{a}}{4{(2+\sqrt{a})}^{2}}{g}_{N}\left(\displaystyle \frac{1}{2}-\lambda \right){\sigma }_{{\rm{c}}}.\end{eqnarray} \tag{ 9 }$$

As a result, we recognize the leading term in ${\eta }_{\mathrm{sq}}^{(N,\lambda )}$, independent of ${\sigma }_{{\rm{c}}}$, as the Rezek–Kosloff efficiency shared by single-particle heat engines under sudden driving [11]. The leading quantum ${\sigma }_{{\rm{c}}}$-contribution to the efficiency at optimal work increases with the number of particles $N$ for $\lambda \lt 1/2$. On the contrary, if $\lambda \gt 1/2$ quantum effects lower the efficiency and vanish for $\lambda =1/2$ or $N=1$. A similar result can be found for the optimal output power. We note that for an adiabatic protocol quantum corrections of this order vanish and the engine operates at the Curzon–Ahlborn efficiency [37], see details in appendix C. Equation (9) suggests that for a QHE with a small particle number, quantum corrections can lead to quantum supremacy. Yet, ${r}_{{\rm{sq}}}^{(N,\lambda )}$ and ${\rho }_{{\rm{sq}}}^{(N,\lambda )}$ vary as $1+O({\sigma }_{{\rm{c}}})$ where the corrections depend strongly on λ and $N$. Numerically, one can see that quantum deviations predicted in equation (9) amount to a few percents with respect to the single-particle case.

In the opposite limit of very low temperature ${\sigma }_{{\rm{c}}}\gg N$ of the cold reservoir (keeping at high temperature the hot reservoir ${\sigma }_{{\rm{h}}}\ll N$), the thermal energy of the equilibrium state A does not contribute to the total work. One has ${({\omega }_{1}/{\omega }_{2})}_{\mathrm{sq}}\approx {({\kappa }_{N,\lambda }^{2}{\hslash }{\beta }_{{\rm{h}}}{\omega }_{1})}^{1/4}$ (see appendix D for details), where ${\kappa }_{N,\lambda }=\sqrt{(1+(N-1)\lambda )}$. The efficiency at optimal output power reads

$$\begin{eqnarray}&&{\eta }_{\mathrm{sq}}^{(N,\lambda )}=\displaystyle \frac{1-{\kappa }_{N,\lambda }\sqrt{{\hslash }{\beta }_{{\rm{h}}}{\omega }_{1}/2}}{2+{\kappa }_{N,\lambda }\sqrt{{\hslash }{\beta }_{{\rm{h}}}{\omega }_{1}/2}}.\end{eqnarray} \tag{ 10 }$$

For $N=1$, equation (10) reduces to the efficiency reported for a single-particle quantum Otto cycle [8]. In this regime, the efficiency at optimal output power decreases drastically with respect to the high temperature regime, see equation (9), because ${\hslash }{\beta }_{{\rm{h}}}{\omega }_{1}\gg {\beta }_{{\rm{h}}}/{\beta }_{{\rm{c}}}$ by assumption. Therefore, it is not worth running the QHE in the very low temperature limit of the cold reservoir because the working substance is effectively 'frozen' in the ground state, reducing the output work per cycle. Numerical simulations show that a similar conclusion holds in the regime ${\sigma }_{{\rm{c}}}\sim N$ and for ${\sigma }_{{\rm{c}}},{\sigma }_{{\rm{h}}}\geqslant N$.

We have seen that estimate (9) is valid for small number of particles (e.g., $N\sim 1-30$ for ${\sigma }_{{\rm{c}}}/N=0.01$) and values of the ratio ${\beta }_{{\rm{h}}}/{\beta }_{{\rm{c}}}$ consistent with the assumption ${\sigma }_{{\rm{h}}}\ll 1$. In this regime, the enhancement of the efficiency is limited to few percents of the single-particle case. In addition, for very low temperatures, where one could expect significant quantum corrections, the efficiency decreases drastically. Yet, many-particle quantum thermodynamics exhibits a novel, intermediate regime accessible for a large number of particle and characterized by $1\leqslant {\sigma }_{{\rm{c}}}\ll N$. This is a low temperature regime but admitting enough thermal excitations to prevent the working medium from being effectively 'frozen' in a single many-particle eigenstate (the ground state). In this regime, the energy ${\langle H\rangle }_{A(C)}\approx (N/{\beta }_{{\rm{c(h)}}}){\mu }_{\lambda }({\sigma }_{{\rm{c(h)}}})$ where ${\mu }_{\lambda }(\sigma )={\sigma }^{-1}{\int }_{0}^{\sigma }s{({{\rm{e}}}^{s}-1)}^{-1}{\rm{d}}s+\lambda \sigma /2$, ${\sigma }_{{\rm{c(h)}}}=N{\hslash }{\beta }_{{\rm{c(h)}}}{\omega }_{1(2)}$, and ${\mu }_{\lambda }({\sigma }_{{\rm{c(h)}}})$ represents the relative deviation to the classical value of the mean energy $N/{\beta }_{{\rm{c(h)}}}$ for an equilibrium state in the canonical ensemble, see appendix B.1. Using (4) for a sudden quench driving we explore numerically the optimization of the output power $P=-\langle W\rangle /\tau$. It is clear that for large particle number, the efficiency at optimal output power deviates from the Rezek–Kosloff efficiency given by the leading term in equation (9). In this scenario, nonadiabatic effects during sudden-quench driving can be used to achieve quantum supremacy. In particular, they lead to a many-particle enhancement of the efficiency at optimal output power by up to ∼50% of the single-particle value (for ${\sigma }_{{\rm{c}}}\geqslant 10$, typically), see appendix E for a derivation of this upper bound. In figure 3 we observe that for $\lambda =0$ (free bosons) and for $N\leqslant 300$, the optimal output power (and efficiency at optimal output power) of a many-particle QHE surpasses the value of $N$ independent single-particle QHEs. For $N\geqslant 300$ the relative efficiency at optimal output power increases with $N$ while the many-particle optimal output power can be improved (with respect to $N$ single-particle QHEs) for ${\beta }_{{\rm{h}}}/{\beta }_{{\rm{c}}}$ less than a critical value (∼0.2 for $N=500$). For a weakly interacting Bose gas (typically $\lambda \leqslant 0.2$) the optimal output power and efficiency are boosted for a range of $N$ below a critical value depending on the interaction strength λ and also on temperatures ${\beta }_{{\rm{h}}}$ and ${\beta }_{{\rm{c}}}$.

We consider a confined quantum fluid as a working medium with Hamiltonian

$$\begin{eqnarray}&&H=\displaystyle \sum _{i=1}^{N}\left[-\displaystyle \frac{{{\hslash }}^{2}}{2m}\displaystyle \frac{{\partial }^{2}}{\partial {z}_{i}^{2}}+\displaystyle \frac{1}{2}m\omega {(t)}^{2}{z}_{i}^{2}\right]+\displaystyle \sum _{i\lt j}V({z}_{i}-{z}_{j}),\end{eqnarray} \tag{ A1 }$$

where the inter-particle interactions are described by a potential that exhibits the scaling property $V(\lambda z)={\lambda }^{-2}V(z)$. The dynamical group for the Hamiltonian (A1) is ${SU}(1,1)$ [22]. For simplicity, in this appendix we consider a one-dimensional system but one would obtain the same results for higher dimension. Remarkably, for the broad family of systems described by Hamiltonian (A1) the exact quantum dynamics of an arbitrary eigenstate under a modulation of the trapping frequency $\omega (t)$ is described by a scaling law according to which the time-evolving state can be written in terms of the state at t = 0 [23, 24]

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\{{m}_{i}\}}({z}_{1},\ldots ,{z}_{N},t)={b}^{-\tfrac{N}{2}}\exp \left[{\rm{i}}\displaystyle \frac{m\dot{b}}{2{\hslash }b}\displaystyle \sum _{i=1}^{N}{z}_{i}^{2}-{\rm{i}}{\int }_{0}^{t}\displaystyle \frac{E(\{{m}_{i}\})}{{\hslash }b{(t^{\prime} )}^{2}}{\rm{d}}t^{\prime} \right]{{\rm{\Psi }}}_{\{{m}_{i}\}}\left(\displaystyle \frac{{z}_{1}}{b},\ldots ,\displaystyle \frac{{z}_{N}}{b},t=0\right),\end{eqnarray} \tag{ A2 }$$

where $E(\{{m}_{i}\})$ is the eigenvalue of (A1) associated with the multi-index $\{{m}_{i}| i=1,\ldots ,N\}$ defining a complete set of quantum numbers. Here, $b=b(t)\gt 0$ is the scaling factor that can be obtained as a solution of the Ermakov differential equation

$$\begin{eqnarray}&&\ddot{b}+\omega {(t)}^{2}b={\omega }_{0}^{2}{b}^{-3}.\end{eqnarray} \tag{ A3 }$$

The boundary conditions $b(0)=1$ and $\dot{b}(0)=0$ follow from the requirement for the scaling law to reduce to the initial eigenstate at t = 0 at the beginning of the process. For the Calogero–Sutherland model, equation (A2) generalizes the scaling law previously reported for the ground state (m_i = 0 for all $i=1,\cdots ,N)$ in [38, 39] but follows from the exact coherent states under scale-invariant driving found for a broad family of many-body systems in [24, 40].

In what follows, we shall give a general derivation of the relation between the nonadiabatic and the adiabatic mean energies given before the equation (3) in the main body of the article, and that holds for all quantum fluids described by Hamiltonian (A1). The scaling dynamics (A2) is associated with the ${SU}(1,1)$ dynamical symmetry group [22], shared by the familiar single-particle time-dependent harmonic oscillator [41]. We follow Lohe [41] and describe the time-evolution in terms of the action of two spatial unitary transformations, elements of ${SU}(1,1)$. To this end, we introduce

$$\begin{eqnarray}&&{T}_{{\rm{dil}}}=\exp \left[-{\rm{i}}\displaystyle \frac{\mathrm{log}b}{2{\hslash }}\displaystyle \sum _{i=1}^{N}({z}_{i}{p}_{i}+{p}_{i}{z}_{i})\right],\end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray}&&{T}_{z}\,=\,\exp \left[{\rm{i}}\displaystyle \frac{m\dot{b}}{2{\hslash }b}\displaystyle \sum _{i=1}^{N}{z}_{i}^{2}\right],\end{eqnarray} \tag{ A5 }$$

to derive the invariant of motion

$$\begin{eqnarray}&&{\text{I}}={{TH}}_{0}{T}^{\dagger }=\displaystyle \sum _{i=1}^{N}\left[\displaystyle \frac{{b}^{2}{p}_{i}^{2}}{2m}+\displaystyle \frac{1}{2}{{mz}}_{i}^{2}\left(\displaystyle \frac{\omega {(0)}^{2}}{{b}^{2}}+{\dot{b}}^{2}\right)-\displaystyle \frac{b\dot{b}}{2}\ ({z}_{i}{p}_{i}+{p}_{i}{z}_{i})\right]+{b}^{2}\displaystyle \sum _{i\lt j}V({z}_{i}-{z}_{j}),\end{eqnarray} \tag{ A6 }$$

where $T={T}_{z}{T}_{{\rm{dil}}}$, and where H₀ denotes the Hamiltonian (A1) at t = 0. This invariant is an instance of the many-body invariants known under scale-invariant driving, see equation (56) in [42]. It has time-dependent eigenfunctions ${{\rm{\Phi }}}_{\{{m}_{i}\}}({z}_{1},\ldots ,{z}_{N},t)=T{{\rm{\Psi }}}_{\{{m}_{i}\}}({z}_{1},\ldots ,{z}_{N},0)$ and time-independent eigenvalues $E(\{{m}_{k}\})$. It also means that the time-evolving state (A2) can be written as

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\{{m}_{i}\}}({z}_{1},\ldots ,{z}_{N},t)={{\rm{e}}}^{{\rm{i}}{\alpha }_{t}}T{{\rm{\Psi }}}_{\{{m}_{i}\}}({z}_{1},\ldots ,{z}_{N},0),\end{eqnarray} \tag{ A7 }$$

where ${\alpha }_{t}={\int }_{0}^{t}\tfrac{E(\{{m}_{i}\})}{{\hslash }b{(t^{\prime} )}^{2}}{\rm{d}}t^{\prime}$ is a time-dependent phase. As a function of it, Hamiltonian (A1) can be rewritten as

$$\begin{eqnarray}&&H=\displaystyle \frac{1}{{b}^{2}}{\text{I}}\,+{\rm{i}}{\hslash }\displaystyle \frac{\partial T}{\partial t}{T}^{\dagger },\end{eqnarray} \tag{ A8 }$$

where

$$\begin{eqnarray}&&{\rm{i}}{\hslash }\displaystyle \frac{\partial T}{\partial t}{T}^{\dagger }=\displaystyle \sum _{i=1}^{N}\left[-\displaystyle \frac{1}{2}{{mz}}_{i}^{2}\left(\displaystyle \frac{\ddot{b}}{b}+\displaystyle \frac{{\dot{b}}^{2}}{{b}^{2}}\right)+\displaystyle \frac{\dot{b}}{2b}\ ({z}_{i}{p}_{i}+{p}_{i}{z}_{i})\right].\end{eqnarray} \tag{ A9 }$$

By (A7) we find that

$$\begin{eqnarray}&&\langle {{\rm{\Psi }}}_{\{{m}_{i}\}}(t)| H| {{\rm{\Psi }}}_{\{{m}_{i}\}}(t)\rangle =\langle {{\rm{\Psi }}}_{\{{m}_{i}\}}(0)| {T}^{\dagger }{HT}| {{\rm{\Psi }}}_{\{{m}_{i}\}}(0)\rangle =\langle {{\rm{\Psi }}}_{\{{m}_{i}\}}(0)| \left(\displaystyle \frac{{H}_{0}}{{b}^{2}}+{\rm{i}}{\hslash }{T}^{\dagger }\displaystyle \frac{\partial T}{\partial t}\right)| {{\rm{\Psi }}}_{\{{m}_{i}\}}(0)\rangle ,\end{eqnarray} \tag{ A10 }$$

where we use

$$\begin{eqnarray}&&{\rm{i}}{\hslash }{T}^{\dagger }\displaystyle \frac{\partial T}{\partial t}=\displaystyle \sum _{i=1}^{N}\left[\displaystyle \frac{1}{2}({\dot{b}}^{2}-\ddot{b}b){z}_{i}^{2}+\displaystyle \frac{\dot{b}}{2b}({z}_{i}{p}_{i}+{p}_{i}{z}_{i})\right],\end{eqnarray} \tag{ A11 }$$

which is derived from equation (A9). As a result, the mean nonadiabatic energy is given by

$$\begin{eqnarray}\langle {{\rm{\Psi }}}_{\{{m}_{i}\}}(t)| H| {{\rm{\Psi }}}_{\{{m}_{i}\}}(t)\rangle & = & \displaystyle \frac{E(\{{m}_{i}\})}{{b}^{2}}-\displaystyle \frac{m}{2}(b\ddot{b}-{\dot{b}}^{2})\langle {{\rm{\Psi }}}_{\{{m}_{i}\}}(0)| \displaystyle \sum _{i=1}^{N}{z}_{i}^{2}| {{\rm{\Psi }}}_{\{{m}_{i}\}}(0)\rangle \\ & & +\displaystyle \frac{\dot{b}}{2b}\langle {{\rm{\Psi }}}_{\{{m}_{i}\}}(0)| \displaystyle \sum _{i=1}^{N}({z}_{i}{p}_{i}+{p}_{i}{z}_{i})| {{\rm{\Psi }}}_{\{{m}_{i}\}}(0)\rangle .\end{eqnarray} \tag{ A12 }$$

During time-evolution, the mean-energy of an initial thermal state

$$\begin{eqnarray}&&\rho =\displaystyle \sum _{\{{m}_{i}\}}{p}_{\{{m}_{i}\}}| {{\rm{\Psi }}}_{\{{m}_{i}\}}(0)\rangle \langle {{\rm{\Psi }}}_{\{{m}_{i}\}}(0)| ,\end{eqnarray} \tag{ A13 }$$

with ${p}_{\{{m}_{i}\}}\geqslant 0$ and ${\sum }_{\{{m}_{i}\}}{p}_{\{{m}_{i}\}}=1$ is simply given by

$$\begin{eqnarray}&&\langle H(t)\rangle =\displaystyle \sum _{\{{m}_{i}\}}{p}_{\{{m}_{i}\}}\langle {{\rm{\Psi }}}_{\{{m}_{i}\}}(t)| H| {{\rm{\Psi }}}_{\{{m}_{i}\}}(t)\rangle .\end{eqnarray} \tag{ A14 }$$

Collecting all the previous equations holding for the general Hamiltonian (A1), the nonadiabatic mean-energy following a change of $\omega (t)$ is found to be

$$\begin{eqnarray}&&\langle H(t)\rangle =\displaystyle \frac{1}{{b}^{2}}\langle H(0)\rangle +\displaystyle \frac{\dot{b}}{2b}\displaystyle \sum _{i=1}^{N}\langle \{{z}_{i},{p}_{i}\}(0)\rangle +\displaystyle \frac{m}{2}({\dot{b}}^{2}-b\ddot{b})\displaystyle \sum _{i=1}^{N}\langle {z}_{i}^{2}(0)\rangle ,\end{eqnarray} \tag{ A15 }$$

where p_i is the momentum of the ith particle and $\{{z}_{i},{p}_{i}\}={z}_{i}{p}_{i}+{p}_{i}{z}_{i}$. As a result of the underlying dynamical symmetry, it suffices to know the evolution of the scaling factor and the initial expectation values at t = 0 of the mean energy, squeezing and position dispersion to characterize the nonadiabatic dynamics. Thanks to the thermal equilibrium property of the initial state, we can compute explicitly these quantities for the general Hamiltonian (A1).

First, it is not difficult to show that at thermal equilibrium the mean squeezing vanishes

$$\begin{eqnarray}&&\langle \{{z}_{i},{p}_{i}\}(0)\rangle =0.\end{eqnarray} \tag{ A16 }$$

Indeed, it suffices to rewrite the squeezing operator as

$$\begin{eqnarray*}&&2{\rm{i}}{\hslash }\{{z}_{i},{p}_{i}\}=[{z}_{i}^{2},{p}_{i}^{2}]=[{z}_{i}^{2},H],\end{eqnarray*}$$

where we used the general form (A1). From (A14) it is clear that

$$\begin{eqnarray*}&&\langle [{z}_{i}^{2},H](0)\rangle =0,\end{eqnarray*}$$

which holds for any density matrix diagonal in the eigenbasis of H, and in particular, for the thermal density matrix ρ (A13).

Secondly, one can show that the initial mean position dispersion is proportional to the mean energy of the system

$$\begin{eqnarray}&&\displaystyle \sum _{i}\langle {z}_{i}^{2}(0)\rangle =\displaystyle \frac{\langle H(0)\rangle }{m{\omega }_{0}^{2}}.\end{eqnarray} \tag{ A17 }$$

Again, this formula is very general and holds for an initial thermal state. To derive (A17) we use the Hellmann–Feynman theorem ${\sum }_{i}\langle {z}_{i}^{2}(0)\rangle ={(m{\omega }_{0})}^{-1}\partial \langle H(0)\rangle /\partial {\omega }_{0}$. Using the scaling invariance of the potential V, we have $V(\sqrt{{\omega }_{0}}z)={\omega }_{0}^{-1}V(z)$. It follows that $H(0)={\omega }_{0}\tilde{H}(0)$ with $\tilde{H}(0)={\sum }_{i=1}^{N}\left(-\tfrac{{{\hslash }}^{2}}{2m}\tfrac{{\partial }^{2}}{\partial {q}^{2}}+\tfrac{m}{2}{q}^{2}\right)+{\sum }_{i\lt j}V({q}_{i}-{q}_{j})$ where we rescaled the position as ${q}_{i}=\sqrt{{\omega }_{0}}{z}_{i}$. Therefore, we find that $\partial \langle H(0)\rangle /\partial {\omega }_{0}={\omega }_{0}^{-1}\langle H(0)\rangle$.

Gathering equations (A15), (A16), and (A17), we obtain

$$\begin{eqnarray}&&\langle H(t)\rangle =\displaystyle \frac{{Q}^{* }(t)}{{b}_{{\rm{ad}}}^{2}}\langle H(0)\rangle ,\end{eqnarray} \tag{ A18 }$$

where ${b}_{{\rm{ad}}}={[{\omega }_{0}/\omega (t)]}^{1/2}$ and where the nonadiabatic factor ${Q}^{* }(t)$ that accounts for the amount of energy excitations over the adiabatic dynamics is given by

$$\begin{eqnarray}&&{Q}^{* }(t)={b}_{{\rm{ad}}}^{2}\left[\displaystyle \frac{1}{2b{(t)}^{2}}+\displaystyle \frac{\omega {(t)}^{2}b{(t)}^{2}}{2{\omega }_{0}^{2}}+\displaystyle \frac{\dot{b}{(t)}^{2}}{2{\omega }_{0}^{2}}\right].\end{eqnarray} \tag{ A19 }$$

Equations (A18) and (A19) give the general scaling law of the nonadiabatic mean energy of a class of systems governed by Hamiltonian (A1). We further recognize

$$\begin{eqnarray}&&{\langle H(t)\rangle }_{{\rm{ad}}{\rm{}}}=\displaystyle \frac{1}{{b}_{{\rm{ad}}}^{2}}\langle H(0)\rangle \,\end{eqnarray} \tag{ A20 }$$

as the mean energy in the adiabatic limit. As a result, ${Q}^{* }(t)$ is simply the ratio between the nonadiabatic and adiabatic mean energies.

In this section we prove that the nonadiabatic factor ${Q}^{* }(t)$ introduced in equation (A19) is equal to that introduced by Husimi for the single-particle time-dependent harmonic oscillator ${Q}_{{\rm{H}}}^{* }(t)$ [25],

$$\begin{eqnarray}&&{Q}_{{\rm{H}}}^{* }(t)={Q}^{* }(t).\end{eqnarray} \tag{ A21 }$$

As discussed, the nonadiabatic factor ${Q}^{* }(t)$ provides the ratio between the nonadiabatic and adiabatic mean energies for an infinite family of harmonically confined quantum fluids following a variation of the trapping frequency [24], whenever the dynamics is scale-invariant and the initial state at t = 0 is a thermal state (with vanishing expectation value of the squeezing operator). The nonadiabatic factor introduced by Husimi is defined as

$$\begin{eqnarray}&&{Q}_{{\rm{H}}}^{* }(t)=\displaystyle \frac{({\dot{{G}_{1}}}^{2}+\omega {(t)}^{2}{G}_{1}{(t)}^{2})+{\omega }_{0}^{2}({\dot{{G}_{2}}}^{2}+\omega {(t)}^{2}{G}_{2}{(t)}^{2})}{2{\omega }_{0}\omega (t)},\end{eqnarray} \tag{ A22 }$$

where ${G}_{1}(t)$ and ${G}_{2}(t)$ are two fundamental solution of the classical harmonic oscillator equation $\ddot{x}(t)+\omega {(t)}^{2}x(t)=0$ that satisfy the following initial conditions:

$$\begin{eqnarray}&&{G}_{1}(0)=1,\dot{{G}_{1}}(0)=0,\end{eqnarray} \tag{ A23a }$$

$$\begin{eqnarray}&&{G}_{2}(0)=0,\dot{{G}_{1}}(0)=1.\end{eqnarray} \tag{ A23b }$$

The Wronskian of these two functions ${W}_{{G}_{1},{G}_{2}}(t)=\dot{{G}_{1}}(t){G}_{2}(t)-{G}_{1}(t)\dot{{G}_{2}}(t)$ is shown to be constant

$$\begin{eqnarray}&&{W}_{{G}_{1},{G}_{2}}(t)={W}_{{G}_{1},{G}_{2}}(0)=-1.\end{eqnarray} \tag{ A24 }$$

As a preliminary result, we recall the derivation of a simple and useful formula previously derived in [43] for the scaling factor b(t) solution of the Ermarkov equation (A3)

$$\begin{eqnarray}&&b(t)=\sqrt{{G}_{1}{(t)}^{2}+{\omega }_{0}^{2}{G}_{2}{(t)}^{2}}.\end{eqnarray} \tag{ A25 }$$

To prove (A25), we first compute the first time-derivative of b(t)

$$\begin{eqnarray}&&\dot{b}(t)=\displaystyle \frac{\dot{{G}_{1}}(t){G}_{1}(t)+{\omega }_{0}^{2}\dot{{G}_{2}}(t){G}_{2}(t)}{\sqrt{{G}_{1}{(t)}^{2}+{\omega }_{0}^{2}{G}_{2}{(t)}^{2}}}=\displaystyle \frac{\dot{{G}_{1}}(t){G}_{1}(t)+{\omega }_{0}^{2}\dot{{G}_{2}}(t){G}_{2}(t)}{b(t)},\end{eqnarray} \tag{ A26 }$$

as well as its second derivative

$$\begin{eqnarray*}\ddot{b}(t) & = & \displaystyle \frac{\ddot{{G}_{1}}(t){G}_{1}(t)+{\omega }_{0}^{2}\ddot{{G}_{2}}(t){G}_{2}(t)}{b(t)}+\displaystyle \frac{\dot{{G}_{1}}{(t)}^{2}+{\omega }_{0}^{2}\dot{{G}_{2}}{(t)}^{2}}{b(t)}-\dot{b}(t)\displaystyle \frac{\dot{{G}_{1}}(t){G}_{1}(t)+{\omega }_{0}^{2}\dot{{G}_{2}}(t){G}_{2}(t)}{b{(t)}^{2}},\\ & = & -\omega {(t)}^{2}\displaystyle \frac{{G}_{1}{(t)}^{2}+{\omega }_{0}^{2}{G}_{2}^{2}}{b(t)}+\displaystyle \frac{(\dot{{G}_{1}}{(t)}^{2}+{\omega }_{0}\dot{{G}_{2}}{(t)}^{2})({G}_{1}{(t)}^{2}+{\omega }_{0}^{2}{G}_{2}^{2})}{b{(t)}^{3}}\\ & & -\displaystyle \frac{{(\dot{{G}_{1}}(t){G}_{1}(t)+{\omega }_{0}^{2}\dot{{G}_{2}}(t){G}_{2}(t))}^{2}}{b{(t)}^{3}}.\end{eqnarray*}$$

The first term of the latter equation reduces to

$$\begin{eqnarray}&&-\omega {(t)}^{2}b(t).\end{eqnarray} \tag{ A27 }$$

The second term is given by

$$\begin{eqnarray*}&&\displaystyle \frac{\dot{{G}_{1}}{(t)}^{2}{G}_{1}{(t)}^{2}+{\omega }_{0}^{2}\dot{{G}_{1}}{(t)}^{2}{G}_{2}{(t)}^{2}+{\omega }_{0}^{2}\dot{{G}_{2}}{(t)}^{2}{G}_{1}{(t)}^{2}+{\omega }_{0}^{4}\dot{{G}_{2}}{(t)}^{2}{G}_{2}{(t)}^{2}}{b{(t)}^{3}},\end{eqnarray*}$$

while the third term gives

$$\begin{eqnarray*}&&-\displaystyle \frac{\dot{{G}_{1}}{(t)}^{2}{G}_{1}{(t)}^{2}+2{\omega }_{0}^{2}\dot{{G}_{1}}(t){G}_{1}(t)\dot{{G}_{2}}(t){G}_{2}(t)+{\omega }_{0}^{4}\dot{{G}_{2}}{(t)}^{2}{G}_{2}{(t)}^{2}}{b{(t)}^{3}}.\end{eqnarray*}$$

Thus, after adding these two terms, we find

$$\begin{eqnarray}&&{\omega }_{0}^{2}\displaystyle \frac{{(\dot{{G}_{1}}(t){G}_{2}(t)-\dot{{G}_{2}}(t){G}_{1}(t))}^{2}}{b{(t)}^{3}}={\omega }_{0}^{2}\displaystyle \frac{{({W}_{{G}_{1},{G}_{2}}(t))}^{2}}{b{(t)}^{3}}=\displaystyle \frac{{\omega }_{0}^{2}}{b{(t)}^{3}},\end{eqnarray} \tag{ A28 }$$

where we used the property (A24). Gathering (A27) and (A28) one gets the Ermarkov equation. One can easily check the consistency of the initial conditions using (A23), (A25) and (A26).

To prove the identity (A21), using (A25) we rewrite (A22) and (A19) as

$$\begin{eqnarray*}{Q}_{{\rm{H}}}^{* }(t) & = & \displaystyle \frac{\omega {(t)}^{2}{G}_{1}{(t)}^{2}+{\omega }_{0}^{2}\omega {(t)}^{2}{G}_{2}{(t)}^{2}}{2{\omega }_{0}\omega (t)}+\displaystyle \frac{{\dot{{G}_{1}}}^{2}+{\omega }_{0}^{2}{\dot{{G}_{2}}}^{2}}{2{\omega }_{0}\omega (t)}=\displaystyle \frac{\omega (t)b(t)}{2{\omega }_{0}}+\displaystyle \frac{{\dot{{G}_{1}}}^{2}+{\omega }_{0}^{2}{\dot{{G}_{2}}}^{2}}{2{\omega }_{0}\omega (t)},\\ {Q}^{* }(t) & = & \displaystyle \frac{\omega (t){b}^{2}}{{\omega }_{0}}+\displaystyle \frac{1}{2{\omega }_{0}\omega (t){b}^{2}}({\omega }_{0}^{2}+{\dot{b}}^{2}{b}^{2}),\end{eqnarray*}$$

From (A26) and (A24), we have

$$\begin{eqnarray*}{\dot{b}}^{2}{b}^{2} & = & {(\dot{{G}_{1}}{G}_{1}+{\omega }_{0}^{2}\dot{{G}_{2}}{G}_{2})}^{2}={\dot{{G}_{1}}}^{2}{G}_{1}^{2}+2{\omega }_{0}^{2}\dot{{G}_{1}}{G}_{1}\dot{{G}_{2}}{G}_{2}+{\omega }_{0}^{4}{\dot{{G}_{2}}}^{2}{G}_{2}^{2},\\ {\omega }_{0}^{2} & = & {\omega }_{0}^{2}{({W}_{{G}_{1},{G}_{2}}(t))}^{2}={\omega }_{0}^{2}({\dot{{G}_{1}}}^{2}{G}_{2}^{2}+{\dot{{G}_{2}}}^{2}{G}_{1}^{2}-2\dot{{G}_{1}}{G}_{1}\dot{{G}_{2}}{G}_{2}).\end{eqnarray*}$$

Thus,

$$\begin{eqnarray*}&&{\omega }_{0}^{2}+{\dot{b}}^{2}{b}^{2}={G}_{1}^{2}({\dot{{G}_{1}}}^{2}+{\omega }_{0}^{2}{G}_{2}^{2})+{\omega }_{0}^{2}{\dot{{G}_{2}}}^{2}({\dot{{G}_{1}}}^{2}+{\omega }_{0}^{2}{G}_{2}^{2})={b}^{2}({\dot{{G}_{1}}}^{2}+{\omega }_{0}^{2}{\dot{{G}_{2}}}^{2}).\end{eqnarray*}$$

Gathering the last equalities we obtain the identity (A21).

Although the Calogero–Sutherland model described by the Hamiltonian in equation (7) is a truly interacting many-body system, to account for its thermodynamics it is useful to exploit its mapping to effectively free harmonic oscillators [28, 35, 44, 45]. Instead of using the set of quantum numbers $\{{m}_{i}| i=1,\ldots ,N\}$, it is convenient to use the set of occupation numbers $\{{n}_{j}| j=1,\ldots ,\infty \}$ in each renormalized harmonic j-mode of energy ${\varepsilon }_{j}=j{\hslash }\omega$, with the identification ${n}_{j}={\sum }_{i=1}^{N}{\delta }_{j,{m}_{i}}$ where ${\delta }_{l,m}$ is the Kronecker delta. As shown in the main body of the article, the canonical partition function of the (bosonic) Calogero–Sutherland gas with interaction strength λ is

$$\begin{eqnarray}&&{Z}_{N}^{(\lambda )}=\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{\{{n}_{j}\}=1}{{\displaystyle \sum }_{j=1}^{\infty }{n}_{j}=N}}^{\infty }{{\rm{e}}}^{-\beta E(\{{n}_{j}\})}=\displaystyle \prod _{k=1}^{N}\displaystyle \frac{{{\rm{e}}}^{-\beta \tfrac{{\hslash }\omega }{2}[1+\lambda (N-1)]}}{1-{{\rm{e}}}^{-\beta k{\hslash }\omega }},\end{eqnarray} \tag{ B1 }$$

where $\{{n}_{j}\}$ is the occupation number in the jth mode of energy $j{\hslash }\omega$, satisfying the normalization ${\sum }_{j=1}^{\infty }{n}_{k}=N$, and $E(\{{n}_{j}\})=\tfrac{{\hslash }\omega }{2}N[1+\lambda (N-1)]+{\sum }_{j=1}^{\infty }j{\hslash }\omega {n}_{j}$.

It follows that the equilibrium mean energy, for a given temperature $(\omega ,\beta )$, is given by

$$\begin{eqnarray}&&{{ \mathcal E }}_{N}(\omega ,\beta )=\displaystyle \frac{{\hslash }\omega N}{2}[1+\lambda (N-1)]+\displaystyle \sum _{k=1}^{N}\displaystyle \frac{k{\hslash }\omega }{{{\rm{e}}}^{k{\hslash }\beta \omega }-1},\end{eqnarray} \tag{ B2 }$$

which can be read as

$$\begin{eqnarray}&&{{ \mathcal E }}_{N}(\omega ,\beta )={{ \mathcal E }}^{(0)}(\omega )+{{ \mathcal E }}^{(\mathrm{th})}(\omega ,\beta ),\end{eqnarray} \tag{ B3 }$$

where the ground state energy ${{ \mathcal E }}_{N}^{(0)}(\omega )$ and the thermal energy ${{ \mathcal E }}_{N}^{(\mathrm{th})}(\omega ,\beta )$ correspond to the first and second term of (B2), respectively. Depending on the value of $\sigma =N{\hslash }\beta \omega$, three different regimes can be distinguished:

• At very low temperature $\sigma \gg N$, the mean energy is approximately given by the ground state energy
  $$\begin{eqnarray}&&{{ \mathcal E }}_{N}(\omega ,\beta )\approx {{ \mathcal E }}_{N}^{(0)}(\omega )=\displaystyle \frac{N{\hslash }\omega }{2}[1+\lambda (N-1)].\end{eqnarray} \tag{ B4 }$$
• At high temperature $\sigma \ll 1$, we find corrections to the thermal energy due to quantum fluctuations
  $$\begin{eqnarray}&&{{ \mathcal E }}_{N}(\omega ,\beta )\approx \displaystyle \frac{N}{\beta }\left[1+\displaystyle \frac{\sigma }{2}{g}_{N}\left(\lambda -\displaystyle \frac{1}{2}\right)+{\sigma }^{2}\displaystyle \frac{(N+1)(2N+1)}{72{N}^{2}}\right],\end{eqnarray} \tag{ B5 }$$
  where ${g}_{N}=(N-1){N}^{-1}$ and we have disregarded corrections of order $O({\sigma }^{3})$. Formula (B5) is derived using the fact that for $x\lt 1$, ${({{\rm{e}}}^{x}-1)}^{-1}\approx \tfrac{1}{x}-\tfrac{1}{2}+\tfrac{x}{12}$ so that
  $$\begin{eqnarray*}&&{{ \mathcal E }}_{N}(\omega ,\beta )={{ \mathcal E }}_{N}^{(0)}(\omega )+\displaystyle \frac{1}{\beta }\displaystyle \sum _{k=1}^{N}\left(1-\displaystyle \frac{k\sigma }{2N}+\displaystyle \frac{{k}^{2}{\sigma }^{2}}{12{N}^{2}}\right).\end{eqnarray*}$$
  The sums can be carried out explicitly ${\sum }_{k=1}^{N}k=N(N+1)/2$ and ${\sum }_{k=1}^{N}{k}^{2}=N(N+1)(2N+1)/6$.
• For many-particles, a novel regime absent in the single particle case is found for intermediate temperatures such that $1\lt \sigma \ll N$ (which requires $N$ to be large compared to 1)
  $$\begin{eqnarray}&&{{ \mathcal E }}_{N}(\omega ,\beta )\approx {{ \mathcal E }}_{N}^{(0)}(\omega )+\displaystyle \frac{1}{{\hslash }\omega {\beta }^{2}}{\int }_{0}^{\sigma }{\rm{d}}s\displaystyle \frac{s}{{{\rm{e}}}^{s}-1}\approx \displaystyle \frac{N}{\beta }{\mu }_{\lambda }(\sigma ),\end{eqnarray} \tag{ B6 }$$
  where
  $$\begin{eqnarray}&&{\mu }_{\lambda }(\sigma )=\displaystyle \frac{1}{\sigma }{\int }_{0}^{\sigma }{\rm{d}}s\displaystyle \frac{s}{{{\rm{e}}}^{s}-1}+\sigma \displaystyle \frac{\lambda }{2}.\end{eqnarray} \tag{ B7 }$$
  To obtain the last equation, we estimated the Riemann sum in (B2) by a Riemann integral since ${\hslash }\beta \omega \ll 1$ (the corrections are bounded by $N{({\hslash }\beta \omega )}^{2}/2$ and so $N$ can not be too large). The integral in (B7) can be found in closed form
  $$\begin{eqnarray*}&&{\int }_{0}^{\sigma }{\rm{d}}s\displaystyle \frac{s}{{{\rm{e}}}^{s}-1}=\displaystyle \frac{{\pi }^{2}}{6}+\sigma \mathrm{log}(1-{{\rm{e}}}^{-\sigma })-{\mathrm{Li}}_{2}({{\rm{e}}}^{-\sigma }),\end{eqnarray*}$$
  in terms of the standard polylogarithm function ${\mathrm{Li}}_{n}(z)={\sum }_{j=1}^{\infty }{j}^{-n}{z}^{j},n\geqslant 1$ [46]. The behavior of this function as a function of σ for different values of λ is shown in figure B1 . We will see in the following sections that this function plays an important role in the optimization of the total work.
• The very high temperature limit, $N{\hslash }\beta \omega \ll 1$, corresponds to the classical limit (quantum fluctuation are negligible in this regime) where the mean energy reads
  $$\begin{eqnarray}&&{{ \mathcal E }}_{N}(\omega ,\beta )\approx \displaystyle \frac{N}{\beta }.\end{eqnarray} \tag{ B8 }$$
  We note that to derive (B8) we use formula (B5) where we neglect the quantum corrections $O({\hslash }\beta \omega )$.

Zoom In Zoom Out Reset image size

Consider now a given protocol $\omega (t)$ for the compression stroke such that $\omega (0)={\omega }_{1},\omega (\tau )={\omega }_{2}$ (similarly, for the compression stroke, let $\widetilde{\omega }(t)$ satisfy $\widetilde{\omega }(0)={\omega }_{2},\widetilde{\omega }(\tau )={\omega }_{1}$). We denote the nonadiabatic factor ${Q}_{{AB}}^{* }(\tau )$ (respectively ${Q}_{{CD}}^{* }(\tau )$) corresponding to the step 1 of the Otto cycle $A\to B$ (respectively the step 3 $C\to D$), see formula (A22). Let us assume that for any $\tau \gt 0$ the nonadiabatic factor is continuously differentiable with respect to ${\omega }_{1}$ and ${\omega }_{2}$. We recall that for the adiabatic case, ${Q}_{{AB}}^{* }={Q}_{{CD}}^{* }=1$ while in the nonadiabatic sudden-quench limit ${Q}_{{AB}}^{* }={Q}_{{CD}}^{* }=\tfrac{{\omega }_{1}^{2}+{\omega }_{2}^{2}}{2{\omega }_{1}{\omega }_{2}}=\tfrac{{x}^{2}+1}{2x}$ with $x={\omega }_{1}/{\omega }_{2}$.

For the Otto cycle, the total work is usually defined as $W=\langle {W}_{1}\rangle +\langle {W}_{3}\rangle$, where $\langle {W}_{1}\rangle \gt 0$ is the work corresponding to the compression (${\omega }_{1}\to {\omega }_{2}\gt {\omega }_{1}$) and $\langle {W}_{3}\rangle \lt 0$ is the work corresponding to the expansion (${\omega }_{2}\to {\omega }_{1}\lt {\omega }_{2}$). The engine produces work when W is negative. Alternatively, we define

$$\begin{eqnarray*}&&\tilde{W}=-(\langle {W}_{1}\rangle +\langle {W}_{3}\rangle ).\end{eqnarray*}$$

We recall the equations given in the main body of the article:

$$\begin{eqnarray}&&\langle {W}_{1}\rangle ={\langle H\rangle }_{B}-{\langle H\rangle }_{A}=\left({Q}_{{AB}}^{* }\displaystyle \frac{{\omega }_{2}}{{\omega }_{1}}-1\right){\langle H\rangle }_{A},\end{eqnarray} \tag{ B9a }$$

$$\begin{eqnarray}&&\langle {W}_{3}\rangle ={\langle H\rangle }_{D}-{\langle H\rangle }_{C}=\left({Q}_{{CD}}^{* }\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}}-1\right){\langle H\rangle }_{C},\end{eqnarray} \tag{ B9b }$$

$$\begin{eqnarray}&&\langle {Q}_{2}\rangle ={\langle H\rangle }_{C}-{\langle H\rangle }_{B}={\langle H\rangle }_{C}-{Q}_{{AB}}^{* }\displaystyle \frac{{\omega }_{2}}{{\omega }_{1}}{\langle H\rangle }_{A},\end{eqnarray} \tag{ B9c }$$

$$\begin{eqnarray}&&\langle {Q}_{4}\rangle ={\langle H\rangle }_{D}-{\langle H\rangle }_{A}={Q}_{{CD}}^{* }\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}}{\langle H\rangle }_{C}-{\langle H\rangle }_{A},\end{eqnarray} \tag{ B9d }$$

where ${\langle H\rangle }_{A}={{ \mathcal E }}_{N}({\omega }_{1},{\beta }_{{\rm{c}}})$ and ${\langle H\rangle }_{C}={{ \mathcal E }}_{N}({\omega }_{2},{\beta }_{{\rm{h}}})$ denote the energies at the equilibrium.

As a consequence, the total work of the engine is given by

$$\begin{eqnarray}&&\tilde{W}=\left(1-{Q}_{{AB}}^{* }\displaystyle \frac{{\omega }_{2}}{{\omega }_{1}}\right){\langle H\rangle }_{A}+\left(1-{Q}_{{CD}}^{* }\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}}\right){\langle H\rangle }_{C},\end{eqnarray} \tag{ B10 }$$

which is also also equal to $\tilde{W}=\langle {Q}_{2}\rangle +\left({\langle H\rangle }_{A}-{Q}_{{CD}}^{* }\tfrac{{\omega }_{1}}{{\omega }_{2}}{\langle H\rangle }_{C}\right)$. The efficiency can be written in the following form

$$\begin{eqnarray}&&\eta =\displaystyle \frac{\tilde{W}}{\langle {Q}_{2}\rangle }=1-\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}}\displaystyle \frac{{Q}_{{CD}}^{* }{\langle H\rangle }_{C}-\tfrac{{\omega }_{2}}{{\omega }_{1}}{\langle H\rangle }_{A}}{{\langle H\rangle }_{C}-{Q}_{{AB}}^{* }\tfrac{{\omega }_{2}}{{\omega }_{1}}{\langle H\rangle }_{A}}.\end{eqnarray} \tag{ B11 }$$

In the adiabatic case (${Q}^{* }=1$) the efficiency reduces to the Otto efficiency

$$\begin{eqnarray}&&{\eta }_{\mathrm{ad}}=1-\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}},\end{eqnarray} \tag{ B12 }$$

which is independent of $N$ and λ. However, for an adiabatic driving the total work (B10) depends on $N$ (and not on λ as we shall see in the section B.2). Thus, one can expect that the frequency ratio ${\omega }_{1}/{\omega }_{2}$ at which the work output is maximum depends on $N$. The efficiency at optimal work inherits a nontrivial dependence on $N$, as discussed in the section C.

For a finite-time protocol, the nonadiabatic factors ${Q}_{{AB}({CD})}^{* }\geqslant 1$ and the efficiency is bounded from above by the adiabatic Otto efficiency, which constitutes a universal limit as it does not depend on $N,\lambda$ or the specific driving protocol. This bound can be too conservative for finite-time thermodynamics. For example, for a sudden quench driving the efficiency can actually be upper bounded by half the adiabatic Otto efficiency. Hence, it is convenient to derive a tighter upper bound for the efficiency. To this end, let us rewrite the efficiency as follows

$$\begin{eqnarray*}&&\eta =1-{Q}_{{CD}}^{* }\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}}\displaystyle \frac{{\langle H\rangle }_{C}-{({Q}_{{CD}}^{* })}^{-1}\tfrac{{\omega }_{2}}{{\omega }_{1}}{\langle H\rangle }_{A}}{{\langle H\rangle }_{C}-{Q}_{{AB}}^{* }\tfrac{{\omega }_{2}}{{\omega }_{1}}{\langle H\rangle }_{A}}.\end{eqnarray*}$$

Since ${({Q}_{{CD}}^{* })}^{-1}\leqslant 1$ and ${Q}_{{AB}}^{* }\geqslant 1$, we find the following upper bound

$$\begin{eqnarray}&&\eta \leqslant 1-{Q}_{{CD}}^{* }\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}},\end{eqnarray} \tag{ B13 }$$

where ${Q}_{{CD}}^{* }\tfrac{{\omega }_{1}}{{\omega }_{2}}$ is nothing but $\tfrac{{\langle H\rangle }_{D}}{{\langle H\rangle }_{C}}$, see (B9d). Notice that this limit depends only on the factor ${Q}_{{CD}}^{* }$ and not on ${Q}_{{AB}}^{* }$. An upper bound depending on ${Q}_{{AB}}^{* }$ can as well be obtained, $\eta \leqslant 1-{({Q}_{{AB}}^{* })}^{-1}\tfrac{{\omega }_{1}}{{\omega }_{2}}=1-\tfrac{{\langle H\rangle }_{A}}{{\langle H\rangle }_{B}}$, but this is not interesting as it is weaker than the (adiabatic) Otto efficiency.

In what follows, we shall provide a derivation of the efficiency at optimal output power $P=\tilde{W}/\tau$, where τ is the running time of a cycle and $\tilde{W}$ is the absolute value of the output work per cycle. To this end, we first find the optimal frequency ratio corresponding to the maximal value of the total work per cycle using the exact expression given by equation (B10), and then we insert this optimal value in equation (B11) to obtain the efficiency at optimal output power.

In this section we discuss the effect of the interaction strength λ on the total work and the efficiency, for arbitrary dynamics. We consider a finite-time protocol for the Otto cycle as described in the previous section.

General statement:

• (i)  
  The total work (as well as the optimal work) and the efficiency are two monotonically decreasing functions of λ.
• (ii)  
  In the adiabatic limit the total work and the efficiency do not depend on λ.

We prove this important statement. We introduced λ as a superscript, explicitly. From the total work (B9a)–(B9d), one can easily show that

$$\begin{eqnarray*}{\tilde{W}}^{\lambda } & = & {\tilde{W}}^{0}+\displaystyle \frac{{\hslash }}{2}({\omega }_{1}-{Q}_{{AB}}^{* }{\omega }_{2}){C}_{N,\lambda }+\displaystyle \frac{{\hslash }}{2}({\omega }_{2}-{Q}_{{CD}}^{* }{\omega }_{1}){C}_{N,\lambda },\\ & = & {\tilde{W}}^{0}+\displaystyle \frac{{\hslash }}{2}[{\omega }_{1}-{\omega }_{2}-({Q}_{{AB}}^{* }-1){\omega }_{2}){C}_{N,\lambda }+\displaystyle \frac{{\hslash }}{2}({\omega }_{2}-{\omega }_{1}-({Q}_{{CD}}^{* }-1){\omega }_{1}]{C}_{N,\lambda },\\ & = & {\tilde{W}}^{0}+\displaystyle \frac{{\hslash }}{2}[(1-{Q}_{{AB}}^{* }){\omega }_{2}+(1-{Q}_{{CD}}^{* }){\omega }_{1}]{C}_{N,\lambda },\end{eqnarray*}$$

where ${C}_{N,\lambda }=N(N+1)\lambda$ and ${\tilde{W}}^{0}$ is the total work for $\lambda =0$. Knowing that ${Q}_{{AB}}^{* }\geqslant 1$ and ${Q}_{{CD}}^{* }\geqslant 1$, it follows that ${\partial }_{\lambda }{\tilde{W}}^{\lambda }\leqslant 0$. In addition, for the adiabatic case one finds ${\tilde{W}}^{\lambda }={\tilde{W}}^{0}$ (since ${Q}^{* }=1$) and the total work does not depends on λ. Similarly, one can show that ${\langle {Q}_{2}\rangle }^{\lambda }={\langle {Q}_{2}\rangle }^{0}+\tfrac{{\hslash }{\omega }_{2}}{2}(1-{Q}_{{AB}}^{* }){C}_{N,\lambda }$ so ${\partial }_{\lambda }{\langle {Q}_{2}\rangle }^{\lambda }\leqslant 0$. In regard to the efficiency:

$$\begin{eqnarray*}&&{\eta }^{\lambda }=\displaystyle \frac{{\tilde{W}}^{\lambda }}{{\langle {Q}_{2}\rangle }^{\lambda }}=\displaystyle \frac{{\tilde{W}}^{0}+\tfrac{{\hslash }}{2}[(1-{Q}_{{AB}}^{* }){\omega }_{2}+(1-{Q}_{{CD}}^{* }){\omega }_{1}]{C}_{N,\lambda }}{{\langle {Q}_{2}\rangle }^{0}+\tfrac{{\hslash }{\omega }_{2}}{2}(1-{Q}_{{AB}}^{* }){C}_{N,{\rm{\Lambda }}}}.\end{eqnarray*}$$

The first derivative of the above quantity with respect to λ is given by

$$\begin{eqnarray*}{\partial }_{\lambda }{\eta }^{\lambda } & = & \displaystyle \frac{\tfrac{{\hslash }}{2}[(1-{Q}_{{AB}}^{* }){\omega }_{2}+(1-{Q}_{{CD}}^{* }){\omega }_{1}]N(N+1)}{{\langle {Q}_{2}\rangle }^{\lambda }}-\displaystyle \frac{{\hslash }{\omega }_{2}}{2}(1-{Q}_{{AB}}^{* })N(N+1)\displaystyle \frac{{\tilde{W}}^{\lambda }}{{({\langle {Q}_{2}\rangle }^{\lambda })}^{2}},\\ & = & \displaystyle \frac{{\hslash }{\omega }_{1}(1-{Q}_{{CD}}^{* })N(N+1)}{2{\langle {Q}_{2}\rangle }^{\lambda }}+\displaystyle \frac{{\hslash }{\omega }_{2}(1-{Q}_{{CD}}^{* })N(N+1)}{2{\langle {Q}_{2}\rangle }^{\lambda }}(1-{\eta }^{\lambda })\leqslant 0.\end{eqnarray*}$$

It follows that for the adiabatic case (${Q}_{{AB}}^{* }={Q}_{{CD}}^{* }=1$) the efficiency does not depends on λ, i.e., ${\eta }^{\lambda }={\eta }^{0}$.

Next we show that the optimal work is a monotonically increasing function of λ. For the sake of simplicity, we introduce the parameter $x={\omega }_{1}/{\omega }_{2}$. We have ${\tilde{W}}^{\lambda }(x)={\tilde{W}}^{0}(x)+\tfrac{{\hslash }{\omega }_{1}}{2}\left[(1-{Q}_{{AB}}^{* })\tfrac{1}{x}+(1-{Q}_{{CD}}^{* })\right]{C}_{N,\lambda }$. We proved that $\forall x\in [0,1]$, ${\partial }_{\lambda }{\tilde{W}}^{\lambda }\leqslant 0$ which means that $\forall x\in [0,1]$ and $\forall \lambda ,\lambda ^{\prime}$ such that $\lambda \leqslant \lambda ^{\prime}$, ${\tilde{W}}^{\lambda }(x)\geqslant {\tilde{W}}^{\lambda ^{\prime} }(x)$, and so $\forall x\in [0,1],{\max }_{y\in [\mathrm{0,1}]}{\tilde{W}}^{\lambda }(y)\geqslant {\tilde{W}}^{\lambda }(x)$, thus if there exists ${x}_{{\rm{opt}}}^{\lambda }$ and ${x}_{{\rm{opt}}}^{\lambda ^{\prime} }$ such that ${\tilde{W}}^{\lambda }({x}_{{\rm{opt}}}^{\lambda })={\max }_{y\in [\mathrm{0,1}]}{\tilde{W}}^{\lambda }(y)\lt \infty$ and ${\tilde{W}}^{\lambda ^{\prime} }({x}_{{\rm{opt}}}^{\lambda ^{\prime} })={\max }_{x\in [\mathrm{0,1}]}{\tilde{W}}^{\lambda ^{\prime} }(x)\lt \infty$, then ${\tilde{W}}_{{\rm{opt}}}^{\lambda }\geqslant {\tilde{W}}_{{\rm{opt}}}^{\lambda ^{\prime} }$.

We introduce two parameters

$$\begin{eqnarray}&&x=\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}},a=\displaystyle \frac{{\beta }_{{\rm{h}}}}{{\beta }_{{\rm{c}}}},\end{eqnarray} \tag{ B14 }$$

with $a\lt 1,x\lt 1$.

To quantify the boost in the performance of the QHE resulting from many-particle effects, we introduce the following quantities

$$\begin{eqnarray}&&{r}_{N,\lambda }(x,a)\equiv \displaystyle \frac{{W}_{N,\lambda }(x,a)}{N\times {W}_{1}(x,a)},\end{eqnarray} \tag{ B15a }$$

$$\begin{eqnarray}&&{\rho }_{N,\lambda }(x,a)\equiv \displaystyle \frac{{\eta }_{N,\lambda }(x,a)}{{\eta }_{1}(x,a)}.\end{eqnarray} \tag{ B15b }$$

The first one is the the ratio between the total output work ${W}_{N,\lambda }(x,a)$ of a many-particle QHE and the corresponding value of $N$ independent single-particle QHEs given by $N\times {W}_{1}(x,a)$ under the same resources, that is, for the same values of $N,\lambda ,{\beta }_{{\rm{c,h}}},{\omega }_{\mathrm{1,2}}$. The ratio ${\rho }_{N,\lambda }(x,a)$ between the corresponding efficiencies is defined analogously. We note that the ratio between output power is equal to ${\rho }_{N,\lambda }(x,a)$ as the cycle time is the same for both systems.

Assuming that the temperature of the cold reservoir is high ${\hslash }{\omega }_{1}{\beta }_{{\rm{c}}}\ll 1$ the total output work as well as the efficiency of a single-particle QHE are dominated by the classical values, see equation (B5) (taking $N=1$) and (9), as well as the appendix D for a sudden quench protocol. This is no longer the case for the many-particle QHE as the Planck constant is scaled up by the number of particles and quantum effects become comparable with the classical ones, see equation (B6). In this large-$N$ regime, we observe numerically that both ratios in equation (B15) are greater than one. In other words, the performance of the many-particle QHE can surpass that of an ensemble of single-particle QHEs with the same resources. We call this phenomenon quantum supremacy. In the main body of this article we look at an alternative definition of quantum supremacy considering the many-particle QHE and the series of single-particle QHEs with the same number of particles $N$, temperatures ${\beta }_{{\rm{c,h}}}$, trap frequency ${\omega }_{1}$, and interaction strength λ but with a frequency ratio $x={\omega }_{1}/{\omega }_{2}$ that corresponds to the value optimizing the output power of the respective engines.

Here we show that for the same resources (i.e. when the value of ${\omega }_{2}$ is also fixed) the many-particle efficiency is always greater than the classical single-particle efficiency for a given range of the interaction strength λ. While no such general statement holds for the total output work, regimes exhibiting quantum supremacy can be found, showing the robustness of this phenomenon. We consider ${\hslash }{\omega }_{1}{\beta }_{{\rm{c}}}\ll 1$ and ${\sigma }_{{\rm{c}}},{\sigma }_{{\rm{h}}}\gt 1$, meaning that $N$ is large compared to 1. Then, from equations (B6), (B9)–(B11) we find

$$\begin{eqnarray}&&{W}_{N,\lambda }=\displaystyle \frac{N}{{\beta }_{{\rm{c}}}}\left[\left(1-\displaystyle \frac{{Q}_{{AB}}^{* }}{x}\right){\mu }_{\lambda }({\sigma }_{{\rm{c}}})+(1-{{xQ}}_{{CD}}^{* })\displaystyle \frac{1}{a}{\mu }_{\lambda }({\sigma }_{{\rm{h}}})\right],\end{eqnarray} \tag{ B16a }$$

$$\begin{eqnarray}&&N\times {W}_{1}=\displaystyle \frac{N}{{\beta }_{{\rm{c}}}}\left[\left(1-\displaystyle \frac{{Q}_{{AB}}^{* }}{x}\right)+(1-{{xQ}}_{{CD}}^{* })\displaystyle \frac{1}{a}\right],\end{eqnarray} \tag{ B16b }$$

and

$$\begin{eqnarray}&&{\eta }_{N,\lambda }=\displaystyle \frac{\left(1-\tfrac{{Q}_{{AB}}^{* }}{x}\right){\mu }_{\lambda }({\sigma }_{{\rm{c}}})+(1-{{xQ}}_{{CD}}^{* })\tfrac{1}{a}{\mu }_{\lambda }({\sigma }_{{\rm{h}}})}{\tfrac{1}{a}{\mu }_{\lambda }({\sigma }_{{\rm{h}}})-\tfrac{{Q}_{{AB}}^{* }}{x}{\mu }_{\lambda }({\sigma }_{{\rm{c}}})},\end{eqnarray} \tag{ B17a }$$

$$\begin{eqnarray}&&{\eta }_{1}=\displaystyle \frac{\left(1-\tfrac{{Q}_{{AB}}^{* }}{x}\right)+(1-{{xQ}}_{{CD}}^{* })\tfrac{1}{a}}{\tfrac{1}{a}-\tfrac{{Q}_{{AB}}^{* }}{x}}.\end{eqnarray} \tag{ B17b }$$

As a result

$$\begin{eqnarray}&&{r}_{N,\lambda }(x,a)={\mu }_{\lambda }({\sigma }_{{\rm{c}}})\displaystyle \frac{{W}_{1}(x,{{af}}_{\lambda })}{{W}_{1}(x,a)},\end{eqnarray} \tag{ B18a }$$

$$\begin{eqnarray}&&{\rho }_{N,\lambda }(x,a)=\displaystyle \frac{{\eta }_{1}(x,{{af}}_{\lambda })}{{\eta }_{1}(x,a)},\end{eqnarray} \tag{ B18b }$$

where

$$\begin{eqnarray}&&{f}_{\lambda }=\displaystyle \frac{{\mu }_{\lambda }({\sigma }_{{\rm{c}}})}{{\mu }_{\lambda }({\sigma }_{{\rm{h}}})},\end{eqnarray} \tag{ B19 }$$

with ${\sigma }_{{\rm{h}}}=\tfrac{a}{x}{\sigma }_{{\rm{c}}}\lt {\sigma }_{{\rm{c}}}$. In equations (B18a) and (B18b) we see that the total work and the efficiency of a many-particle QHE can be written by rescaling the ratio of temperatures $a={\beta }_{{\rm{h}}}/{\beta }_{{\rm{c}}}$ in the single-particle quantities W₁ and ${\eta }_{1}$. One finds that

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial a}{W}_{1}(x,a)=-\displaystyle \frac{(1-{{xQ}}_{{CD}}^{* })}{{a}^{2}}\lt 0,\mathrm{for}\ \mathrm{all}\ a,x,\end{eqnarray} \tag{ B20a }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial a}{\eta }_{1}(x,a)=-\displaystyle \frac{({Q}_{{AB}}^{* }{Q}_{{CD}}^{* }-1){x}^{2}}{{(x-{{aQ}}_{{AB}}^{* })}^{2}}\lt 0,\mathrm{for}\ \mathrm{all}\ a,x.\end{eqnarray} \tag{ B20b }$$

From equations (B20b) and (B19) we find that the efficiency ratio (B18b) exceeds unity if the function ${\mu }_{\lambda }(\sigma )$ defined by (B7) is decreasing. We know that it is monotonically decreasing for $\lambda =0$ and monotonically increasing for $\lambda \geqslant 1/2$. This means that for the ideal Bose gas the many-particle efficiency is always greater than the corresponding value for an ensemble of single-particle QHEs in the classical regime. For $0\lt \lambda \lt 1/2$, the function ${\mu }_{\lambda }(\sigma )$ decreases if σ is less than a critical value that increases with λ, see figure B1. Consequently for $\sigma \leqslant {\sigma }_{{\rm{c}}}$, ${\mu }_{\lambda }(\sigma )$ is a decreasing function of σ if the interaction strength is below a critical value ${\lambda }_{{\rm{c}}}$ which is typically less than 0.2 (we obtained this estimate after numerical computations). Thus the ratio ${\rho }_{N,\lambda }(x,a)$ is greater than 1 for $\lambda \leqslant {\lambda }_{{\rm{c}}}$ and less than 1 otherwise. In other words, many-particle quantum effects enhance the efficiency of a QHE for weak interactions.

Given our focus on quantum supremacy, we would also like to show the condition to improve the output power of the engine, i.e., for a ratio (B18a) greater than unity. From equation (B20a) and a similar argument to that used above, we know that the ratio between ${W}_{1}(x,{{af}}_{\lambda })$ and ${W}_{1}(x,a)$ is greater than 1 for $\lambda =0$ and for $0\lt \lambda \leqslant {\lambda }_{{\rm{c}}}$, where ${\lambda }_{{\rm{c}}}$ is the same critical value discussed in the previous paragraph. However, in equation (B18a) we find that this ratio is multiplied by a factor ${\sigma }_{\lambda }({\sigma }_{{\rm{c}}})$ which is less than 1 for $\lambda \leqslant {\lambda }_{{\rm{c}}}$, making the analysis more complex. Yet, we find that for ${\sigma }_{{\rm{c}}}$ not too large (typically ${\sigma }_{{\rm{c}}}\leqslant 5-10$, depending on the values of a and x) the output power of a many-particle QHE surpasses the corresponding value of an ensemble of single-particle QHEs, exhibiting quantum supremacy.

For the sudden-quench driving between ${\omega }_{1}$ and ${\omega }_{2}$ the nonadiabatic factor ${Q}^{* }$ takes the same value in both the compression and expansion strokes and is given by

$$\begin{eqnarray}&&{Q}^{* }=\displaystyle \frac{{\omega }_{1}^{2}+{\omega }_{2}^{2}}{2{\omega }_{1}{\omega }_{2}}.\end{eqnarray} \tag{ D1 }$$

We recall the condition over $a={\beta }_{{\rm{h}}}/{\beta }_{{\rm{c}}}$ and $x={\omega }_{1}/{\omega }_{2}$ (C2) which can be read as $a\leqslant x$. One can derive another condition which is more useful for small x:

$$\begin{eqnarray}&&\displaystyle \frac{a}{x}\leqslant 2\,x,\mathrm{for}\ x\ll 1.\end{eqnarray} \tag{ D2 }$$

To derive this condition we use (B9b) to show that $\langle {Q}_{2}\rangle =\tfrac{{\omega }_{2}}{2}N(1+\lambda (N-1))(1-{Q}^{* })+{\sum }_{k=1}^{N}k{\hslash }{\omega }_{2}\left(\tfrac{1}{{{\rm{e}}}^{k{\hslash }{\beta }_{{\rm{h}}}{\omega }_{2}}-1}-\tfrac{{Q}^{* }}{{{\rm{e}}}^{k{\hslash }{\beta }_{{\rm{c}}}{\omega }_{1}}-1}\right)$. We must have $\langle {Q}_{2}\rangle \geqslant 0$, and the sum in the previous equation must be positive since the first term is negative in the adiabatic case (${Q}^{* }\geqslant 1$). Using $x\leqslant {{\rm{e}}}^{x}-1\leqslant 2x$ for $x\leqslant 1$, one can easily show that an upper bound to the sum is positive if $\tfrac{{\beta }_{{\rm{c}}}{\omega }_{1}}{{\beta }_{{\rm{h}}}{\omega }_{2}}\geqslant 2{Q}^{* }$ which gives (D2) for small x. This condition is weak in the sense that it does not ensure that $\langle {Q}_{2}\rangle$ is positive. However, we know that if this condition is not satisfied then $\langle {Q}_{2}\rangle$ is negative.

Using (B10) and (D1), we find

$$\begin{eqnarray}&&{\tilde{W}}_{\mathrm{sq}}=\left(\displaystyle \frac{{x}^{2}-1}{2{x}^{2}}\right){\langle H\rangle }_{A}+\left(\displaystyle \frac{1-{x}^{2}}{2}\right){\langle H\rangle }_{C},\end{eqnarray} \tag{ D3 }$$

where ${\langle H\rangle }_{A}={{ \mathcal E }}_{N}({\omega }_{1},{\beta }_{{\rm{c}}})$ and ${\langle H\rangle }_{C}={{ \mathcal E }}_{N}({\omega }_{1}/x,a{\beta }_{{\rm{c}}})$.

We want to find a condition on x so that the work is maximal. To this end, we solve this equation

$$\begin{eqnarray}&&{\partial }_{x}{\tilde{W}}_{\mathrm{sq}}=0,\end{eqnarray} \tag{ D4 }$$

which gives the maximal work since the total work is a concave function of x, see (D3).

By (B5), in the high-temperature limit of the hot reservoir ${\sigma }_{{\rm{h}}}=N{\hslash }{\beta }_{{\rm{h}}}{\omega }_{2}\ll 1$, the total work is given by

$$\begin{eqnarray}&&{\tilde{W}}_{\mathrm{sq}}=\left(\displaystyle \frac{{x}^{2}-1}{2{x}^{2}}\right){\langle H\rangle }_{A}+\left(\displaystyle \frac{1-{x}^{2}}{2}\right)\displaystyle \frac{N}{a{\beta }_{{\rm{c}}}}\left(1+\displaystyle \frac{a}{2x}{g}_{N}\left(\lambda -\displaystyle \frac{1}{2}\right){\sigma }_{{\rm{c}}}+O({\sigma }_{{\rm{c}}}^{2})\right),\end{eqnarray} \tag{ D5 }$$

where we have kept only the first order correction term.

If we assume that ${\sigma }_{{\rm{h}}}\ll 1$ and $a={\beta }_{{\rm{h}}}/{\beta }_{{\rm{c}}}$ small, the optimal value for the frequency ratio is given by

$$\begin{eqnarray}&&{\bar{x}}_{\mathrm{sq}}={\left(\displaystyle \frac{a{\beta }_{{\rm{c}}}{\langle H\rangle }_{A}}{N}\right)}^{1/4}={\left(\displaystyle \frac{{\beta }_{{\rm{h}}}{\langle H\rangle }_{A}}{N}\right)}^{1/4}.\end{eqnarray} \tag{ D6 }$$

Thus, one obtains the following expression for the optimal output work and efficiency

$$\begin{eqnarray}&&{\tilde{W}}_{\mathrm{sq}}^{(N,\lambda )}=\displaystyle \frac{N}{2a{\beta }_{{\rm{c}}}}\left[{(1-\sqrt{{\alpha }_{{\rm{sq}}}})}^{2}+\left(\displaystyle \frac{1-\sqrt{{\alpha }_{{\rm{sq}}}}}{2{\alpha }_{{\rm{sq}}}^{1/4}}\right){g}_{N}\left(\lambda -\displaystyle \frac{1}{2}\right)a{\sigma }_{{\rm{c}}}+O({\sigma }_{{\rm{c}}}^{2})\right],\end{eqnarray} \tag{ D7a }$$

$$\begin{eqnarray}&&{\eta }_{\mathrm{sq}}^{(N,\lambda )}=\displaystyle \frac{1-\sqrt{{\alpha }_{\mathrm{sq}}}+\tfrac{a{\sigma }_{{\rm{c}}}}{2{\alpha }_{\mathrm{sq}}^{1/4}}{g}_{N}\left(\lambda -\tfrac{1}{2}\right)}{2+\sqrt{{\alpha }_{\mathrm{sq}}}+\tfrac{a{\sigma }_{{\rm{c}}}}{{\alpha }_{\mathrm{sq}}^{1/4}(1-\sqrt{{\alpha }_{\mathrm{sq}}})}{g}_{N}\left(\lambda -\tfrac{1}{2}\right)},\end{eqnarray} \tag{ D7b }$$

where we have introduced ${\alpha }_{\mathrm{sq}}=\tfrac{a{\beta }_{{\rm{c}}}{\langle H\rangle }_{A}}{N}$ and performed a very similar computation for the heat $\langle {Q}_{2}\rangle$.

Explicit expressions for the optimal work at different regimes of inverse temperature ${\beta }_{{\rm{c}}}$ can be found:

• For ${\sigma }_{{\rm{c}}}\gg N$, one has ${\bar{x}}_{{sq}}\approx {\kappa }_{N,\lambda }\sqrt{{\hslash }{\beta }_{{\rm{h}}}{\omega }_{1}}$, where ${\kappa }_{N,\lambda }=\sqrt{(1+(N-1)\lambda )}$ for $N\geqslant 1$. Then
  $$\begin{eqnarray}&&{\tilde{W}}_{\mathrm{sq}}^{(N,\lambda )}=\displaystyle \frac{N}{2a{\beta }_{{\rm{c}}}}{(1-{\kappa }_{N,\lambda }\sqrt{{\hslash }{\beta }_{{\rm{h}}}{\omega }_{1}/2})}^{2},\end{eqnarray} \tag{ D8a }$$
  $$\begin{eqnarray}&&{\eta }_{\mathrm{sq}}^{(N,\lambda )}=\displaystyle \frac{1-{\kappa }_{N,\lambda }\sqrt{{\hslash }{\beta }_{{\rm{h}}}{\omega }_{1}/2}}{2+{\kappa }_{N,\lambda }\sqrt{{\hslash }{\beta }_{{\rm{h}}}{\omega }_{1}/2}}.\end{eqnarray} \tag{ D8b }$$
• For large temperature ${\sigma }_{{\rm{c}}}\ll 1$,
  • for $N=1$, we have ${\bar{x}}_{\mathrm{sq}}\approx {a}^{1/4}\left(1+\tfrac{1}{48}{\sigma }_{{\rm{c}}}^{2}\right)$, and
    $$\begin{eqnarray}&&{\tilde{W}}_{\mathrm{sq}}^{(1,\lambda )}=\displaystyle \frac{1}{{\beta }_{{\rm{c}}}}\left[\displaystyle \frac{{(1-\sqrt{a})}^{2}}{2a}+\displaystyle \frac{(-1+\sqrt{a}+a-{a}^{3/2})}{24\sqrt{a}}{\sigma }_{{\rm{c}}}^{2}\right],\end{eqnarray} \tag{ D9a }$$
    $$\begin{eqnarray}&&{\eta }_{\mathrm{sq}}^{(1,\lambda )}=\displaystyle \frac{1-\sqrt{a}}{2+\sqrt{a}}-\displaystyle \frac{(3-2{a}^{3/2})\sqrt{a}}{24{(2+\sqrt{a})}^{2}}{\sigma }_{{\rm{c}}}^{2}.\end{eqnarray} \tag{ D9b }$$
  • for $N\gt 1$, ${\bar{x}}_{\mathrm{sq}}\approx {a}^{1/4}\left(1+\tfrac{1}{8}{\sigma }_{{\rm{c}}}{g}_{N}\left(\lambda -\tfrac{1}{2}\right)\right)+O({\sigma }_{{\rm{c}}}^{2})$ and
    $$\begin{eqnarray}&&{\tilde{W}}_{\mathrm{sq}}^{(N,\lambda )}=\displaystyle \frac{N}{{\beta }_{{\rm{c}}}}\left[\displaystyle \frac{{(1-\sqrt{a})}^{2}}{2a}+\left(\displaystyle \frac{1-\sqrt{a}}{4\sqrt{a}}\right)(1-{a}^{1/4})\left(\displaystyle \frac{1}{2}-\lambda \right){g}_{N}{\sigma }_{{\rm{c}}}+O({\sigma }_{{\rm{c}}}^{2})\right],\end{eqnarray} \tag{ D10a }$$
    $$\begin{eqnarray}&&{\eta }_{\mathrm{sq}}^{(N,\lambda )}=\displaystyle \frac{1-\sqrt{a}}{2+\sqrt{a}}+\displaystyle \frac{(3-2{a}^{3/4})\sqrt{a}}{4{(2+\sqrt{a})}^{2}}\left(\displaystyle \frac{1}{2}-\lambda \right){g}_{N}{\sigma }_{{\rm{c}}}+O({\sigma }_{{\rm{c}}}^{2}).\end{eqnarray} \tag{ D10b }$$
• For $N\gt 1$, we consider an intermediate regime corresponding to a large temperature ${\hslash }{\beta }_{{\rm{h}}}{\omega }_{1}\ll 1$ but a relatively small temperature per particle ${\sigma }_{{\rm{c}}}=N{\hslash }{\beta }_{{\rm{c}}}{\omega }_{1}\geqslant 1$, which means that the number of particles is large ${\sigma }_{{\rm{h}}}\ll 1$ (i.e., ${\beta }_{2}\ll {\beta }_{1}$). Using (B6), we find $\tfrac{a{\beta }_{{\rm{c}}}{\langle H\rangle }_{A}}{N}\approx a{\mu }_{\lambda }({\sigma }_{{\rm{c}}})$ where ${\mu }_{\lambda }({\sigma }_{{\rm{c}}})$ is given by (B7). Notice that ${\mu }_{0}({\sigma }_{{\rm{c}}})$ is a decreasing function of ${\sigma }_{{\rm{c}}}$ with the following bounds ${\pi }^{2}/(6N{\hslash }{\omega }_{1}{\beta }_{{\rm{c}}})\leqslant \mu ({\sigma }_{{\rm{c}}})\leqslant 1$. If one considers ${\sigma }_{{\rm{c}}}$ not too large (so that ${\sigma }_{{\rm{h}}}=N{\hslash }{\beta }_{{\rm{h}}}{\omega }_{2}\ll 1$ still holds), by (D6) and (D7b) one finds
  $$\begin{eqnarray}&&{\tilde{W}}_{\mathrm{sq}}^{(N,\lambda )}=\displaystyle \frac{N}{2a{\beta }_{{\rm{c}}}}{(1-\sqrt{{\mu }_{\lambda }({\sigma }_{{\rm{c}}})a})}^{2},\end{eqnarray} \tag{ D11a }$$
  $$\begin{eqnarray}&&{\eta }_{\mathrm{sq}}^{(N,\lambda )}=\displaystyle \frac{1-\sqrt{{\mu }_{\lambda }({\sigma }_{{\rm{c}}})a}}{2+\sqrt{{\mu }_{\lambda }({\sigma }_{{\rm{c}}})a}},\end{eqnarray} \tag{ D11b }$$
  where ${\mu }_{\lambda }({\sigma }_{{\rm{c}}})$ is defined by (B7). Thus, the optimal output power as well as the efficiency at optimal output power can be substantially improved for $\lambda \lt 1/2$ (and for ${\sigma }_{{\rm{c}}}$ small enough) as ${\mu }_{\lambda }\leqslant 1$ and decreases otherwise (i.e., for $\lambda \gt 1/2$), see figures B1–E1. Notice that (D11a) and (D11b) are valid only for $a\ll 1$ as ${\sigma }_{{\rm{h}}}\ll 1$. The case ${\sigma }_{{\rm{h}}}\geqslant 1$ will be studied in the next section.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Numerically, we observe that for an ideal Bose gas (the optimal working medium) the efficiency grows up when $N$ becomes large but saturates when $N{\hslash }{\beta }_{{\rm{h}}}{\omega }_{2}\gt 10$. To understand physically this curious phenomenon, it suffices to write the expression of the total work for large particle number $N\gg {({\hslash }{\beta }_{{\rm{c(h)}}}{\omega }_{1(2)})}^{-1}$ considering ${\sigma }_{{\rm{c(h)}}}\ll N$ (i.e., ${\hslash }{\beta }_{{\rm{c}}}{\omega }_{1}\ll 1$) as we know that in the very-low temperature regime of the cold reservoir (when ${\hslash }{\beta }_{{\rm{c}}}{\omega }_{1}\gg 1$) the engine does not operate efficiently. In this novel regime, we find that the total work does not depend on $N$,

$$\begin{eqnarray}&&{\tilde{W}}_{\mathrm{sq}}=\displaystyle \frac{{\pi }^{2}}{12{\hslash }{\beta }_{{\rm{c}}}^{2}{\omega }_{1}}\left[\displaystyle \frac{{x}^{2}-1}{{x}^{2}}+(1-{x}^{2})\displaystyle \frac{x}{{a}^{2}}\right].\end{eqnarray} \tag{ E1 }$$

To derive this equation we use (B6) and the fact that ${\mu }_{0}({\sigma }_{{\rm{c}}})\approx \tfrac{{\pi }^{2}}{6{\sigma }_{{\rm{c}}}}$ which leads to ${\langle H\rangle }_{A}\approx \tfrac{N{\pi }^{2}}{6{\beta }_{{\rm{c}}}{\sigma }_{{\rm{c}}}}=\tfrac{{\pi }^{2}}{6{\hslash }{\beta }_{{\rm{c}}}^{2}{\omega }_{1}}$ and ${\langle H\rangle }_{C}\approx \tfrac{N{\pi }^{2}}{6N{\hslash }{\beta }_{{\rm{h}}}^{2}{\omega }_{2}}=\tfrac{{\pi }^{2}x}{6a{\hslash }{\beta }_{{\rm{c}}}^{2}{\omega }_{1}}$ as ${\beta }_{{\rm{h}}}=a{\beta }_{{\rm{c}}}$ and ${\omega }_{2}={\omega }_{1}/x$. We have neglected the ground energy contribution assuming that $N{{\hslash }}^{2}{\beta }_{{\rm{c(h)}}}^{2}{\omega }_{1(2)}^{2}\ll 1$. In the case where the ground state energy dominates, it is clear that the medium is not efficient as it is effectively 'frozen' (as previously discussed). As a result, we consider a large particle number $N$ compared to ${({\hslash }{\beta }_{{\rm{c}}}{\omega }_{1})}^{-1}$ but very small compared to ${({\hslash }{\beta }_{{\rm{c}}}{\omega }_{1})}^{-2}$. In figure E2 we take ${\hslash }{\beta }_{{\rm{c}}}{\omega }_{1}=0.01$ and $N=1000$ that corresponds to this regime as in this case $100\ll N\ll 10000$. The work done ${\tilde{W}}_{\mathrm{sq}}$ is of order

$$\begin{eqnarray*}&&\displaystyle \frac{{\pi }^{2}\langle n\rangle }{12{\beta }_{{\rm{h}}}},\end{eqnarray*}$$

where $\langle n\rangle ={({\hslash }{\beta }_{{\rm{h}}}{\omega }_{2})}^{-1}$ is the average number of excited states per particle for the equilibrium state C. We point out that the number of excited states per particle remains small compared to the number of particle in this regime $\langle n\rangle \ll N$ so that the total work is also small compared to the classical one as ${\tilde{W}}_{\mathrm{sq}}\ll \tfrac{N}{2{\beta }_{{\rm{h}}}}$.

In this regime, one can find the frequency ratio for which the work is optimal by solving this polynomial equation

$$\begin{eqnarray*}&&{p}_{a}(x)\equiv 3{x}^{5}-{x}^{3}-2{a}^{2},0\lt a\lt 1,\end{eqnarray*}$$

where the solution only depends on the temperature ratio a. Notice that we expect a to be greater than 0.1 to keep $N\gg {({\hslash }{\beta }_{{\rm{h}}}{\omega }_{2})}^{-1}$. The solution of this equation can not be expressed as a standard solution of a polynomial equation since the degree of the equation is greater than 4. However, numerical analysis shows that an estimate of the real solution takes the form

$$\begin{eqnarray*}&&{\bar{x}}_{\mathrm{sq}}\approx {a}^{s},1/4\leqslant s\leqslant 1/3,\end{eqnarray*}$$

where ${a}^{1/4}$ (respectively ${a}^{1/3}$) gives the upper bound (respectively lower bound) of the root for $a\gt 0.1$.

Hence, after computing the heat $\langle {Q}_{2}\rangle$ we obtain this expression for the efficiency, see figure E2,

$$\begin{eqnarray}&&{\eta }_{\mathrm{sq}}^{(N,0)}\approx (1-{a}^{2s})\displaystyle \frac{1-{a}^{2-3s}}{2-{a}^{2-3s}-{a}^{2-s}},1/4\leqslant s\leqslant 1/3,\end{eqnarray} \tag{ E2 }$$

where the limit case $s=1/4$ (respectively $s=1/3$) gives a lower bound (respectively upper bound) to the numerical exact computations for $N\geqslant 1000$. Numerically (see figure 2) we show that the efficiency can be enhanced up to 50%, see figure C1, and figure E2. We further notice that from (E2), and (D9b), the efficiency in this regime is upper-bounded by 150% of the value of the single-particle efficiency.

While for the ideal Bose gas the λ-term vanishes, for the general case the ground state contribution of λ can reduce the enhancement of the efficiency. This observation follows from the strong dependence of the total work on λ,

$$\begin{eqnarray}&&{\tilde{W}}_{\mathrm{sq}}\approx \displaystyle \frac{{x}^{2}-1}{2{x}^{2}}\displaystyle \frac{N}{{\beta }_{{\rm{c}}}}{\mu }_{\lambda }({\sigma }_{{\rm{c}}})+\displaystyle \frac{1-{x}^{2}}{2}\displaystyle \frac{N}{{\beta }_{{\rm{h}}}}{\mu }_{\lambda }({\sigma }_{{\rm{h}}}),\mathrm{for}\ 1\leqslant {\sigma }_{{\rm{c,h}}}\ll N\end{eqnarray} \tag{ E3a }$$

$$\begin{eqnarray}&&\approx \,\displaystyle \frac{{\pi }^{2}\langle n\rangle }{12{\beta }_{{\rm{h}}}}\left[{a}^{2}\left(\displaystyle \frac{{x}^{2}-1}{{x}^{2}}\right)\left(1+\displaystyle \frac{3\lambda {\sigma }_{{\rm{c}}}^{2}}{{\pi }^{2}}\right)+x(1-{x}^{2})\left(1+\displaystyle \frac{3{a}^{2}\lambda {\sigma }_{{\rm{c}}}^{2}}{{\pi }^{2}{x}^{2}}\right)\right],\mathrm{for}\ 1\ll {\sigma }_{{\rm{c,h}}}\ll N,\end{eqnarray} \tag{ E3b }$$

where $\langle n\rangle ={({\hslash }{\beta }_{{\rm{h}}}{\omega }_{2})}^{-1}$. To derive equation (E3b) we use (B6). In the regime of large number of particles (${\sigma }_{{\rm{c(h)}}}\gg 1$), we estimate ${\mu }_{\lambda }$ as

$$\begin{eqnarray*}&&{\mu }_{\lambda }({\sigma }_{{\rm{c}}})\approx \displaystyle \frac{{\pi }^{2}}{6{\sigma }_{{\rm{c}}}}+\displaystyle \frac{\lambda }{2}{\sigma }_{{\rm{c}}}=\displaystyle \frac{{\pi }^{2}}{6{\sigma }_{{\rm{c}}}}\left(1+\displaystyle \frac{3\lambda {\sigma }_{{\rm{c}}}^{2}}{{\pi }^{2}}\right).\end{eqnarray*}$$

Hence, one finds different regimes depending on the value of λ. For $\lambda \ll {\lambda }_{{\rm{c}}}={\pi }^{2}/(3{\sigma }_{{\rm{c}}}^{2})$ the thermal energy dominates and then the problem can be treated as the case with $\lambda =0$. On the contrary, for $\lambda \gg {\lambda }_{{\rm{c}}}$ the thermal energy becomes negligible, thus the quantum fluid works less efficiently and becomes 'frozen' for $\lambda =1/2$ and 1, see figures C1–E1 .

Zoom In Zoom Out Reset image size

To conclude, for a many-particle QHE in the low temperature regime nonadiabatic effects can enhance the efficiency over the single-particle value for $\lambda =0$ and for small λ (typically $\lambda \ll {\lambda }_{{\rm{c}}}$). However, for $\lambda \geqslant 1/2$ the many-particle QHE is less efficient than the single-particle heat engine. In this case it is better to engineer a series of QHE with lower number of particles.

Under adiabatic driving, the total work is given by

$$\begin{eqnarray}&&{\tilde{W}}_{\mathrm{ad}}=\displaystyle \frac{{\pi }^{2}}{6{\hslash }{\beta }_{{\rm{c}}}^{2}{\omega }_{1}}\left[\displaystyle \frac{x-1}{x}+(1-x)\displaystyle \frac{x}{{a}^{2}}\right].\end{eqnarray} \tag{ E4 }$$

The derivation of equation (E4) is very similar to (E1).

The frequency ratio for which the work is optimal is determined by solving this polynomial equation

$$\begin{eqnarray*}&&{p}_{a}(x)\equiv 2{x}^{3}-{x}^{2}-{a}^{2},0\lt a\lt 1,\end{eqnarray*}$$

whose solution depends only on the temperature ratio a. Here we take a to be greater than 0.1 so that $N\gg {({\hslash }{\beta }_{{\rm{h}}}{\omega }_{2})}^{-1}$. One can find an analytic expression of the real solution of this polynomial equation and provide an estimate for $a\geqslant 0.1$,

$$\begin{eqnarray*}&&{\bar{x}}_{\mathrm{ad}}\approx 1+\displaystyle \frac{a-1}{2}-\displaystyle \frac{{(a-1)}^{2}}{16}.\end{eqnarray*}$$

Inserting this value in the Otto efficiency we obtain

$$\begin{eqnarray}&&{\eta }_{\mathrm{ad}}^{(N)}\approx \displaystyle \frac{9-10a+{a}^{2}}{16}.\end{eqnarray} \tag{ E5 }$$

Combining (C8b) valid for small a and (E5) we find that the efficiency is minimum for $a\sim 0.1$ and is equal to 70% of the single-particle efficiency which is consistent with the numerical exact calculation shown in figure C1.