An Interaction-Driven Many-Particle Quantum Heat Engine: Universal Behavior

A quantum heat engine (QHE) based on the interaction driving of a many-particle working medium is introduced. The cycle alternates isochoric heating and cooling strokes with both interaction-driven processes that are simultaneously isochoric and isentropic. When the working substance is confined in a tight waveguide, the efficiency of the cycle becomes universal at low temperatures and governed by the ratio of velocities of a Luttinger liquid. We demonstrate the performance of the engine with an interacting Bose gas as a working medium and show that the average work per particle is maximum at criticality. We further discuss a work outcoupling mechanism based on the dependence of the interaction strength on the external spin degrees of freedom.

A quantum heat engine (QHE) based on the interaction driving of a many-particle working medium is introduced. The cycle alternates isochoric heating and cooling strokes with both interaction-driven processes that are simultaneously isochoric and isentropic. When the working substance is confined in a tight waveguide, the efficiency of the cycle becomes universal at low temperatures and governed by the ratio of velocities of a Luttinger liquid. We demonstrate the performance of the engine with an interacting Bose gas as a working medium and show that the average work per particle is maximum at criticality. We further discuss a work outcoupling mechanism based on the dependence of the interaction strength on the external spin degrees of freedom. Universality plays a crucial role in thermodynamics, as emphasized by the description of heat engines, that transform heat and other resources into work. In paradigmatic cycles (Carnot, Otto, etc.) the role of the working substance is secondary. The emphasis on universality thus seems to prevent alternative protocols that exploit the many-particle nature of the working substance.
In the quantum domain, this state of affairs should be revisited as suggested by recent works focused on many-particle quantum thermodynamics. The quantum statistics of the working substance can substantially affect the performance of a quantum engine [1]. The need to consider a many-particle thermodynamic cycle arises naturally in an effort to scale-up thermodynamic devices [2][3][4][5] and has prompted the identification of optimal confining potentials [6] as well as the design of superadiabatic protocols [7][8][9][10], first proposed in a singleparticle setting [11,12]. The divergence of energy fluctuations in a working substance near a second-order phase transition has also been proposed to engineer critical Otto engines [13]. In addition, many-particle thermodynamics can lead to quantum supremacy, whereby quantum effects boost the efficiency of a cycle beyond the classical achievable bound [7,14]. The realization of superadiabatic strokes with ultracold atoms using a Fermi gas as a working substance has been reported in [9,10,15].
Quantum technologies have also uncovered novel avenues to design thermodynamic cycles. Traditionally, interactions between particles are generally considered to be "fixed by Nature" in condensed matter. However, a variety of techniques allow to modify interparticle interactions in different quantum platforms. A paradigmatic examples is the use of Feschbach and confinement-induced resonances in ultracold atoms [16]. Digital quantum simulation similarly allows to engineer interactions in trapped ions and superconducting qubits [17].
In this Letter, we introduce a novel thermodynamic cycle that exploits the many-particle nature of the working substance: it consists of four isochoric strokes, alternating heating and cooling with a modulation of the interparticle interactions. This cycle resembles the four-stroke quantum Otto engine in that expansion and compression strokes are substituted by isochoric processes in which the interparticle interactions are modulated in time. Using Luttinger liquid theory, the efficiency of the cycle is shown to be universal in the lowtemperature regime of a one-dimensional (1D) working substance. We further shown that when an interacting Bose gas is used as such, the average work output is maximum at criticality.
Interaction-driven thermodynamic cycle.-We consider a quantum heat engine (QHE) with a working substance consisting of a low-dimensional ultracold gas tightly confined in a waveguide. Ultracold gases have been previously considered based quantum cycles where work is done via expansion and compression processes, both in the non-interacting [6,22] and interacting regimes [7,[23][24][25][26][27]. We proposed the implementation of quantum cycle consisting of four isochoric strokes, in which heating and cooling strokes are alternated with isentropic interaction-driven processes. In the latter, work is done onto and by the working substance by increasing and decreasing the interatomic interaction strength, respectively. This work can be transferred to other degrees of freedom as we shall discuss below. The working substance consists of N particles with interparticle interactions parameterized by the interaction strength c. Both the particle number N and system size L are preserved throughout the cycle and any equilibrium point is parameterized by a point (c, T ) indexed by the temperature T and interaction strength c. Specifically the interactiondriven quantum cycle, shown in Fig. 1 for the 1D Lieb-Linger gas [18,19], involves the following strokes: Interaction driven quantum cycle. The working substance is driven through four sequential isochoric strokes alternating heating and cooling processes at different temperature T with isentropes in which the interparticle interaction strength c is ramped up and down. The dashed lines correspond to the thermal entropy calculated from the thermodynamic Bethe ansatz equation of the Lieb-Liniger gas [18][19][20][21].
(1) Interaction ramp-up isentrope (A → B): The working substance is initially in the thermal state A parameterized by (c A , T A ) and decoupled from any heat reservoir. Under unitary evolution, the interaction strength is enhanced to the value c B and the final state is non-thermal. (2) Hot isochore (B → C): Keeping c B constant, the working substance is put in contact with the hot reservoir at temperature T B and reaches the equilibrium state (c B , T C ).
Interacting Bose gas as a working substance.-Consider as a working substance an ultracold interacting Bose gas tightly confined in a waveguide [28,29], as realized in the laboratory [30][31][32], see Fig. 1. The effective Hamiltonian for N particles is that of the Lieb-Liniger model [18,19] where x jℓ = x j − x ℓ , with 2m =h = 1. The spectral properties of Hamiltonian (2) can be found using coordinate Bethe ansatz [18,19]. We consider a box-like trap [33,34] where any energy eigenvalue can be written as E = ∑ i k 2 i in terms of the ordered the quasimomenta 0 < k 1 < k 2 < · · · < k N . The latter are the (Bethe) roots {k i } of the coupled algebraic equations determined by the sequence of quantum numbers {I i } with i = 1, 2, · · · , N. As a function of c and T , the 1D Bose gas exhibits a rich phase diagram. We first consider the strongly coupling regime and use a Taylor series expansion in 1/c. For a given set of quantum numbers I n = {I (n) i }, the corresponding energy eigenvalue takes the form where the interaction-dependent factor λ c reads The spectrum of a strongly interacting Bose gas is thus characterized by eigenvalues with scale-invariant behavior, i.e., ε n (c)/ε n (c ′ ) = λ c /λ ′ c . In this regime, the work output is thus becomes independent of temperature of the heat reservoirs.
Beyond the strongly-interacting and low-energy regimes, we resort to a numerically-exact solution of Eqs. (3) for finite particle number N. We enumerate all the possible sets I n of quantum numbers for low-energy states and solve Eqs. (3) numerically for a given interaction c. With the resulting quasimomenta {k n,1 , k n,2 , · · · }, and the corresponding energy eigenvalues ε n = ∑ i k 2 n,i , the probability that the Bose gas at temperature T is found with energy ε n is set by the Boltzmann weights p n = e −ε n /T / ∑ m e −ε m /T (with k B = 1). The equilibrium energy of the states A and C is set by thermal averages of the form E = ∑ n p n ε n that in turn yield the expressions for Q 2 and Q 4 . Here, a proper cutoff of the possible sets {I n } can be determined by p n /p G ≪ 1, where p G is the probability for the working substance to be in the ground state. The numerical results for the efficiency η and output work W are shown in Fig. 2 as a function of the interaction strength. For fixed values of T C and T A , the maximum work is studied as a function of c A while keeping c B constant, see Fig. 2. The efficiency is well reproduced by Eq. (6) at strong coupling, that captures as well a monotonic decay with increasing interaction strength. The efficiency is found to be essentially independent of the temperature in the strong interaction regime, whereas the work output is governed by the temperature and interaction strength.
In the thermodynamic limit (where N and L → ∞ with n = N/L being kept constant), the equilibrium state of the 1D Bose gas is determined by the Yang-Yang thermodynamics [20, 21, 35 -37]. The pressure is then given by where the "dressed energy" ε(k) is determined by thermodynamic Bethe Anstatz (TBA) equation The particle density n and entropy density s can be derived from the thermodynamics relations in terms of which the internal energy density reads E = −p + µn + T s. Both interaction-driven strokes are considered to be adiabatic. As a result, the heat absorbed during the hot isochore stroke (B → C) is given by Here T B and T D can be determined from the entropy, by setting s(c A , where T A and T C are the temperature of the cold and hot reservoir, respectively. The efficiency and work can then be obtained by numerically solving the TBA equation. Universal efficiency at low temperature.-The low-energy behaviour of 1D Bose gases is described by the Tomonaga-Luttinger liquid (TLL) theory [38][39][40][41][42], in which the free energy density reads where E 0 is the energy density of ground state, v s is the sound velocity which depends on particle density n and interaction c. The entropy density s can be obtained as the derivative of free energy, s = −∂ F /∂ T = πT /3v s . The expression for the heat absorbed and released are respectively given by where are the entropy density and the sound velocity of the state i, respectively. Using the fact that the strokes (1) and (3) are isentropes, it follows that and as a result, the efficiency and work output are given by As the TLL theory describes the universal low-energy behavior of 1D many-body systems, Eqs. (14)-(16) provide a universal description of the efficiency and work of quantum heat engines with a 1D interacting working substance at low temperatures, which are applicable to any cycles whose working strokes and heat exchanging strokes are separated. In particular, in the strongly interacting regime, the sound velocity of 1D Bose gases is given by v s ≃ 2πn(1 − 4nc −1 + 12n 2 c −2 ) [21]. In this regime, thus the result (14) reduces to Eq. (6). On the other hand, in the weak interaction regime, the sound velocity is given by v s ≃ 2n c n − 1 2π c n 3 2 [21], and thus the efficiency yields indicating an enhancement of the performance with respect to the strongly-interacting case. Quantum critical region.-We next focus on the performance of an interaction-driven QHE at quantum criticality [44]. The 1D Bose gas displays a rich critical behavior in the temperature T and chemical potential µ plane, see Fig. 3. In the region of µ ≪ 0 and T ≫ n 2 , i.e. the mean distance between atoms is much larger than thermal wave length, the system behaves as a classical gas (CG). A quantum critical region (QC) emerges between two critical temperatures (see the white dashed lines in have the same power of the temperature dependence. In the region with µ > 0 and temperatures below the right critical temperature, the quantum and thermal fluctuations can reach equal footing: the TLL region discussed above. For fixed cycle parameters (c A , c B , T A , and T C ), we study the performance of the engine across the QC region by changing n. We numerically calculate the efficiency η and average work W /N by using the TBA equation (8), and show that near the quantum critical region W /N has a maximum value, see Fig. 3. We set c A = 1, c B = 3, T A = 1, and T C = 5 for the heat engine and let the density n increase from 0.1 to 23. The red, green, and blue dashed lines in Fig. 3 correspond to the density n ≃ 0.2, 1.4, and 6.2, respectively. In order to understand the maximum of W /N, we also plot this three engine cycles (A → B → C → D → A) in the phase diagram of specific heat in the T − µ plane in the right panel of Fig. 3. When the engine works near the boundary from QC to TLL, it has the maximum average work.
Discussion.-The modulation of the interaction strength c is associated with the performance of work, that has recently been investigated for different working substances [45][46][47]. An experimentally-realizable work outcoupling mechanism can be engineered whenever the value of the coupling strength specifically depends on the configuration of other degrees of freedom. This is analogous to the standard Carnot or Otto cycles where the working substance is confined in a box-like, harmonic, or more general potential [48]. The latter is endowed with a dynamical degree of freedom that is assumed to be slow (massive) so that in the spirit of the Born-Oppenheimer approximation can be replaced by a parameter. Similarly, the modulation of the coupling strength in an interaction-driven cycle can be associated with the coupling to external degrees of freedom.
As an instance, consider the choice of an interacting SU(2) 1D spinor Fermi gas in a tight waveguide as a working substance. The Hamiltonian can be mapped to the Lieb-Liniger model [49] with an operator-valued coupling strength that depends on the spinor degrees of freedomĉ =ĉ(Ŝ j ·Ŝ ℓ ). The latter can be reduced to a real number c = c( Ŝ j ·Ŝ ℓ ) making use of a variational method [49,50], an approximation corroborated by the exact solution in a broad range of parameters [51]. Thus, the dependence of the interaction strength on the spin degrees of freedom provides a possible work outcoupling mechanism. An alternative relies on the use of confinementinduced resonances [16,52]. The scattering properties of a tightly confined quasi-1D working substance in a waveguide can be tuned by changing the transverse harmonic confinement of the waveguide. The frequency ω ⊥ of the latter directly determines the interaction strength c, i.e., c = c(ω ⊥ ). In this case, the role of ω ⊥ parallels that of the box-size or harmonic frequency in the conventional Carnot and Otto cycles with a confined working substance.
Finally, we note that an interaction-driven cycle can also be used to describe QHEs in which the interaction-driven strokes are substituted by processes involving the transmutation of the particle quantum exchange statistics, e.g., a change of the statistical parameter of the working substance. The Hamiltonian (2) can be used to describe 1D anyons with pair-wise contact interactions with coupling strengthc and statistical parameter θ characterizing the exchange statistics, smoothly interpolating between bosons and fermions [53,54]. The spectral properties of this Lieb-Liniger anyons can be mapped to a bosonic Lieb-Liniger model (2) with coupling strength c =c/ cos(θ /2) [53,54]. The modulation of c can be achieved by the control of the particle statistics, tuning θ as proposed in [55].
Conclusions.-We have proposed an experimental realization of an interaction-driven quantum heat engine that has no single-particle counterpart: It is based on a novel quantum cycle that alternates heating and cooling strokes with processes that are both isochoric and isentropic and in which work is done onto or by the working substance by changing in the in-teratomic interaction strength. This cycle can be realized with a Bose gas in a tight-waveguide as a working substance. Using Luttinger liquid theory, the engine efficiency has been shown to be universal in the low temperature limit, and set by the ratio of the sound velocities in the interaction-driven strokes. The optimal work can be achieved by changing the ratio of the sound velocity, e.g., by tuning the interaction strength. An analysis of the engine performance across the phase diagram of the Bose gas indicates that quantum criticality maximizes the efficiency of the cycle.
Our proposal can be extended to Carnot-like interaction driven cycles in which work and heat are simultaneously exchanged in each stroke. Exploiting effects beyond adiabatic limit may lead to a quantum-enhanced performance [7,14]. The use of non-thermal reservoirs [56][57][58] and quantum measurements [59,60] constitutes another interesting prospect. Our results identify confined Bose gases as an ideal platform for the engineering of scalable many-particle quantum thermodynamic devices.

Consider a system with N identical bosons with a contact interaction confined in a one-dimensional (1D) hard-wall potential. This system is described by the Hamiltonian
where g 1D =h 2 c/m is the 1D coupling constant and c parametrizes the interaction strength. Hereafter, we seth = 2m = 1. This model can be exactly solved by the Bethe ansatz [33,34], and the wavefunction is given by where the sum is taken over all N! permutations P of N integers such that P : (1, 2, · · · , N) → (P 1 , P 2 , · · · P N ), and ε i takes ±1.
Here, k P i is the the wave number which satisfies the Bethe equations and the energy of this state is given by Taking the logarithm of the Bethe equations (20), one obtains where {I i } are integer quantum numbers, which are arranged in ascending order: 1 ≤ I 1 ≤ · · · ≤ I N . Thus, the wave numbers {k i } are also in ascending order, i.e., 0 < k 1 < · · · < k N .

Scaling invariant behavior
We next focus on the low energy excitations which dominate the thermodynamics in the strongly-interacting Bose gas at low temperature. A Taylor expansion of the rhs of Eq. (22) for c ≫ k N yields the following asymptotic solution Note that, in this regime, the dependence of the quasimomentum k i on the other quasimomenta k j with j = i is of higher order in k N /c. In this sense, the algebraic Bethe ansatz equations (3) decouple in this limit. The energy eigenvalue ε n (c) for a given set of quantum numbers I n ≡ {I Thus, up to corrections of order 1/c 2 , the energy eigenvalue can be written in terms of that of free fermions rescaled by an overall factor resulting from the interactions. We refer to this factor as the generalized exclusion statistics parameter, in view of [37], and define it as Using it, the energy eigenvalues of a strongly-interacting 1D Bose gases simply read and exhibit the scale-invariant behavior

ENERGIES OF EQUILIBRIUM AND NON-EQUILIBRIUM STATES
The equilibrium state of the system at temperature T can be characterized by the partition function where the summation is taken over all the quantum states. The average energy of the system is given by with p n (c, T ) being the probability measure of the n-th eigenstate in the Gibbs ensemble given by During the interaction-driven strokes, the system generally deviates from the initial equilibrium state as a result of the modulation of the interaction strength c. However, since the energy gap for low-energy excitations is nonzero in the system with finite N bosons, a sufficiently slow ramping process becomes adiabatic and keeps the probability distribution p n (c, T ) unchanged. Therefore, for the non-equilibrium state resulting from an adiabatic ramping process, the average energy of the system at the interaction parameter c ′ is given by Due to the scale-invariant behavior at strong coupling, the energy of the non-equilibrium state at the end of the interaction-ramp isentropic process is found to be Furthermore, the scale-invariant behavior Eq. (27) indicates that we can relate a non-equilibrium state in adiabatic interactiondriven processes to an equilibrium state with an effective temperature T ′ by noting that where the effective temperature reads

QUANTUM HEAT ENGINE
For the interaction-driven quantum heat engine (see Fig. 1 in the main text), the heat absorbed from the hot reservoir is given by while the heat released from the engine to the cold reservoir reads Therefore, the efficiency η and the work W done by the engine are given by The efficiency η and work W can also be obtained numerically by solving Eq. (3). In our calculations, we take a proper cutoff for excited states whose probability p n given by Eq. (30) is much smaller than that of the ground state.

THERMODYNAMIC LIMIT
In the laboratory, a Bose gas confined in an effectively 1D trap typically consists of thousands of particles. It is generally considered that such system is well described by the thermodynamic limit, i.e., N and L → ∞ with n = N/L being kept constant. The equilibrium states can then be described by the Yang-Yang thermodynamic equation [20]. The pressure p, particle density n, entropy density s and internal energy density E can then be found as described in the text.
Note that we take the thermodynamic limit after the slow ramping limit to guarantee that no diabatic transitions occur during the ramping processes. Thus, the probability p n (c, T ) is unchanged during the ramping processes. In these processes, the entropy should also be constant because no heat is transferred in or out of the working substance. Given the expressions for the heat absorbed during the hot isochore stroke (B to C) and the heat released during the cold isochore stroke (D to A) the efficiency η and work output W are given by Here, T B and T D can be determined using the fact that the interaction-driven strokes are isentropic, namely where T A and T C denote the temperature of the cold and hot reservoir, respectively. The efficiency can thus be obtained numerically by solving Eq. (8).

Universality at low energies
The low-energy physics of 1D Bose gases can be well captured by the Luttinger liquid theory. The free energy density F is given by [61] F where E 0 is the energy density of the ground state, v s is the sound velocity which depends on the particle density n and the interaction strength c. The entropy density s is given by the derivative of the free energy, namely, The heat absorbed from the hot resevoir in the hot isochore stroke (B to C) is given by and the heat released to the cold reservoir in the cold isochore stroke (D to A) is Here, namely, To simplify, we introduce the following two dimensionless parameters: Since T A < T B < T C , we get The work can be extracted from the heat engine is given by An interesting question is what is the optimal work extracted from this heat engine when the temperatures of the two reservoirs, T A and T C , are fixed. It is easy to show that the work output W always has a maximum value (see the right panel of Fig. 4) at ξ = ξ c with  52) and (55), respectively, where ξ for a given c A is obtained by numerically calculating the sound velocities v A s and v B s at zero temperature [21]. The black dashed line corresponds to ξ c ≈ 0.69 for κ = 0.5 given by Eq. (54). In our numerical calculation, we set n = 1, c B = 2, and L = 1.
where a = 27κ 2 1 + (κ 2 /27) − 1 ≈ κ 4 /2. Thus for small κ ≪ 1, we get The efficiency of the heat engine is Especially, for the strong interaction case, the sound velocity is v s ≈ 2πn 1 − 4(n/c) + 12(n/c) 2 [21], and thus the efficiency is given by which agrees with Eq. (37) in the thermodynamic limit. For the weak interaction case, the sound velocity is v s ≈ 2n (c/n) − (2π) −1 (c/n) 3/2 1/2 [21], and the efficiency is given by The work output and the efficiency given by Eqs. (52) and (55) are universal at low energies for 1D system. Here we shall take 1D Bose gases as a platform to test these universal properties. For the interaction-driven engine with fixed particle number N and the length L studied in the present work, the sound velocity is changed during the two isoentropic strokes [(A to B) and (C to D)] by changing the interaction strength. We numerically calculate the work W and efficiency η for κ = 0.5 and n = N/L = 1. The interaction strength c B is fixed at c B = 2, which corresponds to the sound velocity v B s ≈ 2.50. From Eq. (54), the optimal work is obtained at ξ c ≈ 0.69 (58) for κ = 0.5. The numerical results are shown by the blue lines in Fig. 4. In the low temperature region, the numerical results are well explained by the results of the Luttinger liquid theory given by Eqs. (52) and (55), which are shown by the red solid lines.
In the high temperature region, these two results deviate, which indicates the breakdown of the Luttinger liquid theory. FIG. 5. Phase diagram of the 1D Bose gas. The color contour shows the specific heat. With increasing the chemical potential µ, the system undergoes a crossover from the classical gas (CG) to the Tomonaga-Luttinger liquid (TLL) across the quantum critical (QC) region.

Quantum criticality
The 1D Bose gases show rich critical properties in the plane of the temperature T and the chemical potential µ; see Fig. 5. In the so-called the Tomonaga-Luttinger liquid (TLL) region, where µ > 0 and the temperature is below the right critical temperature, the magnitude of quantum and thermal fluctuations are comparable. In this region, the low-energy properties can be well captured by the Luttinger liquid theory and the excitation spectrum is given by E = v s |k|, where v s and k are the sound velocity and the wave vector, respectively. In the region of µ/T ≪ 0, i.e., the mean distance between atoms is much larger than thermal wave length, the system behaves as a classical gas (CG). A quantum critical (QC) region emerges between two critical temperatures (see the white dashed lines in Fig. 5) fanning out from the critical point µ c = 0. In this region, quantum fluctuation and thermal fluctuation have the same power of the temperature dependence. The density n is monotonically increasing with the chemical potential µ for a given temperature T and the interaction strength c. Therefore, one can go across the QC region from the CG to TLL regions only by increasing the density n of the working substance. Our numerical calculation shows W /N takes a maximum value in the crossover region between QC and TLL (see Fig. 3 in the main text).

1D SPINOR FERMI GAS AS WORKING SUBSTANCE
Denoting byP s i j = 1 4 −Ŝ i ·Ŝ j andP t i j = 3 4 +Ŝ i ·Ŝ j are the projectors onto the subspaces of singlet and triplet functions of the spin arguments (σ i , σ j ) for fixed values of all other arguments, the Hamiltonian of the system is given by [49,50] Here, v o is a strong, attractive, zero-range, and odd-wave interaction that is the 1D analog of 3D p-wave interaction. Similarly, g e denotes the even-wave 1D coupling constant arising from 3D s-wave scattering [28]. The Hamiltonian of the spinor Fermi gas can be mapped to the Lieb-Liniger Hamiltonian [49] by promoting the coupling strength to an operator dependent on the spin degrees of freedomĉ where the coupling constants {c o , c e } are set by g e and v o . The operator-valued coupling strengthĉ can however be reduced to a real number c = (3c o + c e )/4 + (c o − c e ) Ŝ i ·Ŝ j making use of a variational method that combines the exact solution of the LL model with that of the 1D Heisenberg models [49,50]. As it turns out, this approximation is corroborated by the exact solution of (59) in a broad range of parameters [51]. Thus, the dependence of the interaction strength on the spin degrees of freedom provides a possible work outcoupling mechanism.