Quantum Supremacy of Many-Particle Thermal Machines

While the emergent field of quantum thermodynamics has the potential to impact energy science, the performance of thermal machines is often classical. We ask whether quantum effects can boost the performance of a thermal machine to reach quantum supremacy, i.e., surpassing both the efficiency and power achieved in classical thermodynamics. To this end, we introduce a nonadiabatic quantum heat engine operating an Otto cycle with a many-particle working medium, consisting of an interacting Bose gas confined in a time-dependent harmonic trap. It is shown that thanks to the interplay of nonadiabatic and many-particle quantum effects, this thermal machine can outperform an ensemble of single-particle heat engines with same resources, demonstrating quantum supremacy in many-particle thermal machines.


I. INTRODUCTION
The interplay of quantum technology and foundations of physics has turned quantum thermodynamics into a blooming field [1,2]. With the miniaturization to the nanoscale, a need has emerged to understand and control the dynamics of thermal machines operating in the presence of both thermal and quantum fluctuations. Quantum heat engines (QHEs) constitute a prominent example, targeting the efficient conversion of heat into mechanical work. The relevance of quantum thermal machines is further strengthened by their use to describe both natural and artificial light harvesting systems [3][4][5][6][7]. While quantum thermodynamics has the potential to revolutionize energy science, quantum thermal machines studied to date often exhibit a classical behavior.
Current efforts towards the realization of a tunable QHE in the laboratory use a single-particle working medium, e.g., a confined ion in a modified Paul trap [8][9][10]. The optimization of this type of single-particle QHE has received considerable attention [8,9,[11][12][13][14][15][16][17]. By contrast, the performance of a QHE with a many-particle working medium remains essentially unexplored [18][19][20][21]. This is however a timely question motivated by the prospects of scaling up QHE and related devices. In particular, an ion-trap realization of a QHE constitutes a natural testbed to explore many-particle effects with well-established quantum technology, e.g., using an ion chain as a working medium.
In this article, we pose the question as to whether there exist scenarios in which the performance of a thermal machine exhibits quantum supremacy, i.e., a superior performance to that achievable in classical thermodynamics. To address this question we introduce a model of a many-particle QHE and show that quantum supremacy can be achieved by exploiting the interplay of nonadiabatic dynamics and many-particle quantum effects. In particular, we identify the conditions for which the use of a manyparticle working medium leads to a simultaneous enhancement of the efficiency and output-power of the many-particle QHE, surpassing the values for an ensemble of single-particle QHEs with same resources.
The QHE we analyze operates an Otto cycle with a many-particle interacting quantum fluid in a time-dependent harmonic trap as a working medium that is alternately coupled to a hot and a cold reservoir, see Fig. 1. The Otto cycle is composed of four strokes shown in Fig. 1: 1) Adiabatic compression: Starting from a thermal state A with inverse temperature β c and decoupled from any thermal reservoir, the working medium is driven by increasing ω(t) from ω 1 to ω 2 until it reaches state B; 2) Hot isochore: Keeping the trap frequency constant, the working medium is coupled to a hot reservoir at inverse temperature β h , relaxing to a thermal state C; 3) Adiabatic expansion: The thermal state C is decoupled from the hot reservoir and driven as ω(t) decreases from ω 2 to ω 1 , reaching state D; 4) Cold isochore: Keeping the trap frequency ω 1 constant, the working medium is coupled to a cold reservoir at inverse temperature β c , relaxing to the initial thermal state A.
The efficiency of a heat engine is defined as the total output work per heat input η = − [ W 1 + W 3 ] / Q 2 . For the Otto cycle the input (output) work of the engine is given respectively by A quantum heat engine with an interacting Bose gas as a working medium. (a) Quantum Otto cycle for a working medium confined in a harmonic trap with frequency varying between ω 1 and ω 2 . States at A and C are thermal while those at B and D are generally nonequilibrium states. (b) Schematic representation of the working medium in a many-particle heat engine compared to an ensemble of independent single-particle QHEs illustrated in (c).

II. RESULTS
A. Nonadiabatic many-particle QHE We consider the working susbtance to be described by a quantum many-body Hamiltonian where N is the total number of particles and ω(t) is the trap frequency of the isotropic harmonic confinement. We assume V (r/b) = b 2 V (r) so that the unitary dynamics generated by (1) during the expansion and compression strokes exhibits scaleinvariance [22]. The time-evolution of an eigenstate Ψ at t = 0 is [23,24] Ψ(r 1 , . . . , r N ,t) = N e iφ Ψ r 1 b , . . . , where the time-dependent phase reads φ = ∑ i mḃr 2 i /(2hb) and we have neglected a global dynamical phase which plays no role in our analysis. The normalization constant is N = b −Nd/2 for a system in d spatial dimensions. The scaling factor b = b(t) > 0 is the solution of the Ermakov differential equationb + ω(t) 2 b = ω 2 0 /b 3 , where we choose the constant ω 0 = ω(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 is the nonadiabatic factor that accounts for the amount of energy excitations over the adiabatic dynamics, and b ad = [ω 0 /ω(t)] 1/2 is the scaling factor in the adiabatic limit, obtained from the Ermakov equation by settingb(t) ≈ 0. In this limit, Q * (t) reduces to unity and the mean energy is given by H(t) ad = H(0) /b 2 ad . From the stationarity of the initial state it follows that b(0) = 1 andḃ(0) = 0. Remarkably, Q * (t) coincides with the nonadiabatic factor introduced by Husimi for the single-particle timedependent harmonic oscillator [25], as we show in Appendix A 2. As a result, at the end of the strokes decoupled from the heat reservoirs, H B = Q * AB b −2 ad H A and H D = Q * CD b 2 ad H C , where b ad = (ω 1 /ω 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 W = W 1 + W 3 reads, From a practical point of view, the output power of the engine is more interesting as it accounts for the work per cycle P ≡ − W /(τ 1 + τ 2 + τ 3 + τ 4 ) , where τ i , i = 1, · · · , 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 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 τ = ∑ 4 i=1 τ i . In this limit (Q * AB(CD) → 1 or equivalently |ω/ω 2 | → 0, see [26]) equation (5) reduces to the Otto efficiency η O = 1 − ω 1 /ω 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 reduced as Q * AB(CD) ≥ 1, 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 η ≤ η 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 non-adiabatic 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.

B. 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], where λ ≥ 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 λ = 0 as well as hard-core bosons (Tonks-Girardeau gas) for λ = 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][33][34]. While the eigenstates are a paradigm of strong correlations with Bijl-Jastrow form, the spectrum is purely linear [28,35], where n j is the occupation number in the j-th mode of energy jhω, satisfying the normalization ∑ ∞ j=1 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 where the partition function for an ideal Bose gas Z (0) 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 H = −∂ β lnZ (λ ) 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) > 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 FIG. 2. Effect of the interactions on the performance of a nonadiabatic many-particle QHE compared to an ensemble of N singleparticle QHEs. (a) Efficiency ratio at optimal output power, see equation (8), as a function of β h /β c with β c = 0.01/(hω 1 ) and N = 200 (σ c = 2) for λ = 0, 0.2, 0.5, 1 from bottom to top. A sudden-quench protocol is considered. (b) Corresponding output power ratio under the same conditions. We observe that the ideal Bose gas λ = 0 as well as for the weakly interacting Bose gas λ = 0.2 one can boost the performance of a many-particle QHE for N large while it is better to engineer an ensemble of single-particle QHEs when the interaction strength λ ≥ 1/2. for nonadiabatic effects, we shall consider in the following a sudden-quench driving of the trap frequency between ω 1 and ω 2 . In this case, the nonadiabatic factor takes the form Q * AB(CD) = (ω 2 1 + ω 2 2 )/(2ω 1 ω 2 ) that is symmetric with respect to the exchange ω 1 ↔ ω 2 . As a result, the efficiency (5) of a realistic QHE with a short time per cycle has as fundamental upper limit (6), η sq ≤ η nad,O = (1 − (ω 1 /ω 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 (N,λ ) sq by optimizing the ratio ω 1 /ω 2 for fixed ω 1 , temperatures β c , β h , number of particles N, and interaction strength λ . We denote the efficiency at optimal output power by η (N,λ ) sq . To characterize many-particle effects, we introduce the following ratios 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 Fig. 2, where it is shown that they can be greater than 1 for λ < 1/2 and optimal for λ = 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 Fig. 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 Figs. 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. In the high-temperature limit of the hot reservoir σ h ≡ Nhβ h ω 2 1 and of the cold reservoir σ c ≡ Nhβ c ω 1 1, we find that the optimal frequency is given by where g N = (N − 1)/N and a ≡ β h /β c denotes the ratio of temperatures (see Appendix D for details), whence it follows that As a result, we recognize the leading term in η (N,λ ) sq , independent of σ c , as the Rezek-Kosloff efficiency shared by single-particle heat engines under sudden driving [11]. The leading quantum σ c -contribution to the efficiency at optimal work increases with the number of particles N for λ < 1/2. On the contrary, if λ > 1/2 quantum effects lower the efficiency and vanish for λ = 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) FIG. 3. Quantum supremacy of a nonadiabatic many-particle QHE compared to an ensemble of N single-particle QHEs. Efficiency ratio (left) and output power ratio (right) at optimal output power, see equation (8) suggests that for a QHE with a small particle number, quantum corrections can lead to quantum supremacy. Yet, r (N,λ ) sq and ρ (N,λ ) sq vary as 1 + O(σ 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 σ c N of the cold reservoir (keeping at high temperature the hot reservoir σ h N), the thermal energy of the equilibrium state A does not contribute to the total work. One has (ω 1 /ω 2 ) sq ≈ (κ 2 N,λh β h ω 1 ) 1/4 (see Appendix D for details), where κ N,λ = (1 + (N − 1)λ ). The efficiency at optimal output power reads 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), becausē hβ h ω 1 β h /β 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 σ c ∼ N and for σ c , σ h ≥ N. We have seen that estimate (9) is valid for small number of particles (e.g., N ∼ 1 − 30 for σ c /N = 0.01) and values of the ratio β h /β c consistent with the assumption σ h 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 ≤ σ c 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 (2) , and µ λ (σ c(h) ) represents the relative deviation to the classical value of the mean energy N/β 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 = − W /τ. 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 σ c ≥ 10, typically), see Appendix E for a derivation of this upper bound. In Fig. 3 we observe that for λ = 0 (free bosons) and for N ≤ 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 ≥ 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 β h /β c less than a critical value (∼ 0.2 for N = 500). For a weakly interacting Bose gas (typically λ ≤ 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 β h and β c .

III. DISCUSSION
In this article, we have shown that quantum thermal machines can exhibit quantum supremacy, i.e., a superior performance to that allowed in classical thermodynamics. To this end, we have introduced a many-particle quantum heat engine with an interacting Bose gas confined in a time-dependent harmonic trap as a working medium. Our analysis shows that quantum supremacy results from the interplay between many-particle quantum effects and nonadiabtic dynamics that boost the finite-time efficiency and output power of a many-particle QHE, surpassing that of an ensemble of single-particle heat engines, matched by the predictions of classical thermodynamics. Our results should motivate future experimental research in scaling up nanoscale thermal machines whereby nonadiabatic many-particle quantum effects are exploited to achieve quantum supremacy.
We consider a confined quantum fluid as a working medium with Hamiltonian where the inter-particle interactions are described by a potential that exhibits the scaling property V (λ z) = λ −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 ω(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], where E({m i }) is the eigenvalue of (A1) associated with the multi-index {m i |i = 1, . . . , N} defining a complete set of quantum numbers. Here, b = b(t) > 0 is the scaling factor that can be obtained as a solution of the Ermakov differential equation The boundary conditions b(0) = 1 andḃ(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, · · · , N) in [38,39] but follows from the exact coherent states under scaleinvariant 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 singleparticle 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 to derive the invariant of motion where T = T z T dil , and where H 0 denotes the Hamiltonian (A1) at t = 0. This invariant is an instance of the many-body invariants known under scale-invariant driving, see Eq. (56) in [42]. It has time-dependent eigenfunctions . It also means that the time-evolving state (A2) can be written as As a function of it, Hamiltonian (A1) can be rewritten as By (A7) we find that where we use which is derived from equation (A9). As a result, the mean nonadiabatic energy is given by During time-evolution, the mean-energy of an initial thermal state Collecting all the previous equations holding for the general Hamiltonian (A1), the nonadiabatic mean-energy following a change of ω(t) is found to be where p i is the momentum of the i-th 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 Indeed, it suffices to rewrite the squeezing operator as where we used the general form (A1). From (A14) it is clear that 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, Again, this formula is very general and holds for an initial thermal state. To derive (A17) we use the Hellmann-Feynman theorem Using the scaling invariance of the potential V , we have V ( Therefore, we find that ∂ H(0) /∂ ω 0 = ω −1 0 H(0) . Gathering equations (A15), (A16), and (A17), we obtain where b ad = [ω 0 /ω(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 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 as the mean energy in the adiabatic limit. As a result, Q * (t) is simply the ratio between the nonadiabatic and adiabatic mean energies.

Equivalence of the nonadiabatic factor under scale-invariant dynamics and the Husimi formula
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 * H (t) [25], 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 where G 1 (t) and G 2 (t) are two fundamental solution of the classical harmonic oscillator equationẍ(t) + ω(t) 2 x(t) = 0 that satisfy the following initial conditions: The Wronskian of these two functions 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), To prove (A25), we first compute the first time-derivative of b(t) as well as its second derivativë 3 .
The first term of the latter equation reduces to The second term is given byĠ while the third term gives 3 .
Thus, after adding these two terms, we find 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) and (A25)-(A26).
Appendix B: General results for the many-particle engine

Preliminary
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, . . . , N}, it is convenient to use the set of occupation numbers {n j | j = 1, . . . , ∞} in each renormalized harmonic j-mode of energy ε j = jhω, with the identification n j = ∑ N i=1 δ j,m i where δ 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 where {n j } is the occupation number in the j-th mode of energy jhω, satisfying the normalization ∑ ∞ j=1 n k = N, and E({n j }) = hω 2 N[1 + λ (N − 1)] + ∑ ∞ j=1 jhωn j . It follows that the equilibrium mean energy, for a given frequency and temperature (ω, β ), is given by which can be read as where the ground state energy E N (ω, β ) correspond to the first and second term of (B2), respectively. Depending on the value of σ = Nhβ ω, three different regimes can be distinguished: • At very low temperature σ N, the mean energy is approximately given by the ground state energy • At high temperature σ 1, we find corrections to the thermal energy due to quantum fluctuations where g N = (N − 1)N −1 and we have disregarded corrections of order O(σ 3 ). Formula (B5) is derived using the fact that for The sums can be carried out explicitly ∑ N k=1 k = N(N + 1)/2 and ∑ N k=1 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 < σ N (which requires N to be large compared to 1) where µ λ (σ ) = 1 To obtain the last equation, we estimated the Riemann sum in (B2) by a Riemann integral sincehβ ω 1 (the corrections are bounded by N(hβ ω) 2 /2 and so N can not be too large). The integral in (B7) can be found in closed form in terms of the standard polylogarithm function Li n (z) = ∑ ∞ j=1 j −n z j , n ≥ 1 [46]. The behavior of this function as a function of σ for different values of λ is shown in Fig. 4. 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, Nhβ ω 1, corresponds to the classical limit (quantum fluctuation are negligible in this regime) where the mean energy reads We note that to derive (B8) we use formula (B5) where we neglect the quantum corrections O(hβ ω).
Consider now a given protocol ω(t) for the compression stroke such that ω(0) = ω 1 , ω(τ) = ω 2 (similarly, for the compression stroke, letω(t) satisfyω(0) = ω 2 ,ω(τ) = ω 1 ). We denote the nonadiabatic factor Q * AB (τ) (resp. Q * CD (τ)) corresponding to the step 1 of the Otto cycle A → B (resp. the step 3 C → D), see formula (A22). Let us assume that for any τ > 0 the nonadiabatic factor is continuously differentiable with respect to ω 1 and ω 2 . We recall that for the adiabatic case, Q * AB = Q * CD = 1 while in the nonadiabatic sudden-quench limit Q * AB = Q * CD = ω 2 1 +ω 2 2 2ω 1 ω 2 = x 2 +1 2x with x = ω 1 /ω 2 . For the Otto cycle, the total work is usually defined as W = W 1 + W 3 , where W 1 > 0 is the work corresponding to the compression (ω 1 → ω 2 > ω 1 ) and W 3 < 0 is the work corresponding to the expansion (ω 2 → ω 1 < ω 2 ). The engine produces work when W is negative. Alternatively, we defineW We recall the equations given in the main body of the article: where H A = E N (ω 1 , β c ) and H C = E N (ω 2 , β h ) denote the energies at the equilibrium. As a consequence, the total work of the engine is given bỹ which is also also equal toW = Q 2 + H A − Q * CD ω 1 ω 2 H C . Then, the efficiency can be written in the following form In the adiabatic case (Q * = 1) the efficiency reduces to the Otto efficiency 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 ω 1 /ω 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) ≥ 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, λ 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 Since (Q * CD ) −1 ≤ 1 and Q * AB ≥ 1, we find the following upper bound , 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,

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 =W /τ, where τ is the running time of a cycle andW 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.

Effect of λ on the total work and on the efficiency
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 thatW where C N,λ = N(N + 1)λ andW 0 is the total work for λ = 0. Knowing that Q * AB ≥ 1 and Q * CD ≥ 1, it follows that ∂ λW λ ≤ 0. In addition, for the adiabatic case one findsW λ =W 0 (since Q * = 1) and the total work does not depends on λ . Similarly, one can show that Q 2 λ = Q 2 0 +¯h ω 2 2 (1 − Q * AB )C N,λ so ∂ λ Q 2 λ ≤ 0. In regard to the efficiency: The first derivative of the above quantity with respect to λ is given by It follows that for the adiabatic case (Q * AB = Q * CD = 1) the efficiency does not depends on λ , i.e., η λ = η 0 . Next we show that the optimal work is a monotonically increasing function of λ . For the sake of simplicity, we introduce the We proved that ∀x ∈ [0, 1], ∂ λW λ ≤ 0 which means that ∀x ∈ [0, 1] and ∀λ , λ such that λ ≤ λ ,W λ (x) ≥W λ (x), and so ∀x

Quantum Supremacy under identical resources
We introduce two parameters with a < 1, x < 1.
To quantify the boost in the performance of the QHE resulting from many-particle effects, we introduce the following quantities The first one is the the ratio between the total output work W N,λ (x, a) of a many-particle QHE and the corresponding value of N independent single-particle QHEs given by N ×W 1 (x, a) under the same resources, that is, for the same values of N, λ , β c,h , ω 1,2 . The ratio ρ N,λ (x, a) between the corresponding efficiencies is defined analogously. We note that the ratio between output power is equal to ρ N,λ (x, a) as the cycle time is the same for both systems.
Assuming that the temperature of the cold reservoir is highhω 1 β c 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 β c,h , trap frequency ω 1 , and interaction strength λ but with a frequency ratio x = ω 1 /ω 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 ω 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 considerhω 1 β c 1 and σ c , σ h > 1, meaning that N is large compared to 1. Then, from equations (B6), (B9), (B10), and (B11) we find and As a result, where with σ h = a x σ c < σ 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 = β h /β c in the single-particle quantities W 1 and η 1 . One finds that From equations (B20b) and (B19) we find that the efficiency ratio (B18b) exceeds unity if the function µ λ (σ ) defined by (B7) is decreasing. We know that it is monotonically decreasing for λ = 0 and monotonically increasing for λ ≥ 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 singleparticle QHEs in the classical regime. For 0 < λ < 1/2, the function µ λ (σ ) decreases if σ is less than a critical value that increases with λ , see Fig. 4. Consequently for σ ≤ σ c , µ λ (σ ) is a decreasing function of σ if the interaction strength is below a critical value λ c which is typically less than 0.2 (we obtained this estimate after numerical computations). Thus the ratio ρ N,λ (x, a) is greater than 1 for λ ≤ λ 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, a f λ ) and W 1 (x, a) is greater than 1 for λ = 0 and for 0 < λ ≤ λ c , where λ 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 σ λ (σ c ) which is less than 1 for λ ≤ λ c , making the analysis more complex. Yet, we find that for σ c not too large (typically σ c ≤ 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.
Appendix C: Total work and efficiency in the adiabatic limit Given that the nonadiabatic factor reduces to unity Q * = 1 under slow driving, using (B2) the adiabatic heat can explicitly written as Under nonadiabatic driving, Q * is generally greater than 1 so Q 2 > 0 obeys a stronger condition. We will discuss this condition for the sudden quench in the next section.
Using the asymptotic expansion of the thermal energy, one can easily find the asymptotic expressions for the optimal work and the corresponding efficiency in different regimes distinguished by the value of σ c : • For σ c N, E (th) ≈ 0 so E N,1/2 (ω 1 , β c ) ≈¯h ω 1 4 N(N + 1) and then α ad ≈¯h ω 1 β h (N+1)) 4 which gives • For N > 1, we consider an intermediate regime corresponding to a large temperature σ c N but a relatively small temperature per particle σ c ≥ 1, which means that the particle number is large, N 1 but keeping σ h 1 (i.e., β 2 β 1 ). Using (B6) we find that E N,1/2 (ω 1 , β c ) ≈ N β c µ 1/2 (σ c ) where σ c = Nhβ c ω 1 , and µ 1/2 (σ c ) is defined by (B7). As a result, by (C4) and (C5),W FIG. 5. Performance of an adiabatic many-particle QHE compared to N single-particle QHE in series. The adiabatic efficiency at optimal output power of a N-particle QHE normalized by the single-particle value (a), and the adiabatic optimal output power divided by N times the single-particle value (b) (see equation (B15)) are displayed for N = 20, 200, 500 and 2, 000 from top to bottom, under an adiabatic driving with β c = 0.01/(hω 1 ). In the adiabatic limit, the efficiency becomes independent of the interaction strength λ and many-particle quantum effects are detrimental.
where µ 1/2 (σ c ) = µ(σ c ) is greater than unity and increases monotonically with σ c , as shown in Fig. 4. Consequently, in this regime the efficiency is lesser than the Curzon-Ahlborn efficiency 1 − √ a associated with a classical heat engine operated in the adiabatic limit [37], see Fig. 5 as well as Fig. 2. We emphasize that equations (C8a) and (C8b) are valid only for a 1.
Appendix D: Total work and efficiency for the sudden quench

Preliminary
For the sudden-quench driving between ω 1 and ω 2 the nonadiabatic factor Q * takes the same value in both the compression and expansion strokes and is given by We recall the condition over a = β h /β c and x = ω 1 /ω 2 (C2) which can be read as a ≤ x. One can derive another condition which is more useful for small x: To derive this condition we use (B9b) to show that We must have Q 2 ≥ 0, and the sum in the previous equation must be positive since the first term is negative in the adiabatic case (Q * ≥ 1). Using x ≤ e x − 1 ≤ 2x for x ≤ 1, one can easily show that an upper bound to the sum is positive if β c ω 1 β h ω 2 ≥ 2Q * which gives (D2) for small x. This condition is weak in the sense that it does not ensure that Q 2 is positive. However, we know that if this condition is not satisfied then Q 2 is negative.

Optimal output work and efficiency
We want to find a condition on x so that the work is maximal. To this end, we solve this equation 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 σ h = Nhβ h ω 2 1, the total work is given bỹ where we have kept only the first order correction term. If we assume that σ h 1 and a = β h /β c small, then the optimal value for the frequency ratio is given by Thus, one obtains the following expression for the optimal output work and efficiencỹ where we have introduced α sq = aβ c H A N and performed a very similar computation for the heat Q 2 . Explicit expressions for the optimal work at different regimes of inverse temperature β c can be found: • For σ c N, one has x sq ≈ κ N,λ hβ h ω 1 , where κ N,λ = (1 + (N − 1)λ ) for N ≥ 1. Then, • For large temperature σ c 1, for N = 1, we have x sq ≈ a 1/4 (1 + 1 48 σ 2 c ), and • For N > 1, we consider an intermediate regime corresponding to a large temperaturehβ h ω 1 1 but a relatively small temperature per particle σ c = Nhβ c ω 1 ≥ 1, which means that the number of particles is large σ h 1 (i.e., β 2 β 1 ). Using (B6), we find aβ c H A N ≈ aµ λ (σ c ) where µ λ (σ c ) is given by (B7). Notice that µ 0 (σ c ) is a decreasing function of σ c FIG. 6. Sudden quench and adiabatic efficiency at optimal work of a many-particle engine. The efficiency at optimal work as a function of β h /β c for β c = 0.01/(hω 1 ) is plotted. The two top curves represent the adiabatic efficiency for N = 1 (dotted line) and for N = 500 (continuous line). The four curves below correspond to the efficiency for a sudden quench driving for N = 1 (dotted line) and for N = 500 (continuous line for λ = 0, dotted-dashed line for λ = 1/2 and dashed line for λ = 1). We observe that the efficiency for an adiabatic driving is above the sudden quench efficiency as expected. The novelty concerns the variation of the efficiency with N and λ . It is clear that for λ < 1/2 (between bosons and semions) the efficiency of the many-particle QHE is enhanced while it decreases for λ ≥ 1/2 (semions, fermions and strongly correlated bosons). In the adiabatic case, the efficiency at optimal work does not depend on the exclusion statistics (i.e., on λ ) and decreases as a function of the particle number N, see Fig. 2 in the main body of the manuscript.
FIG. 7. Efficiency at optimal work under sudden-quench driving for a many-particle QHE and a working medium for interaction strength larger than 1. The efficiency is displayed for N = 200, β c = 0.01/(hω 1 ), and for λ = 1 (continuous line), λ = 2 (dotted-dashed line), λ = 4 (dashed line), λ = 8 (dotted line). This figure illustrates the general result we proved concerning the fact that the efficiency is maximal for bosons and decreases for λ > 0. For λ > 1 the Calogero-Sutherland gas describes strongly-correlated bosons, similar to a super-Tonks-Girardeau gas [47]. For large λ the working medium becomes a Wigner crystal with spatially separated particles. In this case it is clear that the many-particle QHE is not efficient. with the following bounds π 2 /(6Nhω 1 β c ) ≤ µ(σ c ) ≤ 1. If one considers σ c not too large (so that σ h = Nhβ h ω 2 1 still holds), by (D6) and (D7b) one findsW where µ λ (σ c ) is defined by (B7). Thus, the optimal output power as well as the efficiency at optimal output power can be substantially improved for λ < 1/2 (and for σ c small enough) as µ λ ≤ 1 and decreases otherwise (i.e., for λ > 1/2), see Figs. 4, 5, 6, and 7. Notice that (D11a) and (D11b) are valid only for a 1 as σ h 1. The case σ h ≥ 1 will be studied in the next section.
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 ≥ 1000. Numerically (see Fig. 2) we show that the efficiency can be enhanced up to 50%, see Fig. 5, and Fig. 8. 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 λ , ≈ π 2 n 12β h a 2 x 2 − 1 x 2 1 + 3λ σ 2 c π 2 + x(1 − x 2 ) 1 + 3a 2 λ σ 2 c π 2 x 2 , for 1 σ c,h N , FIG. 8. Efficiency at optimal work of a QHE with an ideal Bose gas as working medium. In this figure we show the lower and upper bound for the efficiency at optimal work in the large temperature regime 1 ≤ σ c(h) N. The top continuous line is the numerical exact efficiency for N = 1000, and the bottom continuous line is for N = 1. The dotted line (respectively dashed line) is the lower bound (respectively upper bound) given by (E2). We remark that the efficiency can be improved drastically and a similar plot for the relative efficiency shows that it is enhanced between 25% (lower bound) and 50% (upper bound). We notice that for β h /β c 0.4 the upper bound is sharper than the lower bound and provides an accurate estimate of the efficiency.
where n = (hβ h ω 2 ) −1 . To derive equation (E3b) we use (B6). In the regime of large number of particles (σ c(h) 1), we estimate µ λ as Hence, one finds different regimes depending on the value of λ . For λ λ c = π 2 /(3σ 2 c ) the thermal energy dominates and then the problem can be treated as the case with λ = 0. On the contrary, for λ λ c the thermal energy becomes negligible, thus the quantum fluid works less efficiently and becomes "frozen" for λ = 1/2 and 1, see Figs. 5, 6 and 7.
To conclude, for a many-particle QHE in the low temperature regime nonadiabatic effects can enhance the efficiency over the single-particle value for λ = 0 and for small λ (typically λ λ c ). However, for λ ≥ 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.

Adiabatic driving
Under adiabatic driving, the total work is given bỹ 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 whose solution depends only on the temperature ratio a. Here we take a to be greater than 0.1 so that N (hβ h ω 2 ) −1 . One can find an analytic expression of the real solution of this polynomial equation and provide an estimate for a ≥ 0.1, Inserting this value in the Otto efficiency we obtain η (N) ad ≈ 9 − 10a + a 2 16 . (E5)