Quantum Statistical Enhancement of the Collective Performance of Multiple Bosonic Engines

TL;DR

This paper demonstrates that quantum statistical effects in ensembles of indistinguishable bosonic engines can enhance collective thermodynamic performance, unlike distinguishable systems, due to permutation symmetry.

Contribution

It introduces a quantum statistical mechanism for collective enhancement in bosonic quantum engines, highlighting the role of permutation symmetry.

Findings

Quantum statistical effects lead to performance enhancement in bosonic engines.

Permutation symmetry is key to the collective thermodynamic advantage.

Enhancement is absent in systems with distinguishable resources.

Abstract

We consider an ensemble of indistinguishable quantum machines and show that quantum statistical effects can give rise to a genuine quantum enhancement of the collective thermodynamic performance. When multiple indistinguishable bosonic work resources are coupled to an external system, the internal energy change of the external system exhibits an enhancement arising from permutation symmetry in the ensemble, which is absent when the latter consists of distinguishable work resources.

Equations

$$H(t)=H_R(t)+H_C(t)+H_S,$$

$$H_C(t)=g_C(t)V_R\otimes V_S,$$

$$p_i=\mathrm{Tr}{}_R\left[{}_S⟨i|U^{(I)}(T,0){\rho }_0{U^{(I)}}^{†}(T,0)|i{⟩}_S\right].$$

$$p_i\simeq {\int }_0^{T}dt{\int }_0^{T}d{t}^{\prime }$$

$$⟨V_R^{(I)}(t^{\prime })V_R^{(I)}(t){⟩}_{{\rho }_0^R},$$

$$H_C(t)=g_C(t)2S_x\otimes V_S,$$

$$H_E(t)=2\Omega (t)S_z+2\Delta S_x,$$

$$U_R(t,t_0)\approx \sum _{m=-N/2}^{N/2}|m,{\theta }_t{⟩}_E{}_E⟨m,{\theta }_{t_0}|e^{-im\varphi (t,t_0)},$$

$$⟨V_R^{(I)}(t^{\prime })V_R^{(I)}(t){⟩}_{{\rho }_{t_0}^R}=4\cos {\theta }_t\cos {\theta }_{t^{\prime }}⟨m^2{⟩}_{t_0}$$

$$+\sin {\theta }_t\sin {\theta }_{t^{\prime }}\sum _{\sigma =\pm }{e}^{-i\sigma \varphi (t^{\prime },t)}\left[\frac{N}{2}\left(\frac{N}{2}+1\right)-F_{\sigma }(N,{\beta }_{t_0}{E}_{t_0})\right],$$

$$⟨V_R^{(I)}(t){⟩}_{{\rho }_{t_0}^R}=2\cos {\theta }_t⟨m{⟩}_{t_0}.$$

$$⟨V_R^{(I)}(t^{\prime })V_R^{(I)}(t){⟩}_{{\rho }_{t_0}^R}=$$

$$+\frac{N}{2}\frac{\sin {\theta }_t\sin {\theta }_{t^{\prime }}}{\cosh ({\beta }_{t_0}{E}_{t_0})}\sum _{\sigma =\pm }{e}^{\sigma \left[i\varphi (t^{\prime },t)-{\beta }_{t_0}{E}_{t_0}\right]},$$

$$⟨V_R^{(I)}(t){⟩}_{{\rho }_{t_0}^R}=-N\cos {\theta }_t\tanh ({\beta }_{t_0}{E}_{t_0}),$$

$$⟨\Delta U{⟩}_N\simeq g^2{⟨[V_R^{(I)}(t_1){]}^2⟩}_{{\rho }_0^R}\sum _{i\ne 0}{ϵ}_i^S{\left|{}_S⟨i|V_S^{(I)}(t_1)|0{⟩}_S\right|}^2.$$

$$p_i^{\mathrm{indist}}$$

$$p_i^{\mathrm{dist}}$$

$$|m,{\theta }_t{⟩}_E=e^{-i({\theta }_t+\frac{\pi }{2})S_y}|m{⟩}_z,$$

$$V_R^{(I)}(t)$$

$$=2\sum _{m,m^{\prime }}{e}^{i(m^{\prime }-m)\varphi (t,t_0)}{}_z⟨m^{\prime }|e^{i({\theta }_t+\frac{\pi }{2})S_y}{S}_x{e}^{-i({\theta }_t+\frac{\pi }{2})S_y}|m{⟩}_z|m^{\prime },{\theta }_{t_0}{⟩}_E{}_E⟨m,{\theta }_{t_0}|$$

$$=2\sum _{m,m^{\prime }}{e}^{i(m^{\prime }-m)\varphi (t,t_0)}{}_z⟨m^{\prime }|(-\sin {\theta }_t{S}_x+\cos {\theta }_t{S}_z)|m{⟩}_z|m^{\prime },{\theta }_{t_0}{⟩}_E{}_E⟨m,{\theta }_{t_0}|$$

$$=2\cos {\theta }_t{S}_z(t_0)-\sin {\theta }_t\left[S_{+}(t_0)e^{i\varphi (t,t_0)}+S_{-}(t_0)e^{-i\varphi (t,t_0)}\right],$$

$$\varphi (t,t_0)={\int }_{t_0}^{t}d{t}^{\prime }2E_{t^{\prime }}.$$

$$⟨V_R^{(I)}(t){⟩}_{{\rho }_{t_0}^R}$$

$$=2\cos {\theta }_t\frac{{\sum }_{m}m{e}^{-2m{\beta }_{t_0}{E}_{t_0}}}{Z_{{\beta }_{t_0}}}=2\cos {\theta }_t⟨m{⟩}_{t_0},$$

$$V_R^I(t^{\prime })V_R^I(t)=$$

$$+\sin {\theta }_t\sin {\theta }_{t^{\prime }}\left[S_{+}^2(t_0)e^{i\varphi (t,t_0)+i\varphi (t^{\prime },t_0)}+h.c.\right]-2\cos {\theta }_{t^{\prime }}\sin {\theta }_t{S}_z(t_0)\left[S_{+}(t_0)e^{i\varphi (t,t_0)}+h.c.\right]$$

$$-2\cos {\theta }_t\sin {\theta }_{t^{\prime }}\left[S_{+}(t_0)e^{i\varphi (t,t_0)}+h.c.\right]S_z(t_0),$$

$$⟨m^2{⟩}_{t_0}\equiv \frac{{\sum }_m{m}^2{e}^{-2m{\beta }_{t_0}{E}_{t_0}}}{Z_{{\beta }_{t_0}}},$$

$$H_E(t)=\sum _{j=1}^N{H}_{E,j}(t)=\sum _{j=1}^N\left[\Omega (t){\sigma }_{j,z}+\Delta {\sigma }_{j,x}\right],$$

$$U_R(t,t_0)\approx \underset{j=1}{\overset{N}{⨂}}\sum _{m=-1/2}^{1/2}|m,{\theta }_t{⟩}_j{}_j⟨m,{\theta }_{t_0}|e^{-im\varphi (t,t_0)}.$$

$$|m,{\theta }_t{⟩}_j=e^{-i({\theta }_t+\pi /2){\sigma }_{j,y}/2}|m{⟩}_{jz},$$

$$V_R^{(I)}(t)=\sum _{j=1}^N\left\{\cos {\theta }_t{\sigma }_{j,z}(t_0)-\sin {\theta }_t\left[{\sigma }_{j,+}(t_0)e^{i\varphi (t,t_0)}+{\sigma }_{j,-}(t_0)e^{-i\varphi (t,t_0)}\right]\right\},$$

$${\rho }_{t_0}^E=\frac{{\prod }_{j=1}^N{e}^{-{\beta }_{t_0}{E}_{t_0}{\sigma }_{j,z}(t_0)}}{{\cosh }^N({\beta }_{t_0}{E}_{t_0})},$$

$$⟨V_R^{(I)}(t){⟩}_{{\rho }_{t_0}^R}=\frac{{\sum }_{j=1}^N\mathrm{Tr}[\cos {\theta }_t{\sigma }_{j,z}(t_0)e^{-{\beta }_{t_0}{E}_{t_0}{\sigma }_{j,z}(t_0)}]}{{\cosh }^N({\beta }_{t_0}{E}_{t_0})}.$$

Full text

Introduction.— Technological advances have allowed us to miniaturize thermal machines to the nanoscale and beyond, where quantum effects can play an important role QTroadmap18 . A paradigmatic instance of a thermal machine is a quantum heat engine (QHE). First conceived in the late 1950s, a QHE is a quantum system that serves as the working fuel of a thermodynamic cycle scovil59 ; alicki79 ; kosloff84 ; bender00 ; quan07 ; kosloff14 . More recently, a synergy between technology and progress on the foundations of quantum thermodynamics Talkner07 ; Campisi11 ; Goold16 ; vinjanampathy16 ; strasberg17 ; alicki18 has led to a surge of activity on the study of quantum machines Scully03 ; Scully11 ; Abah12 ; Abah14 ; Rossnagel14 ; Rossnagel16 ; Brantut13 ; maslennikov19 ; Klatzow19 , consolidating it as an active area of research ghosh18 . A prominent challenge in this context is the identification of situations where quantum effects govern and lead to an enhanced performance with no classical counterpart Scully11 ; jaramillo16 ; Yi17 ; Ding18 . One strategy to identify such situations is to consider thermal machines composed of multiple components Sisman01 ; Saygin01 ; gong14 ; zheng15 ; binder15 ; hardal15 ; Beau16 ; jaramillo16 ; uzdin16 ; campaioli17 ; hardal18 ; niedenzu18 ; andolina19 ; gelbwaser19 ; Klatzow19 ; manatuly18 ; chen18 ; myers20 described by collective quantum states with mutual coherence. In this setting, the enhancement requires either interaction among the components or the performance of collective unitary operations on the constituents. Furthermore, the natural process of extracting work from a QHE by outcoupling it to drive another quantum system levy16 ; watanabe17 or even the process of storing work in a quantum system binder15 ; andolina19 can also lead to the manifestation of genuine quantum effects. In this Letter, we identify a third route, in which quantum statistics leads to a genuine quantum enhancement of the performance as a result of the statistical indistinguishability of the constituent work resources. Specifically, we consider multiple work resources, each composed of a single QHE with an individual piston Abah12 ; gelbwaser14 , coupled to a single external system and show that the internal energy change of the external system displays quantum enhancement when the QHEs are indistinguishable. We note that such a setting is fundamentally different from a single QHE with a working fluid consisting of multiple particles Kim11 ; zheng15 ; Beau16 ; jaramillo16 ; Bengtsson18 ; hardal18 ; Deng18 ; chen18 .

Setup.— Consider a collective work resource $R$ made of $N$ heat engines $E_{1}$, $\cdots$, $E_{N}$ interacting with two heat baths ($B_{1}$ and $B_{2}$) and an external quantum system $S$ on which the work is performed. The coupling between $R$ and $S$ is solely established via the heat engines; see Fig. 1. The global Hamiltonian of the whole system is the sum of that of the work resources, the external system, and the coupling $C$ between them:

[TABLE]

where the external system is assumed to be time independent. If the work resources are QHEs, $H_{R}(t)$ collectively represents the Hamiltonian for $N$ engines and the two common baths. For simplicity, we consider the following form of the coupling Hamiltonian $H_{C}$:

[TABLE]

where $g_{C}(t)$ is a time-dependent coupling constant, and $V_{R}$ and $V_{S}$ are operators of the work resource and of the external system, respectively. In the analytical treatment below, we assume a sufficiently weak coupling between the work resources and the external system justifying a perturbative treatment note:general .

In this Letter, we characterize the work performed by the work resources on the system $S$ by its internal energy change $\Delta U$. A more in-depth discussion as to what extent $\Delta U$ represents work will be provided in the conclusions. $\Delta U$ is evaluated by energy measurements on the external system $S$ at the beginning and the end of the cycle at $t=0$ and $T$, respectively watanabe17 . For simplicity, we turn off the coupling $g_{C}(t)$ at $t=0$ and $T$, and, thus, $[H_{S},H(t)]=0$ at these moments in time. Consequently, measurements of the system energy $H_{S}$ at these two times do not affect the state of the work resources $R$. The external system is initially prepared in its ground state $|0\rangle_{S}$, and the initial state $\rho_{0}$ of the total system is $\rho_{0}=\rho^{R}_{0}\otimes|0\rangle_{S}{{}_{S}}\langle 0|$ with $\rho^{R}_{0}$ being the initial state of the work resources.

Average energy delivered by outcoupled work resources.— We consider the average internal energy change of the system $S$, $\langle\Delta U\rangle_{N}$, caused by the coupling to the $N$ work resources. In the rotating frame with respect to $H_{0}(t)\equiv H_{R}(t)+H_{S}$, the propagator in the interaction picture is $U^{(I)}(t,0)=\mathcal{T}\exp{[-i\int_{0}^{t}H_{C}^{(I)}(t^{\prime})\,dt^{\prime}]}$, where $\mathcal{T}$ is the time-ordering operator and $H_{C}^{(I)}(t)\equiv U_{0}^{\dagger}(t,0)\,H_{C}(t)\,U_{0}(t,0)$ with $U_{0}(t,0)\equiv\mathcal{T}\exp{[-i\int_{0}^{t}H_{0}(t^{\prime})\,dt^{\prime}]}$. Regarding the coupling Hamiltonian $H_{C}$ in Eq. (2), $H_{C}^{(I)}(t)$ reads $H_{C}^{(I)}(t)=g_{C}(t)\,V_{R}^{(I)}(t)\otimes V_{S}^{(I)}(t)$ with $V_{R}^{(I)}(t)\equiv U_{R}^{\dagger}(t,0)\,V_{R}\,U_{R}(t,0)$, $V_{S}^{(I)}(t)\equiv e^{iH_{S}t}\,V_{S}\,e^{-iH_{S}t}$, and $U_{R}(t,0)\equiv\mathcal{T}\exp{[-i\int_{0}^{t}H_{R}(t^{\prime})\,dt^{\prime}]}$.

Setting the energy of the ground state $|0\rangle_{S}$ of the system to be zero without loss of generality, the average internal energy change is given by $\langle\Delta U\rangle_{N}=\sum_{i\neq 0}\epsilon_{i}^{S}\,p_{i}$, where the probability $p_{i}$ for measuring the $i$th eigenvalue $\epsilon_{i}^{S}$ of $H_{S}$ as an outcome of the energy measurement at $t=T$ reads

[TABLE]

Here, $\mathrm{Tr}\,_{R}[\cdots]$ is the trace over the Hilbert space of the work resources and $|i\rangle_{S}$ is the $i$th eigenvector of $H_{S}$. To gain an analytical insight, we first resort to the weak coupling regime where $\int_{0}^{T}g_{C}(t)\,dt\ll 1$. In this limit, expanding the propagator to leading order as $U^{(I)}(T,0)\approx I-i\int_{0}^{T}dt\,g_{C}(t)V_{R}^{(I)}(t)\otimes V_{S}^{(I)}(t)$ in Eq. (3), the excitation probability of the system reduces to

[TABLE]

with $\langle\cdots\rangle_{\rho^{R}_{0}}\equiv\mathrm{Tr}\,_{R}[\cdots\rho^{R}_{0}]$.

Quantum statistical enhancement.— To demonstrate the genuinely quantum mechanical advantage of indistinguishable bosons in comparison to distinguishable particles as the work resources, we consider $N$ QHEs, each performing an Otto cycle with the two lowest internal energy levels of a bosonic atom prepared in its center of mass (COM) ground state as a working fluid, i.e., the temperature $\beta_{\rm COM}^{-1}$ of the COM degrees of freedom is set to be zero. As sketched in Fig. 1, the work resources $R$ contain these engines, together with the hot and cold heat baths.

The four strokes of the Otto cycle are performed as follows: (0) Initial state.— In the absence of the coupling to the external system $S$, $g_{C}(0)=0$, all two-level atoms are prepared in thermal equilibrium with the common cold bath at inverse temperature $\beta_{c}$. Thus, the initial reduced density matrix $\rho^{E}_{0}\equiv\mathrm{Tr}\,_{B}{\,\rho^{R}_{0}}$ of the engine part is $\rho^{E}_{0}=Z_{\beta_{c}}^{-1}\,\exp{[-\beta_{c}H_{E}(0)]}$ with $Z_{\beta_{c}}\equiv\mathrm{Tr}\,_{E}{\,\exp{[-\beta_{c}H_{E}(0)]}}$, where $\mathrm{Tr}\,_{E}[\cdots]$ and $\mathrm{Tr}\,_{B}[\cdots]$ are the trace over the Hilbert space of the engines and that of the baths, respectively. The baths are assumed to be time independent and in the canonical state of $H_{B}$ throughout the cycle. (1) Isentropic compression.— From $0<t<T/2$, all the engines are decoupled from the baths, $H_{EB}=0$, and the level distance of all the two-level atoms is slowly increased in the same manner. (2) Hot isochore.— At $t=T/2$, setting $g_{C}=0$, all the two-level atoms are brought into weak contact with a common hot bath and thermalized at inverse temperature $\beta_{h}$. At the end of this process, the state of the engine is given by $\rho^{E}_{T/2}=\mathrm{Tr}\,_{B}{\,\rho^{R}_{T/2}}=Z_{\beta_{h}}^{-1}\,\exp{[-\beta_{h}H_{E}(T/2)]}$ with $Z_{\beta_{h}}\equiv\mathrm{Tr}\,_{E}{\,\exp{[-\beta_{h}H_{E}(T/2)]}}$. (3) Isentropic expansion.— From $T/2<t<T$, all engines are decoupled from the baths, $H_{EB}=0$, and the energy separation of each two-level atoms is decreased slowly in the same way. (4) Cold isochore.— At $t=T$, setting $g_{C}=0$, and all the two-level atoms are brought into contact with the common cold bath again and quickly return to the initial state.

First, we focus on the case of indistinguishable atoms. We choose $V_{R}=2S_{x}$ in the coupling Hamiltonian (2):

[TABLE]

with $S_{x}\equiv(a^{\dagger}b+b^{\dagger}a)/2$, where $a^{\dagger}$ and $a$ are creation and annihilation operators of the ground-state atoms in the lowest COM level, and $b^{\dagger}$ and $b$ are those of the excited-state atoms, respectively. While we keep the external system general in this discussion, we note that if the external system is a harmonic oscillator (HO) and $V_{S}=c^{\dagger}+c$ with $c$ ($c^{\dagger}$) denoting the annihilation (creation) operator of the HO, Eq. (5) reduces to the standard dipole coupling between an ensemble of atoms and a single-mode HO. In order to compute $\langle\Delta U\rangle_{N}$ using the probability in Eq. (4), we now choose the following engine Hamiltonian:

[TABLE]

with $S_{z}\equiv(a^{\dagger}a-b^{\dagger}b)/2$. In the Otto cycle, both the compression and expansion strokes are done without coupling to the heat baths, and, hence, they are described by unitary dynamics governed by the engine Hamiltonian (6). For quasistatic changes of $\Omega(t)$, the propagator can be written as

[TABLE]

with $\phi(t,t_{0})=\int_{t_{0}}^{t}dt^{\prime}\,2E_{t^{\prime}}$. Here, we denote by $|m,\theta_{t}\rangle_{E}$ the eigenstate of the instantaneous engine Hamiltonian $H_{E}(t)$ with eigenvalue $2E_{t}m$, where $E_{t}\equiv\sqrt{\Omega(t)^{2}+\Delta^{2}}$. Furthermore, $\theta_{t}$ is defined by $\tan{\theta_{t}}=-\Omega(t)/\Delta$. The initial time is $t_{0}=0$ for the isentropic compression and $t_{0}=T/2$ for the isentropic expansion strokes. With this adiabatic propagator, we obtain the autocorrelation function of the operator $V_{R}$ (see Supplement for details), when $t_{0}\leq\{t,t^{\prime}\}\leq t_{0}+T/2$, in Eq. (4) as note:autocorrfunc

[TABLE]

with $\beta_{0}\equiv\beta_{c}$, $\beta_{T/2}\equiv\beta_{h}$, and $F_{\pm}(N,\beta_{t_{0}}E_{t_{0}})\equiv\langle m^{2}\rangle_{t_{0}}\pm\langle m\rangle_{t_{0}}$, where the expectation values are defined with respect to the thermal state of $H_{E}(t_{0})$ at $\beta_{t_{0}}$. On the other hand, if $t$ and $t^{\prime}$ are separated by the thermalization process at $t=T/2$, for instance $t<T/2$ and $t^{\prime}>T/2$, the autocorrelation function takes the factorized form $\langle V_{R}^{(I)}(t^{\prime})\,V_{R}^{(I)}(t)\rangle_{\rho^{R}_{0}}=\langle V_{R}^{(I)}(t^{\prime})\rangle_{\rho^{R}_{T/2}}\langle V_{R}^{(I)}(t)\rangle_{\rho^{R}_{0}}$ with

[TABLE]

Next, we examine the case in which all the atoms are distinguishable. In this case, the coupling Hamiltonian (5) reduces to $H_{C}(t)=g_{C}(t)\,2\sum_{j=1}^{N}(\sigma_{j,\,x}/2)\otimes V_{S}$ and the engine Hamiltonian (6) to $H_{E}(t)=2\Omega(t)\sum_{j=1}^{N}(\sigma_{j,\,z}/2)+2\Delta\sum_{j=1}^{N}(\sigma_{j,\,x}/2)$, where $\sigma_{j,\,x}$ and $\sigma_{j,\,z}$ are the Pauli matrices of the $j$th atom. Assuming quasistatic changes of $\Omega(t)$ like in the indistinguishable case, we find the following autocorrelation function:

[TABLE]

and the average

[TABLE]

allowing the calculation of the probability in Eq. (4) for the distinguishable case. We emphasize that in the distinguishable case, while taking the trace over the engine states, all the possible $2^{N}$ configurations of the atomic pseudospins have to be considered, while in the indistinguishable case the trace is taken over only $N+1$ symmetrized eigenstates of $S_{z}$. Thus, the collective nature of the latter set of states gives rise to an enhanced coupling with the external system and results in the enhancement of $\Delta U$ that we demonstrate next.

To clearly evidence that Bose statistics leads to a quantum advantage, we specialize the coupling protocol to the impulse form $g_{C}(t)=g\delta(t-t_{1})$, with $0<t_{1}<T/2$. In this case, the expressions for the probability (4) are simplified greatly, and the average internal energy change equals

[TABLE]

Thus, for the impulse form of the coupling, $\langle\Delta U\rangle_{N}$ for the indistinguishable and distinguishable cases differ by the value of the variance $\langle[V_{R}^{(I)}(t_{1})]^{2}\rangle_{\rho^{R}_{0}}$ that can be evaluated from Eqs. (8) and (10). Remarkably, we find that, in the indistinguishable case, this variance is larger than or equal to that of the distinguishable case for any values of the parameters $N\geq 1$, $\beta_{c}E_{0}$, and $\theta_{t_{1}}$ (see Supplement ). This fact guarantees that $\langle\Delta U\rangle_{N}^{\text{indist}}$ due to $N$ indistinguishable bosonic engines is always larger than $\langle\Delta U\rangle_{N}^{\text{dist}}$ arising from the same number of distinguishable engines. The resulting enhancement can be quantified by the ratio $\mathcal{E}=\langle\Delta U\rangle_{N}^{\text{indist}}/\langle\Delta U\rangle_{N}^{\text{dist}}$.

In Fig. 2 (a), we compare our analytical result for the enhancement $\mathcal{E}$, for different values of $\theta_{t_{1}}$ and $N$, with numerical simulations for a specific choice of the system as a HO with frequency $\omega$, i.e., $H_{S}=\omega c^{\dagger}c$ and $V_{S}=c^{\dagger}+c$ in the coupling Hamiltonian (5). For the engine Hamiltonian $H_{E}(t)$, we consider linear sweeps of $\Omega(t)$ as $\Omega(t)=\Omega(0)+vt$ and $\Omega(t)=\Omega(0)+v(T-t)$ for the isentropic compression and expansion strokes, respectively. Hereafter, let us focus on the situation with $\Delta=0$ (i.e., $\theta_{t}=-\pi/2$), where the difference between the indistinguishable and distinguishable cases is most prominent note:adiabaticity . When $\Delta=0$, using Eq. (10) we see that $\langle\Delta U\rangle_{N}\propto N$ for the distinguishable case, while in the indistinguishable case using Eq. (8) the dependence on $N,\beta_{c}E_{0}$ is more involved (see Supplement ). In general, although $\langle[V_{R}^{(I)}(t_{1})]^{2}\rangle_{\rho^{R}_{0}}$ contains both terms proportional to $N^{2}$ and $N$, we find that $\langle\Delta U\rangle_{N}$ shows $N^{2}$ scaling for moderate values of $N$ with $N\beta_{c}\Omega(0)\lesssim 1$ note:dicke . We see this behavior in Fig. 2 (b), where we plot $\sqrt{\langle\Delta U\rangle_{N}/\langle\Delta U\rangle_{1}}$ to bring out the quadratic scaling. We also note that, for sufficiently large $N$, the $N^{2}$ scaling of $\langle[V_{R}^{(I)}(t_{1})]^{2}\rangle_{\rho^{R}_{0}}$ for the indistinguishable case turns into a linear scaling with an enhanced slope of $\coth{\beta_{c}\Omega(0)}>1$, while it is unity in the distinguishable case.

Considering general coupling protocols $g_{C}(t)$, when $\Delta=0$, we find that the autocorrelation $\langle V_{R}^{(I)}(t^{\prime})\,V_{R}^{(I)}(t)\rangle_{\rho^{R}_{0}}$ is factorized, and hence, vanishes when $t$ and $t^{\prime}$ are separated by a thermalization step with the hot bath. This allows us to simplify (4) and write the probability of excitation of the driven system in the indistinguishable case as

[TABLE]

and in the distinguishable case as

[TABLE]

where the positive, coupling-protocol-dependent terms are determined by the amplitudes $c^{\pm}_{i}(t_{0})=\int_{t_{0}}^{t_{0}+T/2}dt\,g_{C}(t)\,\,{{}_{S}}\langle i|V_{S}^{(I)}(t)|0\rangle_{S}\,e^{\pm i\phi(t,t_{0})}$. Comparing the terms in brackets in Eqs. (13) and (14), we have that $p_{i}^{\mathrm{indist}}\geq p_{i}^{\mathrm{dist}}$. Thus, for $\Delta=0$, our central result of the enhancement of internal energy change for indistinguishable bosonic engines still holds for arbitrary coupling protocols and external system Hamiltonians.

In order to widen the scope of our results, we also consider engine strokes with $\Delta\neq 0$ as described in more detail in Ref. Supplement . There, we find that the enhancement persists for small values of $N$ independently of the form of the coupling and the external system Hamiltonian. The fact that the enhancement is guaranteed for small $N$ most likely will be of particular relevance to experiments in the near future that will presumably have access to small $N$. Furthermore, we demonstrate for a given choice of the Hamiltonian how to identify generic parameter regimes and coupling protocols that lead to enhancement. We also extend our results to nonperturbative continuous coupling $g_{C}(t)$ [see Supplement for the exact functional form of $g_{C}(t)$] using numerical simulations for a HO external system. From Fig. 3, we see that, for small values of $N$ and moderate values of the coupling strength, there is enhancement and $\langle\Delta U\rangle_{N}^{\rm indist}$ scales as $N^{2}$. The fact that we are able to find enhancement for a generic set of parameters as in Fig. 3 without fine-tuning suggests the general applicability of our result. Finally, we note that the performance in the case of engines made of $N$ identical noninteracting fermionic two-level atoms, considered in Ref. Supplement , is generally diminished for even $N$ and converges to that by a single engine for odd $N$ in the limit of small $\beta_{\rm COM}^{-1}$ due to the Pauli blocking effect.

In conclusion, we have demonstrated that the statistical indistinguishability of work resources can be exploited to gain a genuine quantum enhancement in quantum thermodynamics. While we identify this enhancement in terms of the internal energy change of an external system coupled to the engines, the question arises as to how much of this change is attributable to the actual action of the engines and how much results from the time dependence of the part of the Hamiltonian describing the interaction between the system and the engines. Our preliminary analysis indicates that an accordingly corrected work contribution of engines also displays enhancement in the parameter regimes considered here Supplement . The predicted enhancement of the energy output from multiple indistinguishable heat engines to a generic external system is readily testable with current or near-future experimental realizations of quantum heat engines, e.g., in nitrogen-vacancy centers Klatzow19 , trapped ions vonlindenfels , and ultracold gases. While we have considered bosonic and fermionic statistics Supplement , exotic fractional statistics khare05 ; chen18 may lead to further interesting results.

.1 Performance of multiple engines: derivation of Eqs. (8)–(11)

Let us begin by considering the indistinguishable case. The adiabatic propagator (propagator under quasistatic changes of external parameters) given by Eq. (7) in the main paper has been written in terms of the instantaneous energy eigenstates, which can be expressed as

[TABLE]

with $S_{z}|m\rangle_{z}=m|m\rangle_{z}$. Note that, in what follows, the variable $m$ takes values between $-N/2\leq m\leq N/2$. For ease of presentation, we will not show the range of $m$ in the summations that appear below. Using the above form for the eigenstates, the operator $V_{R}^{(I)}(t)$ in the coupling Hamiltonian can be written as

[TABLE]

with $S_{\alpha}(t_{0})=e^{-i(\theta_{t_{0}}+\frac{\pi}{2})S_{y}}\,S_{\alpha}\,e^{i(\theta_{t_{0}}+\frac{\pi}{2})S_{y}}$, $S_{\pm}\equiv(S_{x}\pm iS_{y})$, and $\phi(t,t_{0})$ being the dynamical phase acquired during the adiabatic evolution from time $t_{0}$ to $t$ given by

[TABLE]

We can now use the rotated collective spin operators to express the initial Hamiltonian of the engine as $H_{E}(t_{0})=2E_{t_{0}}S_{z}(t_{0})$ and the density matrix of the engines at $t_{0}$ is $\rho^{E}_{t_{0}}=e^{-2\beta_{t_{0}}E_{t_{0}}S_{z}(t_{0})}/Z_{\beta_{t_{0}}}$, allowing us to calculate the expectation value of the coupling operator as

[TABLE]

which proves Eq. (9) in the main paper. Note that the partition function in the indistinguishable case is given by $Z_{\beta_{t_{0}}}=e^{N\beta_{t_{0}}E_{t_{0}}}(1-e^{2(N+1)\beta_{t_{0}}E_{t_{0}}})/(1-e^{2\beta_{t_{0}}E_{t_{0}}})$. In order to compute the correlation function given in Eq. (8), consider first the operator product

[TABLE]

with “h.c.” denoting Hermitian conjugate. Since the density matrix at $t_{0}$ is diagonal in the instantaneous basis $|m,\theta_{t_{0}}\rangle$, only the first line of the above equation contributes to the trace. Taking into account the following relations, $S_{+}(t_{0})\,S_{-}(t_{0})=\frac{N}{2}(\frac{N}{2}+1)-[S_{z}^{2}(t_{0})-S_{z}(t_{0})]$ and $S_{-}(t_{0})\,S_{+}(t_{0})=\frac{N}{2}(\frac{N}{2}+1)-[S_{z}^{2}(t_{0})+S_{z}(t_{0})]$, and defining

[TABLE]

we obtain the expression (8) of the main paper for $\langle V_{R}^{I}(t^{\prime})\,V_{R}^{I}(t)\rangle_{\rho^{R}_{t_{0}}}$.

Let us now consider the distinguishable case. As mentioned in the main paper, the time-dependent engine Hamiltonian in this case is given by

[TABLE]

and the resulting propagator in the quasistatic limit is simply given by a direct product over individual engine propagators:

[TABLE]

Here, $|m,{\theta}_{t}\rangle_{j}$ is the instantaneous eigenstate of the $j$th individual engine satisfying $H_{Ej}(t)|m,{\theta}_{t}\rangle_{j}=2E_{t}|m,{\theta}_{t}\rangle_{j}$, which is given by

[TABLE]

with $\sigma_{j,\,z}|m\rangle_{jz}=2m|m\rangle_{jz}$. We can now repeat the calculation in a manner similar to the indistinguishable case. To this end, we first note that the coupling operator in the interaction picture is given by

[TABLE]

with $\sigma_{j,\alpha}(t_{0})=e^{-i(\theta_{t_{0}}+\pi/2)\sigma_{j,\,y}/2}\sigma_{j,\alpha}e^{i(\theta_{t_{0}}+\pi/2)\sigma_{j,\,y}/2}$ and $\sigma_{\pm}\equiv(\sigma_{x}\pm i\sigma_{y})/2$. Using the above expression and the fact that the density matrix for the distinguishable case is given by

[TABLE]

we obtain the average of the coupling operator (noting that only the first term from Eq. (S9) contributes as the density operator is diagonal in the $|m,{\theta}_{t}\rangle_{j}$ basis)

[TABLE]

This leads to Eq. (11) of the main paper since each term in the sum is the same and the trace over everything but the $j$th engine cancels with the partition function in the denominator to give $\mathrm{Tr}\,{[\sigma_{j,\,z}(t_{0})\,\rho^{E}_{t_{0}}]}=-\tanh{(\beta_{t_{0}}E_{t_{0}})}$. As before, we can again write down the two-time correlation operator for the coupling and many of the terms vanish while taking the average due to the diagonal nature of the density matrix. Finally, we get for the correlation function

[TABLE]

Notice that in the first line of the above equation, the first sum has $N$ equal terms and the second sum has $N(N-1)$ equal terms, and the second line of the equation has $N$ equal terms. Combining this with the fact that $\mathrm{Tr}\,{[\sigma_{j,\,+}(t_{0})\,\sigma_{j,\,-}(t_{0})\,\rho^{E}_{t_{0}}]}=e^{-\beta_{t_{0}}E_{t_{0}}}/[2\cosh{(\beta_{t_{0}}E_{t_{0}})}]$ and $\mathrm{Tr}\,{[\sigma_{j,\,-}(t_{0})\,\sigma_{j,\,+}(t_{0})\,\rho^{E}_{t_{0}}]}=e^{\beta_{t_{0}}E_{t_{0}}}/[2\cosh{(\beta_{t_{0}}E_{t_{0}})}]$ leads to Eq. (10) of the main paper.

.2 Impulse-type coupling between the engines and the system

The average internal energy change with impulse-type coupling at time $t_{1}$ (assuming $0<t_{1}<T/2$) is discussed in the main paper, see Eq. (12). This discussion carries over to the case $T/2<t_{1}<T$ up to the fact that the variance has to be taken with respect to the work resource state $\rho^{R}_{T/2}$ and is given by

[TABLE]

An important simplification arising from the impulse coupling becomes evident in the above equation: the average internal energy change is a product of a term that depends only on the engine or work-resource operators and a term that contains the system operators alone. The engine-dependent quantity is the second moment of the coupling $\langle[V_{R}^{(I)}(t_{1})]^{2}\rangle_{\rho^{R}_{0}}$, and it can be obtained using the general expression for the two-time correlators $\langle V_{R}^{(I)}(t)V_{R}^{(I)}(t^{\prime})\rangle_{\rho^{R}_{0}}$ provided in Eqs. (8) and (10) of the main paper as

[TABLE]

for the indistinguishable case and as

[TABLE]

for the distinguishable case. Let us now introduce the function $f(N,\beta_{t_{0}}E_{0})\equiv\langle m^{2}\rangle_{t_{0}}$ to explicitly account for the functional dependence on the expectation value of $m^{2}$ (recall that $t_{0}=0,T/2$ and $\beta_{0}=\beta_{c}$ and $\beta_{T/2}=\beta_{h}$):

[TABLE]

Now, using the inequalities

[TABLE]

we compare the coefficients of the $\sin^{2}{\theta_{t_{1}}}$ and $\cos^{2}{\theta_{t_{1}}}$ terms in Eqs. (S14) and (S15). It follows that the second moment $\langle[V_{R}^{(I)}(t_{1})]^{2}\rangle$, and thus the average internal energy change as well, are larger for the indistinguishable case, as we stated in the main text.

Let us now consider $\Delta=0$ (i.e., $\cos{\theta_{t_{1}}}=0$) and examine the scaling of the average internal energy change with $N$. To this end, note from Eq. (S13) that it is sufficient to examine the scaling of the second moment with $N$. Nonetheless, to make the comparison with numerical solutions as in Fig. 2(b) of the main paper, we note that when the external system is chosen as a harmonic oscillator with frequency $\omega$, one finds ${{}_{S}}\langle i|V_{S}^{(I)}(t)|0\rangle_{S}=e^{i\omega t}\delta_{i,1}$ and the average internal energy change (S13) becomes

[TABLE]

Turning our attention to the second moment for the distinguishable case, we first see that $\langle[V_{R}^{(I)}(t_{1})]^{2}\rangle_{\rho^{R}_{0}}=N$. To examine the behavior of the variance in Eq. (S14) for the indistinguishable case, consider first Eq. (S16) in the limit of $N\rightarrow\infty$

[TABLE]

that yields

[TABLE]

and scales linearly with $N$ with a slope $\coth(\beta_{c}E_{0})\geq 1$, as stated in the main paper. This term is thus larger in the indistinguishable case than in the distinguishable counterpart and is responsible of the quantum statistical enhancement of the average internal energy change. In the other extreme with $N=1$, since $f(N=1,\,\beta_{c}E_{0})=1/4$, the second moment in the indistinguishable and distinguishable cases are both equal to $1$, as expected. For $N\beta_{c}\Omega(0)\sim 1$ or less we have, $\langle[V_{R}^{(I)}(t_{1})]^{2}\rangle_{\rho^{R}_{0}}\simeq f_{1}(\beta_{c}\Omega(0))\,N^{2}+f_{2}(\beta_{c}\Omega(0))\,N$, with $f_{1}(x)\equiv x(x\coth{x}-1)\cosh{x}/\sinh^{3}{x}$ and $f_{2}(x)\equiv 1-[(x\coth{x}-1)/\sinh^{2}{x}]$. Therefore, $\langle\Delta U\rangle_{N}$ shows $N^{2}$ scaling for moderate values of $N$ with $N\beta_{c}\Omega(0)\lesssim 1$ for the indistinguishable case and is larger than the distinguishable case which scales linearly with $N$.

.3 General coupling between the engines and the system

For a general continuous coupling, the probability of exciting the $i$th eigenstate of the system reads

[TABLE]

Using the two-time correlators and one-time averages stated in the main paper (see Eqs. (8)–(11) there), the excitation probability for the case of many indistinguishable and distinguishable engines can be respectively written as

[TABLE]

Here, the coupling-protocol- and engine-state-dependent amplitudes are given by

[TABLE]

and generalize the expressions presented in the main paper. Similar to the function $f(N,x)$ defined in Eq. (S16), we now define $h(N,\beta_{t_{0}}E_{t_{0}})\equiv\langle m\rangle_{t_{0}}$ to account for the parameter dependence of the average $\langle m\rangle$. This gives

[TABLE]

In Eqs. (S23) and (S24), the common protocol-dependent, and $N$-, $\beta_{t_{0}}$-independent parts of the first three terms given by $|d_{i}(t_{0})|^{2}$ and $|\tilde{c}^{\pm}_{i}(t_{0})|^{2}$ are positive but the third term (on the second line of the equations) proportional to $\Re[d_{i}(0)d_{i}^{*}(T/2)]$ is not necessarily positive. As discussed in the main paper, for $\Delta=0$, $\cos{\theta_{t}}=0$ at all $t$ and the non-positive term vanishes. In this case, one can directly compare the probabilities for the indistinguishable and distinguishable cases in Eqs. (S23) and (S24), by examining the relative size of the $N$- and $\beta_{t_{0}}$-dependent factors. To this end, we use the following relation that holds for any $N>1$ and $x>0$:

[TABLE]

This leads to one of the central results given in the main paper: when $\Delta=0$, for any form of the coupling protocol $g_{C}(t)$, $N>1$, system Hamiltonian and value of $\beta_{c}$ and $\beta_{h}$, one finds that $p_{i}^{\mathrm{indist}}>p_{i}^{\mathrm{dist}}$, i.e., the internal energy change done in the indistinguishable case is larger than the distinguishable case.

Let us now consider $\Delta\neq 0$. While in this case the enhancement of internal energy change for indistinguishable engines is not guaranteed in general, in what follows we identify conditions under which enhancement can be expected. We first note that the inequality Eq. (S18) allows us to compare the first term of Eqs. (S23) and (S24) proportional to $|d_{i}(t_{0})|^{2}$ and see it is larger in the indistinguishable case. Secondly, the $N$- and $\beta_{t_{0}}$-dependent positive factors multiplying $2\Re[d_{i}(0)d_{i}^{*}(T/2)]$ satisfy

[TABLE]

for any $N>1$, $\beta_{c}E_{0}$, and $\beta_{h}E_{T/2}$, with the equality holding for $N=1$. Thus, from the above discussion together with the previous considerations leading up to Eq. (S28), it follows that $p_{i}^{\mathrm{indist}}>p_{i}^{\mathrm{dist}}$ which enhances the average internal energy change, for coupling protocols and external system Hamiltonians that satisfy the constraint $\Re[d_{i}(0)d_{i}^{*}(T/2)]\geq 0$. Note that $h(N,\beta_{c}E_{0})\,h(N,\beta_{h}E_{T/2})$ is an increasing function of $N$ with magnitude comparable to the functions on the left-hand side of the inequalities in Eqs. (S18) and (S28) for small $N$. Thus, for small $N$, the magnitude of the negative term when $\Re[d_{i}(0)d_{i}^{*}(T/2)]<0$ is off-set by the positive terms in Eq. (S23) and enhancement persists. In other words, for small enough $N$, enhancement is generally present, regardless of the form of the coupling function, the external system Hamiltonian, and the temperatures of the baths, even when $\Delta\neq 0$.

We will now demonstrate that we can, without much fine tuning, identify parameter regimes where enhancement is possible for a given choice of the external system Hamiltonian and smooth form for the coupling function. Let us choose for the external system, a harmonic oscillator (HO) of frequency $\omega$ coupled to the engines. For the coupling, we consider a function which is constant for most of the duration of the work stroke and turns on and off smoothly. A suitable choice is the one used in our previous work watanabe17suppl

[TABLE]

where $\alpha$ is the switching rate and $\delta_{t}$ $(0<\delta_{t}<1)$ is the ratio of the duration in which the coupling is turned on from $t_{\rm on}+nT/2$ to $t_{\rm off}+nT/2$ $(n=0$ and $1)$ in the duration of the half of the cycle for $T/2$: $t_{\rm on}=(1-\delta_{t})T/4$ and $t_{\rm off}=t_{\rm on}+\delta_{t}T/2$. We note that this is the coupling form used in the numerical calculations generating the results in Figs. 3 and 4 of the main paper. Let us, for the sake of concreteness, assume that the engine compression and expansion strokes are represented by a linear protocol with velocity $v$ as described in the main paper. We would now like to identify the parameter choices for which $\Re[d_{i}(0)d_{i}^{*}(T/2)]>0$ holds and hence results in enhancement. First, we notice that for the HO external system ${{}_{S}}\langle i|V_{S}^{(I)}(t)|0\rangle_{S}=e^{i\omega t}\delta_{i,1}$, i.e. in the weak coupling, perturbative treatment only the first excited state is occupied. Secondly, for simplicity we take the coupling function to have a fast switching rate with $\delta_{t}\approx 1$ allowing the approximation $g_{C}(t)\approx g/T$. Thus, we have

[TABLE]

When the external system is an oscillator with small frequency $\omega T\ll 1$, we can ignore the oscillating phase factor in the integral and find

[TABLE]

leading to

[TABLE]

Thus, in this limit of small oscillator frequency a quantum enhancement is always present. On the other hand when $\omega T\gg 1$,

[TABLE]

and

[TABLE]

with $\Omega(T/2)=\Omega(0)+vT/2$. The presence of the oscillating factors in the above expression can lead to $\Re[d_{i}(0)d_{i}^{*}(T/2)]<0$. While the dependence on $\Delta$, $\Omega(0)$, and $v$ implied by Eq. (S33) is complicated, we can immediately see that when $\Delta\ll\{\Omega(0),\omega\}$ the magnitude of the term is suppressed. Thus even if $\Re[d_{i}(0)d_{i}^{*}(T/2)]<0$, in this limit this term will not be able to suppress the other positive terms in $p_{i}^{\mathrm{indist}}$. Thus we can anticipate that in the small $\Delta$ regime, in agreement with the idea of continuity with the $\Delta=0$ result, enhancement is exhibited. In the other extreme, when $\Omega(0)=0$ and $\Delta\gg vT$ Eq. (S33) reduces to

[TABLE]

In this regime, varying $\omega T>\pi$ thus unveils the regions in which enhancement appears and disappears as a result of the oscillatory nature of the sign of $\Re[d_{i}(0)d_{i}^{*}(T/2)]$. To summarize and support the above discussion, Fig. S1 shows contour plots, as a function of $\Delta/\Omega(0)$ and $\omega T$, of the regions where enhancement is present (grey) or absent (white) for different fixed values of $2\leq N\leq 20$. Here, the expressions (S23) and (S24) are used to calculate the average internal energy change in the indistinguishable and distinguishable case. The main features are in agreement with the analytical arguments above, namely, enhancement is present at small $N$ independent of other parameter choice (for all the parameter space region explored). Further, regions with no enhancement appear for $\omega T>\pi$ at large enough $N$. Such regions also shrink and vanish as $\Delta/\Omega(0)\rightarrow 0$.

.4 Work content of the energy delivered to the system

The increase in the energy content of the system $\Delta U$ is a consequence of the combined action of the engines which convert heat flowing from the hot to the cold reservoir into work, and the time dependence of the coupling constant $g_{C}$ between the engines and the system $S$. In what follows, we provide one method to isolate the work contribution due to the engines alone and, moreover, show that this exhibits collective enhancement for the parameter regimes considered in the main paper.

As a first task, we break the internal energy change (since the external system $S$ initially starts from the ground state with zero initial energy in our setup, this is a property purely of the final state) into ergotropy $\mathcal{W}$ lenard78suppl ; pusz78suppl ; francica17suppl , which is a measure of the maximum work capacity of the final state of the system with respect to its (uncoupled) Hamiltonian, and the average energy $E_{\mathrm{pas}}$ of a passive state that cannot deliver work under cyclic unitary processes, i.e.,

[TABLE]

Thus, ergotropy provides an important energetic measure and quantifies how much of the energy exchanged with the system can be utilized in a later operation with the external system acting as a battery and we will take this as the quantifier of the total work done on the system. In order to identify the fraction of this work contributed purely by turning the coupling to the system on and off, we need to separate this part from the contribution of the action of the engines through the dynamics of their operator $V_{R}^{(I)}(t)$ in the coupling term $g_{C}(t)\,V_{R}\otimes V_{S}$. To evaluate the former contribution, we focus on the external system $S$ part and consider a time-dependent Hamiltonian $H_{\rm driv}(t)\equiv H_{S}+g_{C}(t)\,\bar{V}_{R,\,t}V_{S}$, which is externally driven solely by the time-dependent parameter $g_{C}(t)$. Here, $\bar{V}_{R,\,t}$ is a piecewise-constant factor during each work stroke given by the time average of the work resource coupling operator $V_{R}$ as

[TABLE]

Here, the driving acts on the system via the same operator $V_{S}$ as in Eq. (2) of the main paper and the averages inside the integrals in Eq. (S35) are taken with respect to the state of the engines in thermal equilibrium with the cold ($\rho_{0}^{R}$, compression stroke, $0<t<T/2$) and hot ($\rho_{T/2}^{R}$, expansion stroke, $T/2<t<T$) baths. We identify the average energy change $E_{\mathrm{driv}}\equiv\langle H_{\rm driv}(T)\rangle-\langle H_{\rm driv}(0)\rangle$ of the external system part caused by the unitary process by the time-dependent Hamiltonian $H_{\rm driv}(t)$ as the energy change due to the time dependence of $g_{C}(t)$. This can be a good estimate provided the time dependence of $V_{R}^{(I)}(t)$ is sufficiently weak. Since this unitary process is cyclic due to the form of the coupling $g_{C}(t)$, the energy change (starting from the passive ground state of the system) will also be equal to the ergotropy $\mathcal{W}_{\mathrm{driv}}$ of the external system’s final state, i.e., $E_{\mathrm{driv}}=\mathcal{W}_{\mathrm{driv}}$. With this, we propose that the difference in ergotropies between the engine-driven case and time-dependent Hamiltonian case, $\Delta\mathcal{W}\equiv\mathcal{W}-\mathcal{W}_{\mathrm{driv}}$, be taken as the proxy of the work delivered to the external system from the engine action alone. We will only consider the continuous coupling case for such estimates of $\Delta\mathcal{W}$ since the impulse-type coupling has a nonzero value only at an instance of time and this makes difficult to separate the contribution of the energy injection due to the time-dependent coupling constant from that of the action of the engines. We next present results for the behavior of both $\mathcal{W}$ and $\Delta\mathcal{W}$ for the settings considered in the main paper.

First, we have performed a numerical calculation of the ergotropy of the external system taken as a harmonic oscillator coupled to $N$ distinguishable and indistinguishable bosonic work resources for the situations presented in Figs. 2 and 3 of the main paper. We show the results in Figs. S2 and S3. From this calculation, it is clear that both impulse-type and continuous coupling of the external system to the work resources results in the external system accumulating ergotropy and not merely passive energy. Moreover, in the both cases, the ergotropy displays enhanced values for the indistinguishable bosonic engines case as opposed to the distinguishable one. Additionally, we note that, in the impulse-type coupling case, the ergotropy strongly depends on the value of the variable $\theta_{t_{1}}$ (decided by the time of the kick) but the enhancement looks to be almost independent of this variable.

Finally, for the continuous coupling case, we consider $\Delta\mathcal{W}$, the estimate of the work output due to the engine action we introduced. In Fig. S4, we present the results from a numerical calculation of $\Delta\mathcal{W}$ for the parameters in Fig. 3 (nonperturbative continuous coupling) of the main paper. Here we clearly find that $\Delta\mathcal{W}$ is positive, i.e., the engine contributes to the work content of the system and, moreover, this quantity also exhibits collective enhancement for the parameters considered in the main paper. We have also included another case with $\Delta\neq 0$ shown in Fig. S5, and see the similar behavior as in Fig. S4. We have checked that, as far as $\Delta/\Omega(0)$ is relatively small, this behavior (positive $\Delta\mathcal{W}$ and enhancement) holds for continuous coupling between the engines and the system provided $g$ is sufficiently large or $\omega\sim\Omega(0)$, both of which allow efficient energy transfer from the engines to the external system for $\Delta/\Omega(0)\ll 1$.

.5 Multiple fermionic engines trapped in a harmonic oscillator potential

We discuss the case in which each engine is made of an identical non-interacting two-level fermionic atom. In this case, the average internal energy change $\langle\Delta U\rangle_{N}$ of the system $S$ displays a pronounced even-odd dependence with respect to $N$. When the temperature $\beta_{\rm COM}^{-1}$ of the center of mass (COM) degrees of freedom is zero, $\langle\Delta U\rangle_{N}=0$ for even $N$, and $\langle\Delta U\rangle_{N}=\langle\Delta U\rangle_{1}$ for odd $N$. This is because the atoms inside the Fermi sphere, in which both the ground and excited internal levels are fully occupied, cannot contribute to internal energy change because of the Pauli blocking: In the case of odd $N$, there is one fermion on the top of the Fermi surface, and this is the only atom which can change the internal states in the engine operation. By contrast, there is no such atom for even $N$. Thus in this regime, the internal energy change of the external system due to $N$ work resources is upper-bounded by $\langle\Delta U\rangle_{1}$ and hence in the limit of perturbative coupling, we can anticipate that this quantity will be diminished with respect to the internal energy change for distinguishable engines which would scale linearly with $N$. At nonzero $\beta_{\rm COM}^{-1}$, $\langle\Delta U\rangle_{N}$ takes non-trivial values due to the contribution of atoms near the Fermi surface, where the COM levels are only partially occupied. Moreover, in the regime where $\beta_{\rm COM}^{-1}$ is much smaller than the energy difference $\omega_{\rm trap}$ between the highest occupied and lowest unoccupied trap levels of the atoms, the ratio $\lambda\equiv\langle\Delta U\rangle_{N}/\langle\Delta U\rangle_{1}$ for an arbitrary trap in general can be well characterized by only $\beta_{\rm COM}\omega_{\rm trap}$. In Fig. S6 we show $\lambda$ numerically calculated for multiple fermionic engines in the Otto cycle for several values of $N$. In this example, we assume that atoms are trapped in a HO potential with frequency $\omega_{\rm trap}$ and another HO with frequency $\omega$ is taken as an external system $S$. At higher COM temperatures $\beta_{\rm COM}^{-1}$ for which $\beta_{\rm COM}\omega_{\rm trap}\lesssim 1$, $\langle\Delta U\rangle_{N}$ depends on details of the system such as the spectrum of the COM degrees of freedom determined by the shape of the trapping potential and the Fermi energy. A detailed study of this dependence is beyond the scope of the present work and will be taken up in the future.

We now derive an analytical estimate for the ratio $\lambda$ when $\beta_{\rm COM}^{-1}\omega_{\rm trap}\gg 1$. We will obtain this estimate by considering isolated engines without coupling to an external system as we show below. The Hamiltonian $H_{E}(t)$ of the engines consists of that of the COM degrees of freedom $H_{\rm COM}$ and that of the internal degrees of freedom $H_{\rm in}(t)$, $H_{E}(t)=H_{\rm COM}+H_{\rm in}(t)$, with

[TABLE]

Here, $\omega_{\rm trap}$ is the frequency of the HO trapping potential whose levels are labeled by $l$. The operators $a_{l}^{\dagger}$ ($a_{l}$) and $b_{l}^{\dagger}$ ($b_{l}$) are creation (annihilation) operators of a spin-up and spin-down atom at $l$th level of the trap, respectively, satisfying the fermionic anti-commutation relations $\{a_{l},\,a_{l^{\prime}}^{\dagger}\}=\{b_{l},\,b_{l^{\prime}}^{\dagger}\}=\delta_{l,\,l^{\prime}}$ and $\{a_{l},\,a_{l^{\prime}}\}=\{a_{l}^{\dagger},\,a_{l^{\prime}}^{\dagger}\}=\{b_{l},\,b_{l^{\prime}}\}=\{b_{l}^{\dagger},\,b_{l^{\prime}}^{\dagger}\}=0$ with $\{A,\,B\}\equiv AB+BA$ being the anti-commutator. The distribution of the population $n_{l}\equiv a_{l}^{\dagger}a_{l}+b_{l}^{\dagger}b_{l}$ of the trap levels is conserved, which is set by the COM temperature $\beta_{\rm COM}^{-1}$ and follows the canonical distribution $\propto\exp{(-\beta_{\rm COM}H_{\rm COM})}$.

Using the internal degrees of freedom of the two-level fermionic atoms described by the above Hamiltonian, we consider the quantum heat engines in the Otto cycle. The average of work done by $N$ isolated engines without outcoupling is

[TABLE]

with $\pm\epsilon_{h}\equiv\pm\sqrt{\Delta^{2}+\Omega_{h}^{2}}$ and $\pm\epsilon_{c}\equiv\pm\sqrt{\Delta^{2}+\Omega_{c}^{2}}$ being the energy eigenvalues of the internal degrees of freedom in the hot and cold isochore process, respectively. Here, $f_{N}$ is a function of $\beta_{\rm COM}\omega_{\rm trap}$ whose detailed form varies with $N$. Since $f_{N=1}=1$, the ratio $\lambda$ is

[TABLE]

which is independent of the temperatures of the heat baths and the spectrum of the internal degrees of freedom. Although the exact form of $f_{N}$ non-trivially depends on $N$ in general, in the low COM temperature limit of $\beta_{\rm COM}\omega_{\rm trap}\gg 1$ the form of $f_{N}$ differs only by the parity of $N$ as

[TABLE]

Interestingly, though Eq. (S40) is obtained for the isolated engines without outcoupling, we find that the numerical result (symbols) for $\lambda$ in Fig. (S6) is well fitted by the above analytical formula (dashed lines). Finally, we add that even though Eq. (S40) is obtained for the HO trapping potential, it is applicable for any type of the trapping potential in the low temperature limit, $\beta_{\rm COM}\omega_{\rm trap}\gg 1$, by replacing $\omega_{\rm trap}$ by the energy difference between the highest occupied and the lowest unoccupied levels of the trap at zero COM temperature.