Interaction-driven thermodynamic cycle

Consider a quantum heat engine (QHE) with a working substance consisting of a low-dimensional ultracold gas tightly confined in a waveguide.25,26 Ultracold gases have previously been explored in quantum cycles where work is done via expansion and compression processes, both in the non-interacting7,27 and interacting regimes.8,28,29,30,31,32 We propose the implementation of a 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 the system size \(L\) are kept constant throughout the cycle and any equilibrium point is parameterized by a point \((c,T)\) indexed by the temperature \(T\) and the interaction strength \(c\). Specifically the interaction-driven quantum cycle, shown in Fig. 1 for the 1D Lieb-Linger gas,21,22 involves the following strokes:

(1) Interaction ramp-up isentrope (A \(\to\) B): The working substance is initially in the thermal state \(A\) parameterized by \(({c}_{{\rm{A}}},{T}_{{\rm{A}}})\) and decoupled from any heat reservoir. Under unitary evolution, the interaction strength is enhanced to the value \({c}_{{\rm{B}}}\) and the final state is non-thermal. (2) Hot isochore (B \(\to\) C): Keeping \({c}_{{\rm{B}}}\) constant, the working substance is put in contact with the hot reservoir at temperature \({T}_{{\rm{C}}}\) and reaches the equilibrium state \(({c}_{{\rm{B}}},{T}_{{\rm{C}}})\). (3) Interaction ramp-down isentrope (C \(\to\) D): The working substance is decoupled from the hot reservoir and performs work adiabatically while the interaction strength decreases from \({c}_{{\rm{B}}}\) to \({c}_{{\rm{A}}}\), reaching a non-thermal state. (4) Cold isochore (D \(\to\) A): The working substance is put in contact with the cold reservoir with temperature \({T}_{{\rm{A}}}\) keeping the interaction strength constant until it reaches the thermal state \(({c}_{{\rm{A}}},{T}_{{\rm{A}}})\). The work extracted from the heat engine is given by \(W={W}_{3}-{W}_{1}={Q}_{2}-{Q}_{4}\) while the efficiency of the heat engine reads

$$\eta =\frac{W}{{Q}_{2}}=1-\frac{{Q}_{4}}{{Q}_{2}}.$$ (1)

Interacting Bose gas as a working substance

Consider as a working substance an ultracold interacting Bose gas tightly confined in a waveguide,33,34 as realized in the laboratory.35,36,37 The effective Hamiltonian for \(N\) particles is that of the Lieb-Liniger model:21,22

$${\hat{H}}_{\text{LL}}\,=\,-{\sum }_{j=1}^{N}{\partial }_{{x}_{j}}^{2}+{\sum }_{1\le j<\ell \le N}2c\delta ({x}_{j\ell }),$$ (2)

where \({x}_{j\ell }={x}_{j}-{x}_{\ell }\), with \(2m=\hslash =1\). The spectral properties of the Hamiltonian of Eq. 2 can be found using coordinate Bethe ansatz.21,22 We consider a box-like trap38,39 where any energy eigenvalue can be written as \(E={\sum }_{i}{k}_{i}^{2}\) in terms of the ordered quasimomenta \(0<{k}_{1}<{k}_{2}<\cdots <{k}_{N}\). The latter are the (Bethe) roots \(\{{k}_{i}\}\) of the coupled algebraic equations

$$L{k}_{i}=\pi {I}_{i}-{\sum }_{j

e i}\left(\arctan \frac{{k}_{i}-{k}_{j}}{c}+\arctan \frac{{k}_{i}+{k}_{j}}{c}\right),$$ (3)

determined by the sequence of quantum numbers \(\{{I}_{i}\}\) with \(i=1,\ 2,\ \cdots \ ,\ N\). As a function of \(c\) and \(T\), the 1D Bose gas exhibits a rich phase diagram. We first consider the regime of strong interactions and use a Taylor series expansion in \(1/c\). For a given set of quantum numbers \({{\mathcal{I}}}_{n}=\{{I}_{i}^{(n)}\}\), the corresponding energy eigenvalue takes the form

$${\epsilon }_{n}(c)\approx \frac{{\pi }^{2}{\lambda }_{c}}{{L}^{2}}{\sum }_{i=1}^{N}{{I}_{i}^{(n)}}^{2},$$ (4)

where the interaction-dependent factor \({\lambda }_{c}\) reads

$${\lambda }_{c}=1-\frac{4(N-1)}{cL}+\frac{12{(N-1)}^{2}}{{c}^{2}{L}^{2}}.$$ (5)

The spectrum of a strongly interacting Bose gas is thus characterized by eigenvalues with scale-invariant behavior, i.e., \({\epsilon }_{n}(c)/{\epsilon }_{n}(c^{\prime} )={\lambda }_{c}/\lambda_{c} ^{\prime}\). This scale invariance in the spectrum allows us to assign an effective temperature along the isentropes; see Supplementary Information Section I. In this regime, the work output is thus set by \(W={Q}_{2}-{Q}_{4}=[1-({\lambda }_{{c}_{{\rm{A}}}}/{\lambda }_{{c}_{{\rm{B}}}})]{Q}_{2}\) and the efficiency

$$\eta =1-\frac{{\lambda }_{{c}_{{\rm{A}}}}}{{\lambda }_{{c}_{{\rm{B}}}}}$$ (6)

becomes independent of temperature of the heat reservoirs, as we show in Supplementary Information. Being the spectrum scale-invariant, the efficiency can be expressed in terms of the scaling factor.40 This feature is as well shared by the well-known harmonic Otto cycle.41 This is a consequence of the fact that both the cycles involve the driving of an isolated quantum system with a scale-invariant energy spectrum. However, the interaction-driven cycle generally involves a change of density of states which can be regarded as the change of the generalized exclusion statistics in an ideal gas;42 also see Supplementary Information Sections II and III. Moreover, the scale invariant character of the energy spectrum in the strongly interacting regime is also present for multicomponent Bose and Fermi gases, preserving the universal efficiency given by Eq. 6.

For weaker interactions, we resort to a numerically-exact solution of Eqs. 3 for finite particle number \(N\). We enumerate all the possible sets \({{\mathcal{I}}}_{n}\) of quantum numbers for low-energy states and solve Eqs. 3 numerically for given \(c\). With the resulting quasi-momenta \(\{{k}_{n,1},{k}_{n,2},\cdots \ \}\), and the corresponding energy eigenvalues \({\epsilon }_{n}={\sum }_{i}{k}_{n,i}^{2}\), the probability that the Bose gas at temperature \(T\) is found with energy \({\epsilon }_{n}\) is set by the Boltzmann weights \({p}_{n}={e}^{-{\epsilon }_{n}/T}/{\sum }_{m}{e}^{-{\epsilon }_{m}/T}\) (with \({k}_{{\rm{B}}}=1\)). The equilibrium energy of the states \(A\) and \(C\) is set by thermal averages of the form \(\langle E\rangle ={\sum }_{n}{p}_{n}{\epsilon }_{n}\) that in turn yield the expressions for \({Q}_{2}\) and \({Q}_{4}\). Here, a proper cutoff of the possible sets \(\{{{\mathcal{I}}}_{n}\}\) can be determined by \({p}_{n}/{p}_{{\rm{G}}}\ll 1\), where \({p}_{{\rm{G}}}\) is the probability for the working substance to be in the ground state. The numerical results for the efficiency \(\eta\) and output work \(W\) are shown in Fig. 2 as a function of the interaction strength. For fixed values of \({T}_{{\rm{C}}}\) and \({T}_{{\rm{A}}}\), the maximum work is studied as a function of \({c}_{{\rm{A}}}\) while keeping \({c}_{{\rm{B}}}\) constant (Fig. 2). The efficiency is well reproduced by Eq. 6 at strong coupling, that captures the 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.

Fig. 2 Work and efficiency as functions of the interaction strength \({c}_{{\rm{A}}}\). a Dependence of the work output \(W\) on \({c}_{{\rm{A}}}\) for different \({T}_{{\rm{A}}}\). b The efficiency \(\eta\) obtained by numerical calculation (solid lines), is compared with the analytical result (black dashed line) in Eq. 6. Here, \(N=5\), \(L=1\), \({T}_{{\rm{C}}}=150\), and \({c}_{{\rm{B}}}=200\) Full size image

Universal efficiency at low temperature

In the thermodynamic limit (where \(N\) and \(L\to \infty\) with \(n=N/L\) being kept constant), the equilibrium state of the 1D Bose gas is determined by the Yang-Yang thermodynamics23,24,42,43,44 (see “Methods”). At low energy, excitations are phonons and the behavior is described by the Tomonaga-Luttinger liquid (TLL) theory,45,46,47,48,49 in which the free energy density reads

$${\mathcal{F}}={{\mathcal{E}}}_{0}-\frac{\pi {T}^{2}}{6{v}_{{\rm{s}}}},$$ (7)

where \({{\mathcal{E}}}_{0}\) is the energy density of the ground state, \({v}_{{\rm{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=-\partial {\mathcal{F}}/\partial T=\pi T/3{v}_{{\rm{s}}}\). The expression for the heat absorbed and released are respectively given by energy differences

$${Q}_{2}=\frac{\pi L}{6{v}_{{\rm{s}}}^{{\rm{B}}}}({T}_{{\rm{C}}}^{2}-{T}_{{\rm{B}}}^{2}),$$ (8)

$${Q}_{4}=\frac{\pi L}{6{v}_{{\rm{s}}}^{{\rm{A}}}}({T}_{{\rm{D}}}^{2}-{T}_{{\rm{A}}}^{2}),$$ (9)

where \({s}_{i}=\pi {T}_{i}/3{v}_{s}^{i}\) and \({v}_{{\rm{s}}}^{i}\) with \({v}_{{\rm{s}}}^{{\rm{B}}}={v}_{{\rm{s}}}^{{\rm{C}}}\) and \({v}_{{\rm{s}}}^{{\rm{A}}}={v}_{{\rm{s}}}^{{\rm{D}}}\) denote 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

$$\xi \equiv \frac{{T}_{{\rm{A}}}}{{T}_{{\rm{B}}}}=\frac{{T}_{{\rm{D}}}}{{T}_{{\rm{C}}}}=\frac{{v}_{{\rm{s}}}^{{\rm{A}}}}{{v}_{{\rm{s}}}^{{\rm{B}}}},$$ (10)

where \({T}_{{\rm{B}}}\) and \({T}_{{\rm{D}}}\) are the temperatures at states B and D, respectively. As a result, the efficiency and work output are given by

$${\eta }_{{\rm{TLL}}}=1-\frac{{v}_{{\rm{s}}}^{{\rm{A}}}}{{v}_{{\rm{s}}}^{{\rm{B}}}}=1-\xi ,$$ (11)

$${W}_{{\rm{TLL}}}=\frac{\pi L{T}_{{\rm{C}}}^{2}}{6{v}_{s}^{B}}(1-\xi )\left(1-\frac{{\kappa }^{2}}{{\xi }^{2}}\right),$$ (12)

where \(\kappa ={T}_{{\rm{A}}}/{T}_{{\rm{C}}}\). Since \({T}_{{\rm{A}}}\;<\;{T}_{{\rm{B}}}\;<\;{T}_{{\rm{C}}}\), we have \(0\;<\;\kappa \;<\;\xi \;<\;1\). This efficiency (11) is universal regardless the microscopic physics of the working substance with a low-energy TLL description. The work output for fixed \({T}_{{\rm{A}}}\) and \({T}_{{\rm{C}}}\) is maximized at \(\xi ={\xi }_{{\rm{c}}}\) as shown in Supplementary Information Section IV,

$${\xi }_{{\rm{c}}}\simeq {(2{\kappa }^{2})}^{\frac{1}{3}}[1-{\left(\kappa /2\right)}^{\frac{2}{3}}/3]\ .$$ (13)

As the TLL theory describes the universal low-energy behavior of 1D many-body systems, Eqs. 11–13 provide a universal description of the efficiency and work of QHE with a 1D interacting working substance at low temperatures, which are applicable to any cycle in which work and heat exchange occurs in different strokes. In particular, in the strongly interacting regime, the sound velocity of 1D Bose gases is given by \({v}_{{\rm{s}}}\simeq 2\pi n(1-4n{c}^{-1}+12{n}^{2}{c}^{-2})\).24 In this regime, the result 11 thus reduces to Eq. 6. On the other hand, in the weak interaction regime, \({v}_{{\rm{s}}}\simeq 2n\sqrt{\frac{c}{n}-\frac{1}{2\pi }{\left(\frac{c}{n}\right)}^{\frac{3}{2}}}\),24 and thus the efficiency reaches the asymptotic value

$$\eta \simeq 1-\sqrt{\frac{{c}_{{\rm{A}}}}{{c}_{{\rm{B}}}}},$$ (14)

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 across the phase diagram of the 1D Bose gas and characterized the role of quantum criticality.50 The 1D Bose gas displays a rich critical behavior as a function of the temperature \(T\) and chemical potential \(\mu\) (Fig. 3). In the region of \(\mu /T\ll 0\), i.e. the mean distance between atoms is much larger than the thermal wave length \({\lambda }_{{\rm{th}}}=\sqrt{2\pi {\hslash }^{2}/(m{k}_{{\rm{B}}}T)}\), and the system behaves as a classical gas (CG). A quantum critical region (QC) emerges between two critical temperatures (see the white dashed lines in Fig. 3c) fanning out from the critical point \({\mu }_{{\rm{c}}}=0\). In this region, the energy gap closes universally as \(\Delta \sim | \mu -{\mu }_{{\rm{c}}}{| }^{z

u }\), in terms of the correlation length and dynamic critical exponents \(

u =1/2\) and \(z=2\). The TLL region characterized by a phononic spectrum is found with \(\mu >0\) and \(T\) below the right critical temperature.

Fig. 3 Performance across the phase diagram. a Efficiency \(\eta\) as a function of the density \(n\). b Average work output \(W/N\) as a function of density \(n\). c Phase diagram of a Bose gas given by a contour plot of the specific heat. The three different polygons delineate interaction-driven thermodynamic cycles in different regimes of the phase diagram corresponding to the three dashed lines in the other panels. The equilibrium states A and C are chosen with fixed \(c\) and \(T\) for the different cycles. The efficiency saturates at \(0.42\) given by Eq. 14 in the TLL regime. (d) The energy derivative \(\partial E/\partial c\) per particle as a function of the interaction strength \(c\). The areas of these three polygons correspond to the works per particle \(W/N\) for the three cycles in figure (c), respectively Full size image

For fixed cycle parameters (\({c}_{{\rm{A}}}\), \({c}_{{\rm{B}}}\), \({T}_{{\rm{A}}}\), and \({T}_{{\rm{C}}}\)), we study the performance of the engine across the QC region by changing \(n\). We numerically calculate the efficiency \(\eta\) and average work \(W/N\) by using the TBA equation 16, and show that near the QC region \(W/N\) has a maximum value; see Fig. 3. We set \({c}_{{\rm{A}}}=1\), \({c}_{{\rm{B}}}=3\), \({T}_{{\rm{A}}}=1\), and \({T}_{{\rm{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 densities \(n\simeq 0.3\), \(1.4\), and \(3.0\), respectively. In order to understand the maximum of \(W/N\), we also plot these three engine cycles (A \(\to\) B \(\to\) C \(\to\) D \(\to\) A) in the phase diagram of specific heat in the \(T-\mu\) plane; see Fig. 3 (c). The performance of the cycle is optimal at the critical region (green dashed line). Specifically, when the stroke A \(\to\) B is near the TLL boundary and the stroke C \(\to\) D is in the critical region, the average work output per particle is maximized. In this scenario, the two adiabatic processes (A \(\to\) B and C \(\to\) D) pass through zones in which the change in the energy dispersion is maximized, as shown in Fig. 3c. Qualitatively, the process A \(\to\) B passes though the region where the specific heat is maximum, whereas the C \(\to\) D evolves through the region with lowest specific heat. In the adiabatic process, the particle number \(N\) and entropy \(S\) are all fixed. Figure 3d shows the work for the different cycles, where the work for each adiabatic process is given by \(W={\int }_{{c}_{{\rm{initial}}}}^{{c}_{{\rm{final}}}}\left(\partial E/\partial c\right)dc\). Even if the process is not rigorously adiabatic,51 this maximum work output still holds, as further discussed in Supplementary Information Section IV.