### Set-up and sequence

Table of Contents

We prepare a degenerate two-component fermionic ^{6}Li gas in the two lowest-lying Zeeman substates of the electronic ground state ^{2}*S*_{1/2} confined in an elongated quasi-harmonic trap comprising a magnetic and an optical dipole trap (ODT) operating at a wavelength of 1,070 nm; for details of the experimental set-up, see ref. ^{32}. We cool the gas using evaporative cooling until reaching a temperature of 120 nK in a trap with geometric mean trap frequency \(\bar{\omega }={({\omega }_{x}{\omega }_{y}{\omega }_{z})}^{1/3}\), where typical values used for the engine performance are given in Extended Data Table 1. Evaporation takes place on the BEC side of the crossover at a magnetic field of 763.6 G. After evaporation, we hold the cloud for 150 ms to let it equilibrate. By varying the loading time of the magneto-optical trap before evaporation, we adjust the number of atoms per spin state *N*^{i} of the cloud from 1.75 × 10^{5} up to 3.0 × 10^{5}.

The quantum Pauli cycle is implemented by alternating changes of the trap frequency through the dipole trap laser and of the interaction strength through the magnetic field close to the Feshbach resonance. Extended Data Fig. 1 illustrates the experimental sequence of a single cycle. After initial preparation of the molecular BEC, we compress the trap adiabatically, increasing the laser power of the ODT during 300 ms from *P*_{A} to *P*_{B}. The trap frequencies in both radial directions increase with the square root of the laser power of the ODT whereas the trap frequency along the axial direction remains irrelevant compared with the trapping frequency of the magnetic field for all the powers used in this work. The resulting geometric mean trap frequencies are given in Extended Data Table 1. After a waiting time of 150 ms with constant trap frequency, we reach cycle point B. Then, we adiabatically increase the magnetic field strength linearly from *B*_{A} to *B*_{C} to change the interaction strength during the stroke B → C. Changing the magnetic field strength does not substantially alter the trap frequency in the axial direction. After a waiting time of 150 ms with constant magnetic field, we reach point C. In the next step, we ramp the ODT laser power to the initial value *P*_{A}, expanding the gas (C → D) in the trap. Finally, we close the cycle with a second magnetic field ramp back to *B*_{A}. After running an entire cycle, we verify that the measurement point A_{2} is equivalent to point A. In all of the measurement series, we have *B*_{A} < *B*_{C}, *P*_{A} < *P*_{B} and, therefore, \({\bar{\omega }}_{{\rm{A}}} < {\bar{\omega }}_{{\rm{B}}}\). Extended Data Table 1 summarizes the experimental parameters for the three types of cycles that we run.

### Atom-number measurement

To determine the number of atoms from in situ absorption pictures, we image the energetically higher-lying spin state *m*_{I} = 1 by standard absorption imaging on a charge-coupled device camera. This procedure is identical for all the regimes considered. Since the number of atoms of the spin state *m*_{I} = 0 is the same as for *m*_{I} = 1, the total number of atoms *N* (sum of atoms with spin up and down) is twice the number of atoms of the measured spin state *N*^{i}. Importantly, in the case of a molecular BEC, the number of molecules is half the number of total atoms. We extract the column density from the measured absorption pictures. Adding up the density pixelwise and including the camera’s pixel area and the imaging system’s magnification, we obtain the number of atoms. We include additional laser-intensity-dependent corrections^{51}. The imaging system has a resolution of 2.2 μm. The optical absorption cross-section changes its value throughout the BEC–BCS crossover, because it depends on the optical transition frequency^{26} which is, in turn, a function of the magnetic field strength. We determine the absorption cross-section for a magnetic field of 763.6 G and we use this value for further magnetic field strengths. This leads to the fact that the measured number of atoms on resonance is higher than the actual number of atoms. Therefore, we determine a correction factor to compensate this imperfection in independent measurements at different magnetic field strengths (Extended Data Fig. 2). We prepare an atomic cloud of varying atom number at *B*_{A} = 763.6 G. We then measure the atom number for this field at points A and B. We repeat the measurement for *B*_{C} and, to exclude atom loss, we ramp the field from *B*_{B} = *B*_{A} adiabatically to *B*_{C} and back to *B*_{B} before measuring. The atom numbers for *B*_{B} with and without the additional ramp to C are almost identical, from which we conclude that the number measured at *B*_{C} has to be the same. The correction factor is then calculated from the difference in atom number between these two data sets.

Extended Data Fig. 3 shows the measured number of atoms for one spin state *N*^{i} for the different points of the cycle of the Pauli engine as a function of the number of atoms of one spin state \({N}_{{\rm{A}}}^{i}\) in point A. The correction factors for different magnetic fields are included. In this way, we independently verify that we experimentally have a constant number of atoms during the cycle from A → D through B and C. Only the last stroke D → A_{2} suffers from atom losses of about 10%. These losses primarily occur when the quantum statistics are changed from Fermi–Dirac to Bose–Einstein in a shallow trap. The thermodynamics of adiabatic transfers through the crossover has been theoretically studied in ref. ^{39}. During the change of statistics, bosonic molecules are formed. We attribute the losses to an excess-energy transfer between molecules during the ramping. Owing to the relatively low trap depth (needed for a high compression ratio), even small kinetic energies are sufficient to remove molecules from the trap.

The measured atom number has an uncertainty due to several mechanisms. First, the produced quantum gases feature a fluctuating number of particles due to uncontrolled technical fluctuations in, for example, cooling- and trapping-laser intensity or magnetic-field currents. But it is also due to physical statistical fluctuations, such as Poissonian fluctuations of the particle number during laser cooling or fluctuations occurring during the phase transition to a quantum fluid. These mechanisms cause fluctuations of the particle number of the produced ultracold clouds. Second, even for identical atom numbers, an additional uncertainty originates from the measurement process. We deploy resonant high-intensity absorption imaging^{52,53} on the transition \({| }^{2}{{\rm{S}}}_{1/2},{m}_{J}=-\,1/2,{m}_{I}=1\rangle \leftrightarrow {| }^{2}{{\rm{P}}}_{3/2},{m}_{J}=-\,3/2,{m}_{I}=1\rangle \) to acquire the column-integrated density distribution. Determining the accurate atom number from the resulting images requires precise control over the imaging pulse length and power, as well as calibration of the camera counts and effective absorption cross-section. In particular, fluctuations of the laser power after transmission through the imaging system in addition to camera noise are the main origin of statistical density uncertainties and dominate the uncertainty indicated in the measured data points. A systematic uncertainty is due to limited knowledge of the absorption cross-section. It can, in principle, be determined by the transition wavelength, but its value can change owing to optical pumping, the presence of magnetic fields or polarization effects, and was, therefore, determined by comparing theoretical density distributions of a BEC with measured ones^{51}. The correction factor for our specific system was evaluated from a *χ*^{2} fit and has a systematic uncertainty of less than 10%. The resulting statistical atom-number fluctuations dominate the uncertainty for identical experimental parameters. We quantify them by the 1*σ* standard deviation of typically 20 identical realizations and indicate them as bars in equations (2) and (3) as well as in Figs. 2 and 4 around the mean atom number of the respective measurement series.

### Cloud-radii determination

We obtain line densities of the quantum gases by integrating the column-integrated density distributions of the two-dimensional (2D) absorption pictures *n*_{2D}(*x*, *y*) = ∫ *n*_{3D}(*x*, *y*, *z*) d*z* along the *x* and *y* directions separately. We fit these distributions with a 1D-integrated Thomas–Fermi profile. The Thomas–Fermi profiles are different for the interaction regimes throughout the BEC–BCS crossover. In the Thomas–Fermi limit, the in situ density distribution of a molecular BEC has the shape^{9}

$${n}_{{\rm{mBEC}}}\propto {\left(1-\frac{{x}^{2}}{{R}_{{{\rm{mBEC}}}_{x}}^{2}}\right)}^{2},$$

(2)

where \({R}_{{{\rm{mBEC}}}_{x}}\) is the Thomas–Fermi radius of the atomic cloud in the *x* direction. For a resonantly interacting Fermi gas right at \({B}_{{\rm{res}}}\), the density profile \({n}_{{\rm{res}}}\) in the *x* direction can be written as^{34}

$${n}_{{\rm{res}}}\propto {\left(1-\frac{{x}^{2}}{{R}_{{{\rm{res}}}_{x}}^{2}}\right)}^{5/2},$$

(3)

with the rescaled Thomas–Fermi radius \({R}_{{{\rm{res}}}_{x}}\) in the *x* direction. Therefore, we determine the measured cloud radii in both regimes by fitting the measured density profiles with the appropriate line shapes. The extracted radii in the *x* and *y* directions are shown in Extended Data Fig. 4. Our experimental set-up does not allow measurements of the radii *R*_{z} in the *z* direction. For this radial direction, however, we use the measured radius *R*_{x} in the *x* direction (radial) and correct the value with the corresponding ratio of the trap frequencies *R*_{z} = (*ω*_{x}/*ω*_{z})*R*_{x} (refs. ^{46,54}).

### Temperature measurements

We determine the temperature of the molecular BEC by means of a bimodal fit of the density profile (for more details, see ref. ^{32}) at a magnetic field strength of 680 G. This value of *B* is chosen by considering the following trade-off: at lower magnetic fields, losses in the atom clouds are too high for quantitative temperature determination, whereas for higher magnetic fields, the condensate and thermal parts are not well separated, which complicates a bimodal fit.

Owing to three-body recombination losses in the range between 550 G and 750 G (ref. ^{55}), molecule losses in the cloud are already significant at 680 G. To avoid these losses, we choose a field of 763.6 G as the start of the Pauli cycle. To determine the temperature in an independent measurement, we ramp the field to 680 G during 200 ms and determine the temperature there. However, since decreasing the magnetic field value adiabatically throughout the BEC–BCS crossover increases the temperature *T* of the gas^{39,48}, the measured temperature at 680 G is an upper bound for the temperature at point A and beyond. For the work stroke between points A and B, we increase the mean trap frequency. We observe that the reduced temperature *T*/*T*_{F} of the gas does not show significant changes during these work strokes. The temperature *T* does not change after ramping the trapping frequency back and forth. Because of this, we determine the temperature of the cloud, as mentioned above, with a bimodal fit for the two settings. First, we ramp the magnetic field from 763.6 G (value at point A) to 680 G to determine the temperature there. Second, we restart the sequence and run the cycle until point B. Afterwards, we change the direction of the cycle and go directly back to point A. A magnetic field sweep to 680 G allows again a temperature measurement. Extended Data Fig. 5 shows that the temperatures of the cloud after reversal to point A (cyan points) lie within the error (grey shaded area) of the initial temperature at point A (black line). We observe that an increase in temperature of less than 10% takes place for higher ratios of the mean trapping frequency, which can be neglected in the analysis (see below).

The method described above holds for atomic clouds below the critical temperature. For a thermal gas, we can directly determine the temperature of the magnetic field at the cycle point. At finite temperatures, density profiles of thermal clouds can be approximated with a classical Boltzmann distribution. The temperature *T*_{j} can be calculated dependent on direction *j* of the cloud \({T}_{j}=(m{\bar{\omega }}^{2}{{\sigma }}_{j}^{2})/{k}_{{\rm{B}}}\) (ref. ^{46}), with *σ*_{j} as width of a Gaussian determined by fitting the 1D density profiles with \({n}_{{\rm{thermal}}}\propto \exp (-{j}^{2}/2{{\sigma }}_{j}^{2})\). In the case of molecular bosons, we use 2*m* for the mass instead of *m*. Hence, interactions directly influence the width of the cloud; we interpret this temperature for our experiment as an approximated temperature.

### Energy calculation

The calculation of the total energy *E*_{mBEC} of a molecular BEC is based on the Gross–Pitaevskii equation in the Thomas–Fermi limit for a harmonic trapping potential and for zero temperature. This energy consists of two parts, the kinetic energy *E*_{kin} and the energy in the Thomas–Fermi limit *E*_{TF}, which takes into account the oscillator and interaction energies^{46,54}

$$\begin{array}{l}{E}_{{\rm{mBEC}}}\,=\,{E}_{{\rm{TF}}}+{E}_{{\rm{kin}}}\\ \,\,=\,\left(\frac{{\mu }_{{\rm{mBEC}}}}{7}+\frac{{\hbar }^{2}}{4m{R}_{{\rm{mBEC}}}^{2}}\,\left[{\rm{ln}}\frac{{R}_{{\rm{mBEC}}}}{1.3{a}_{{\rm{ho}}}}+\frac{1}{4}\right]\right)\frac{5N}{2},\end{array}$$

(4)

where *R*_{mBEC} is the geometric mean radius of the cloud. The chemical potential is given by

$${\mu }_{{\rm{mBEC}}}=\frac{\hbar \bar{\omega }}{2}{\left(\frac{15\frac{N}{2}{a}_{{\rm{dd}}}}{{a}_{{\rm{ho}}}}\right)}^{2/5},$$

(5)

where *a*_{dd} = 0.6*a* is the *s*-wave scattering length for molecules^{8} and \({a}_{{\rm{ho}}}=\sqrt{\hbar /(2m\bar{\omega })}\) is the oscillator length. The radii and molecule numbers are extracted from the absorption pictures for the considered interaction strengths. When including the molecular energy, which gives the amount of energy needed to dissociate a molecule into two atoms, the total energy of the molecular BEC *E*_{mBEC,m} is given by

$${E}_{{\rm{mBEC}},{\rm{m}}}={E}_{{\rm{mBEC}}}-\frac{N}{2}\frac{{\hbar }^{2}}{m{a}^{2}}.$$

(6)

On resonance, the scattering length diverges and the binding energy vanishes. When we tune the magnetic field to resonance, the system evolves into a strongly interacting unitary Fermi gas. Its total energy \({E}_{{\rm{res}}}\) is related to that of an ideal binary Fermi gas through scaling with the universal constant \(\sqrt{1+\beta }\)

$${E}_{{\rm{res}}}=\frac{3}{4}\sqrt{1+\beta }{E}_{{\rm{F}}}N,$$

(7)

where \({E}_{{\rm{F}}}={(3N)}^{1/3}\hbar \bar{\omega }\) is the Fermi energy^{10,33,34,56}. This zero-temperature expression provides a good approximation for *T*/*T*_{F} < 0.2 as supported by the energy measurements on resonance reported in refs. ^{56,57,58} (the low-*T* measurements are almost independent of *T* in the mentioned range). It should be noted that owing to the different temperature dependency of the entropy of the bosonic (*S*^{bosons} ∝ *T*^{3}) and fermionic (*S*^{fermions} ∝ *T*) trapped gases, a decrease in the temperature is expected when performing an adiabatic sweep from the molecular BEC side to the BCS side and to unitarity. This internal temperature decrease is another indication that the operation of the Pauli engine is not driven by thermal energy. Since the condition *T*/*T*_{F} < 0.2 is fully satisfied for points C and D in all of our experiments the zero-*T* formulas on resonance are highly accurate.

The obtained experimental values were contrasted with those obtained using well-known formulas for BECs^{46,54} in the Thomas–Fermi limit for interacting bosons adapted as previously explained for composite bosons made up of fermions of opposite spin states^{5,6,59}. In these cases, we use the theoretically expected radius and a constant number of particles during the cycle equal to the experimental number of particles at point A. Using equations (4), (5) and \({R}_{{\rm{mBEC}}}={(9{\hbar }^{2}Na/(8{m}^{2}{\bar{\omega }}^{2}))}^{1/5}\) in equation (6), we obtain the following energy for the molecular BEC regime at zero temperature.

$${E}_{{\rm{mBEC}},{\rm{m}}}\,=\,\frac{N}{2}\left\{\frac{5}{14}\hbar \bar{\omega }{\left(\frac{9\frac{N}{2}a}{\sqrt{\frac{\hbar }{2m\bar{\omega }}}}\right)}^{\frac{2}{5}}+\frac{5}{4}\frac{{\hbar }^{2}}{m{\left({\left({\left(\frac{{\hbar }^{2}}{2m\bar{\omega }}\right)}^{2}9\frac{N}{2}a\right)}^{\frac{1}{5}}\right)}^{2}}\left[{\rm{ln}}\left(\frac{10}{13}{\left(\frac{9\frac{N}{2}a}{\sqrt{\frac{\hbar }{2m\bar{\omega }}}}\right)}^{\frac{1}{5}}\right)+\frac{1}{4}\right]\right\}-\frac{N}{2}\frac{{\hbar }^{2}}{m{a}^{2}},$$

(8)

where the first two terms are related to the energy of bosonic particles in a trap including the interaction between bosons (the molecules interact with each other by means of a contact interaction^{46} that can be quantified using the interaction strength^{8} *g* = 1.2π*ħ*^{2}*a*/*m*) and the last term is the contribution of the molecular energy of each of the pairs. For the zero-*T* calculations of the ground-state energy of the bosonic system, we also numerically solve the Gross–Pitaevskii equation with the bosonic interaction in terms of the dimer–dimer scattering length. Owing to the range of the experimental parameters, the Thomas–Fermi approximation holds in all our experiments in the molecular BEC regime; therefore, the numerical results are the same as those obtained when using equation (8) (refs. ^{46,54}). Extended Data Fig. 6 shows the experimental and theoretical energies for each point of the Pauli cycle for the data set of Fig. 3 (point A denotes the first point of the cycle and point A_{2} the last point after a full cycle).

For low but non-zero temperature, the Thomas–Fermi approximation to the Gross–Pitaevskii equation reads^{46}

$$\begin{array}{l}{E}_{{\rm{mBEC}},{\rm{m}}}^{T}\,=\,\frac{N}{2}\hbar \bar{\omega }{\left(\frac{\frac{N}{2}}{\zeta (3)}\right)}^{\frac{1}{3}}\left\{\frac{3\zeta (4)}{\zeta (3)}{\left(\frac{T}{\frac{\hbar \bar{\omega }}{{k}_{{\rm{B}}}}{\left(\frac{\frac{N}{2}}{\zeta (3)}\right)}^{\frac{1}{3}}}\right)}^{4}+\frac{\zeta {(3)}^{\frac{1}{3}}}{14}{\left(9{\left(\frac{N}{2}\right)}^{\frac{1}{6}}\frac{a}{\sqrt{\frac{\hbar }{2m\bar{\omega }}}}\right)}^{\frac{2}{5}}\right.\\ \,\,{\left(1-{\left(\frac{T}{\frac{\hbar \bar{\omega }}{{k}_{{\rm{B}}}}{\left(\frac{\frac{N}{2}}{\zeta (3)}\right)}^{\frac{1}{3}}}\right)}^{3}\right)}^{\frac{2}{5}}\left(5+16{\left(\frac{T}{\frac{\hbar \bar{\omega }}{{k}_{{\rm{B}}}}{\left(\frac{\frac{N}{2}}{\zeta (3)}\right)}^{\frac{1}{3}}}\right)}^{3}\right)\}-\frac{N}{2}\frac{{\hbar }^{2}}{m{a}^{2}}.\end{array}$$

(9)

According to Tan’s generalized virial theorem^{10,42,60,61,62}, the energy for a cigar-shaped trapped Fermi gas is given by \(E=Nm(2{\omega }_{r}^{2}{\langle r}^{2}\rangle \,+\)\({\omega }_{z}^{2}{\langle z}^{2}\rangle )-{\hbar }^{2}{\mathcal{I}}/(8{\rm{\pi }}ma)\), where *r* and *z* stand for the radial and axial coordinates, respectively, and \({\mathcal{I}}\) is the contact (a measure of the probability for two fermions with opposite spins being close together^{10}). This expression is valid for any value of the scattering length *a* and for any temperature *T*. On resonance, that is, at unitarity, we have 1/*k*_{F}*a* = 0 and a finite \({\mathcal{I}}\). Therefore the virial theorem reduces to the case presented by Thomas et al. in refs. ^{40,41}, that is, \(E=Nm(2{\omega }_{r}^{2}{\langle r}^{2}\rangle +{\omega }_{z}^{2}{\langle z}^{2}\rangle )\). In the weak coupling limit for the molecular BEC side, the contact reduces to \({\mathcal{I}}\approx 4\pi N/a\), in which case the last term of the virial expression gives the total molecular energy and the virial energy reduces to equations (8) and (9) (refs. ^{42,46}). All of this means that the expression −*ħ*^{2}/(*m**a*^{2}) for the molecular term holds even for *T*/*T*_{F} ≈ 0.7 as long as the experiments are realized in the deep molecular BEC regime. For our magnetic fields, this condition is fulfilled for the Pauli engine and Feshbach cycle as well as for the thermal cycle. Furthermore, measurements of the molecular energy were presented in ref. ^{4} for *T*/*T*_{F} ≈ 0.15, which show that the formula −*ħ*^{2}/(*m**a*^{2}) holds for *T*/*T*_{F} ≈ 0.2 in a complete *B* sweep from molecular BEC to unitarity.

In the high-*T* regime the interactions are negligible and the distributions can be approximated by Maxwell–Boltzmann distributions. When the thermal energy *k*_{B}*T* is large enough to break all the pairs, the energy can be approximated by that of an ideal classical gas of atomic particles for both magnetic fields. In this case, the temperature in the work strokes does not change and there is no work output. When the thermal energy is below the binding energy, the energy for a magnetic field below resonance corresponds to a classical gas of molecules with mass 2*m* whereas the energy at the resonant field is given by a classical gas of atoms in a harmonic trap. The molecular term cancels in the total work output giving \(W=3{k}_{{\rm{B}}}N\left({T}_{{\rm{A}}}/2-{T}_{{\rm{D}}}+{T}_{{\rm{C}}}-{T}_{{\rm{B}}}/2\right)\). The temperature of the gas is obtained through a Gaussian fitting of the density profile leading to \({k}_{{\rm{B}}}T=m{\bar{\omega }}^{2}{\sigma }^{2}\), where *σ* is the width of the fitted density profile. Since \({\bar{\omega }}_{{\rm{D}}}={\bar{\omega }}_{{\rm{A}}}\) and \({\bar{\omega }}_{{\rm{C}}}={\bar{\omega }}_{{\rm{B}}}\) and the mass of the molecules is twice that of one of the atoms, we find \({T}_{{\rm{A}}}/2={({\sigma }_{{\rm{A}}}/{\sigma }_{{\rm{D}}})}^{2}{T}_{{\rm{D}}}\) and \({T}_{{\rm{B}}}/2={({\sigma }_{{\rm{B}}}/{\sigma }_{{\rm{C}}})}^{2}{T}_{{\rm{C}}}\), that is, \(W=3{k}_{{\rm{B}}}N\left\{\left({\sigma }_{{\rm{A}}}^{2}\,/{\sigma }_{{\rm{D}}}^{2}-1\right){T}_{{\rm{D}}}+\left(1-{\sigma }_{{\rm{B}}}^{2}\,/{\sigma }_{{\rm{C}}}^{2}\,)\right){T}_{{\rm{C}}}\right\}\). The work vanishes for *σ*_{A} = *σ*_{D} and *σ*_{B} = *σ*_{C}. In our experiments, the widths do not show significant statistical differences and therefore the work output vanishes for the engine running in the high-*T* regime.

Let us now detail each one of the theoretical curves presented in the main text. Based on Tan’s generalized virial theorem, we calculate the theoretical trap energy of Fig. 2 by means of equation (9) without the molecular term and using the experimental temperatures at points A and B, whereas for the points C and D we use the energy for zero *T* at unitarity given by equation (7). The vanishing efficiency for the thermal case is expected from the classical gas argument of the previous paragraph. The cyan curves in Fig. 3 arise from numerical calculations of the Gross–Pitaevskii equation. The black dashed curves in Fig. 3e,f were also calculated by solving numerically the Gross–Pitaevskii equation and yield the same results as an ideal molecular gas of point-like bosons having an infinite binding energy, which leads to \({W}_{{\rm{non}}-{\rm{interacting}}}={(3N)}^{4/3}\sqrt{1+\beta }\hbar ({\bar{\omega }}_{{\rm{B}}}-{\bar{\omega }}_{{\rm{A}}})/4\) and *η*_{non-interacting} = 0. The purple curves in Fig. 3e,f correspond to the results obtained for ideal Bose and Fermi gases in the highly degenerate regime and lead to

$${W}_{{\rm{ideal}}}={(3N)}^{4/3}\hbar ({\bar{\omega }}_{{\rm{B}}}-{\bar{\omega }}_{{\rm{A}}})/4$$

(10)

and

$${\eta }_{{\rm{ideal}}}=1-{\bar{\omega }}_{{\rm{B}}}/{\bar{\omega }}_{{\rm{A}}}.$$

(11)

The ideal efficiency is similar to the maximum Otto efficiency^{25}. Compared with these upper limits, we find that the experimental system shows reduced values but of the same order of magnitude. Our model allows us to predict the performance for different magnetic field values in the molecular BEC regime. The grey curves in Fig. 3e,f and the theoretical curves of Fig. 4 follow from equations (8) and (7). Since the latter calculations rely on zero*-T* formulas, we interpret their results as an upper bound or limit for the work output and efficiency. To show that this estimation is accurate, we need to consider two aspects. First, for the two on-resonance points (C, D), we consider the zero-*T* calculations to be a good approximation because of the aforementioned decay in the temperature between the molecular BEC regime and unitarity in the adiabatic sweep. Owing to the dominance of the molecular term on the molecular BEC side, the efficiency obtained when including the temperature overlaps with that obtained within the zero-*T* approximation. For the data of Fig. 4, the zero-*T* calculations give an efficiency of about 0.07(2) while the finite-*T* formulas give an efficiency of about 0.05(3). A null hypothesis test for the mean difference gives a *P* value of 0.1. Setting the usual significance of *α* = 0.05, the null hypothesis of equal efficiencies cannot be rejected. Therefore, we conclude that the finite-*T* corrections have no notable effects on the efficiency of our engine.

Figure 3e,f shows that the work output is reduced with increased initial effective repulsion of the molecular BEC, whereas the efficiency increases. This scaling of the work output is a consequence of the competition between interactions among molecules in the initial molecular BEC state and the effect of changing the quantum statistics. For stronger initial effective repulsion, the molecular BEC cloud already exhibits a relatively large energy in the trap, so that the change of quantum statistics can only contribute a comparatively smaller amount of energy during the Pauli stroke. This suggests an optimal work output for an initially non-interacting molecular BEC. It is important to mention that even though the binding energy has no consequences for the work output (the *s*-wave scattering length does not change during the work strokes) it still has to be provided to the system during the Pauli stroke when dissociating a molecule into two atoms; hence, it must be included in the energy-cost calculation of the engine’s efficiency. This binding energy quickly grows as the magnetic field deviates from the resonance value and the associated energy cost quickly reduces the efficiency. For the experimentally inaccessible case of a non-interacting molecular BEC, this binding energy is so large that the efficiency of the Pauli cycle is essentially zero. These considerations point toward an optimal point of operation which might, additionally, be temperature dependent.

An important result that demonstrates the non-classical drive of the Pauli engine is reproduced by both experimental findings and the theoretical model: the work output as a function of the number of particles *W* ∝ *N*^{α} scales with a fitted exponent *α* = 1.73(18), which is close to the prediction of 1.4 given by theory (Fig. 4a). The exponent of the work output with particle number is different from the 1D case mentioned in the introduction due to a modified density of states, that is, the number of available single-particle states in the potential. Although a similar exponent is expected for the Feshbach-driven stroke which remains always in the molecular BEC regime, the efficiency for this engine is close to zero. Most importantly, the exponent is larger than one, which is the expected exponent for a non-interacting, classical gas^{1,48,63}.

### Pressure and volume calculations

The pressure of the system is calculated as *p* = −∂*E*/∂*V*. As explained in the previous paragraph, the zero-*T* energy constitutes a good approximation for the system and, therefore, we calculate the pressure for *T* = 0. For computing the derivative at the experimentally obtained volumes, we use that, on the molecular BEC side, the theoretical volume is given by

$${V}_{{\rm{m}}{\rm{B}}{\rm{E}}{\rm{C}}}=\frac{4}{3}\pi {\left(\frac{9{\hbar }^{2}Na}{8{m}^{2}{\bar{\omega }}^{2}}\right)}^{\frac{3}{5}},$$

(12)

whereas, on resonance, the volume has a different dependence on the number of atoms and on the geometric mean of the trap frequency, and is given by

$${V}_{{\rm{r}}{\rm{e}}{\rm{s}}}=\frac{4}{3}\pi {(1+\beta )}^{\frac{3}{4}}\sqrt{24\frac{{\hbar }^{3}N}{{m}^{3}{\bar{\omega }}^{3}}},$$

(13)

(refs. ^{34,46,54,64}). The experimental volume in both cases is computed as the volume of an ellipsoid having the radii obtained by means of equations (2) and (3). Extended Data Fig. 7 shows the corresponding values of volume *V* and pressure *p* as a function of the compression ratio. The *p–**V* diagram in Fig. 2c is obtained using these experimental data (points) and by calculating the corresponding theoretical curves. For completeness, we also checked that the work output obtained from the area inside the *p–**V* diagram −∫ *p* d*V* is in agreement with that shown in Fig. 3e.

A closer analogy with a quantum Otto cycle can be drawn using generalized parameters for *p* and *V* (refs. ^{65,66}). It will be interesting in the future to extend this concept to ensembles with changing quantum statistics and to evaluate the performance of the Pauli engine according to this description.