The Wave-Particle Duality in a Quantum Heat Engine

According to the wave-particle duality (WPD), quantum systems show both particle- and wave-like behavior, and cannot be described using only one of these classical concepts. Identifying quantum features that cannot be reproduced by any classical means is key for quantum technology. This task is often pursued by comparing the quantum system of interest to a suitable classical counterpart. However, the WPD implies that a comparison to a single classical model is generally insufficient; at least one wave and one particle model should be considered. Here we exploit this insight and contrast a bosonic quantum heat engine with two classical counterparts, one based on waves and one based on particles. While both classical models reproduce the average output power of the quantum engine, neither reproduces its fluctuations. The wave model fails to capture the vacuum fluctuations while the particle model cannot reproduce bunching to its full extent. We find regimes where wave and particle descriptions agree with the quantum one, as well as a regime where neither classical model is adequate, revealing the role of the WPD in non-equilibrium bosonic transport.

Introduction.-Thewave-particle duality (WPD) expresses the coexistence of particle-and wave-like behavior in a single quantum system [1][2][3].This fundamental principle, confirmed in both photon [4] and matter interferometers [3,5,6], is a pillar of our understanding of quantum mechanics.The WPD has been expressed in quantitative terms [7,8], extended to many-body interference [9], and even tested by interferometric scenarios in which neither particle nor wave models can describe the measurement outcomes [10].
Identifying genuine quantum behavior is central both for quantum technologies [11][12][13][14][15] as well as for fundamental aspects of quantum theory [16][17][18][19][20]. Generally, quantum behavior is identified by comparison with classical models.For instance, quantum computers are benchmarked by classical computers to identify a quantum advantage [21].There is however no general recipe to determine the classical models that serve as a benchmark.Often, only a single classical model is considered, e.g., when results from quantum optics experiments are compared to predictions from classical electrodynamics [22].However, the WPD implies that one model is not enough: To identify genuine quantum behavior, a system of bosons, for instance, should be benchmarked by both classical waves and particles.Otherwise, classical wave or particle phenomena may be misinterpreted as quantum signatures.
In this work, we exploit this insight from the WPD and contrast a minimal model of bosonic transport with two classical counterparts.We focus on a setup where quantum coherence is relevant: a pair of harmonic oscillators coupled coherently to each other and to thermal reservoirs at different temperatures, see Fig. 1.This system implements a quantum heat engine [23]; heat flowing from the hot bath to the cold bath gives rise to power output.We compare the quantum heat engine with two classical models, one in which bosons are modeled as waves (using classical Langevin equations) and one where bosons are modeled as particles (using a classical rate equation).Our classical wave and particle models should not be considered merely as approximations to the quantum model.We consider two bosonic modes, described by the time-dependent Hamiltonian in Eq. ( 1), coupled to two thermal reservoirs.(b) Sketch of the wave model: Two classical waves interact by a χ (2) non-linear crystal, resulting in difference frequency generation.(c) Sketch of the particle model: Particles moving from a hot bath to a cold bath produce work by turning a gear, whereby they lose part of their energy.
Instead, they serve as benchmarks, to identify the departure from classical behavior.Remarkably, both classical models reproduce the average power of the quantum model.However, both fail to reproduce fluctuations around this average.While the wave model cannot correctly describe vacuum fluctuations, the particle model does not result in the same amount of bunching as the quantum model.For both models, there are relevant limits where they accurately describe power fluctuations: The wave model becomes accurate in the high-temperature regime, where vacuum fluctuations do not matter; and the particle model becomes accurate for weak coupling, where transport statistics becomes (bi-directional) Poissonian as well as in the high coupling regime, where the two oscillators effectively behave as a single one [24].In these limits, power fluctuations can be described classically using two distinct models for the different limits.Away from these limits, the output power contains signatures of the WPD as neither waves nor particles can capture its fluctuations.Our results showcase that the WPD is a powerful tool to reveal the non-classical features encoded in out-of-equilibrium quantum systems.
Quantum heat engine.-Weconsider a quantum heat engine composed of two bosonic modes [24][25][26], with frequencies Ω h/c , described by the Hamiltonian, in the Schrödinger picture.We work with units ℏ = k B = 1 and ∆ = Ω h − Ω c .Each system mode is connected to a bath with different temperatures, and the heat flow leads to power output, P(t) = −∂ t H(t) [27], our quantity of interest.The setup is depicted in Fig. 1 (a).This model of a quantum heat engine can be realized in a superconducting circuit architecture [26,[28][29][30].In this case, each mode is provided by an LC resonator, the heat baths by transmission lines, and the coupling between the modes is mediated by a Josephson junction.In this case, power is provided by a super-current against a voltage bias, due to photonassisted Cooper-pair tunneling [26].Another possible implementation of this engine is in optomechanical devices [31,32].
We are interested in the average power in the long-time limit, ⟨P⟩ q , and its zero-frequency noise, with δx := x − ⟨x⟩, and at t = 0 we are at the long-time limit (or steady-state in a suitable rotating frame [33]).We also introduced the subscript "q" to distinguish the quantum averages from the averages of the classical models in the following text.Equation ( 2) is directly connected to the variance of work [33].Henceforth, we refer to it simply as noise.
The reduced system dynamics is described by the Lindblad master equation (LME), with Bose-Einstein occupations nα = (e Ω α /T α − 1) −1 , nh ≥ nc , and super-operators We note that due to the coherent coupling, the LME couples diagonal and off-diagonal elements of the density matrix in the particle number basis.
The LME (3) is equivalent to a set of quantum Langevin equations (QLEs), in the input-output formalism [34] where the thermal baths are captured by input fields, b α,in , and with quantum white noise auto-correlation function, ] = δ αβ δ(t ′ − t) and α, β = h, c.This entails classical white noise and vacuum fluctuations, due to the bosonic algebra of the input fields.Moreover, the dynamics of any product of the ladder operators is computed through (ab) = ȧb + a ḃ.In the long-time limit, the average power reduces to ⟨P⟩ q = g∆⟨N h − N c ⟩ q , with N α = a † α a α and explicitly evaluated from a closed set of equations of motion, ⟨d/dt(a † α a β )⟩ q ; α, β = h, c.The same equations of motion are obtained either from the LME (3) or from applying the QLEs (4) and the white-noise auto-correlation functions (5).Armed with the equations of motion, noise (2) is evaluated by employing the quantum regression theorem and Wick's theorem [35], details can be found in the supplemental material [33].
Wave heat engine.-Ourwave model, sketched in Fig. 1 (b), consists of pair of classical fields, externally driven by two thermal white noise sources.Formally, the model is based on the canonical association between the ladder operators and the complex amplitudes for the classical fields, a α ↔ A α .The classical dynamics is given by classical Langevin equations, Above, the input fields encompass classical white noise, with where we indicate the averages of the wave model with "w".Notably, ξ α,in are scalars and commute; thus, the classical fields, A α , are functions of the random inputs, ξ α,in .A similar wave model has also been considered in [12] for unitary dynamics.We can visualize the wave model in a classical optical setting, see Fig. 1 (b).Two cavities with frequencies Ω α are supplied by thermal fluctuations and a χ (2) crystal amounts to difference frequency generation, producing a power-outputfield with frequency ∆ = Ω h − Ω c [36,37].In this case, the average power is given by ⟨P⟩ w = g∆⟨|A h | 2 − |A c | 2 ⟩ w , and the average is taken with respect to classical white noise (7).The evaluation of power statistics closely follows those of the quantum model in the inputoutput formalism.From Eqs. (6) and the white noise relation (7), we compute ⟨d/dt(A * α A β )⟩ w and a (classical) regression theorem [38] combined with Wick-Isserlis' theorem [39] gives the noise.The procedure is carefully addressed in [33].
Particle heat engine.-Inour particle model, particles may reside on two different sites, as sketched in Fig 1 (c).The occupation numbers of those sites is governed by a classical rate equation where p n h ,n c denotes the joint probability for the occupations n h and n c , the intra-system jump rate is given by Γ I = 4g 2 /(κ h + κ c ), and the rate, Γ 0 n h ,n c is s.t.n h ,n c ṗn h ,n c = 0. We note that the rates for particles entering and leaving the system are the same as in the LME (3).Indeed, for g = 0 the rate equation (8) coincides with the evolution of the diagonal elements of ρ given in Eq. (3).In contrast to the quantum model, transport within the system is described by incoherent jump processes, analogous to the jumps between the system and the baths.As sketched in Fig. 1 (c), a particle moving from hot to cold will turn the gear in the orientation of the arrow and produce the work Ω h −Ω c .An opposite and less likely process is also allowed and would decrease the power output.
We find in the long-time limit, ⟨P⟩ p = Γ I ∆⟨n h −n c ⟩ p , where, ⟨x⟩ p = n h ,n c x(n h , n c )p n h ,n c , resembling the behavior of the quantum model.In order to study the power fluctuations, we apply full counting statistics (FCS) [40] to Eq. ( 8).Concretely, we attach counting fields to intra-system transitions and determine the particle current statistics [33].
Average power and noise.-Wefind that the quantum, as well as the wave and particle models, lead to the same average power, where the last equality illustrates an analogy to the addition of three conductances in series.For the wave model, the equality with the quantum one follows since only normal-ordered operators appear in computing the average power; thus, vacuum fluctuations are irrelevant.For the particle model, we note that the steady-state power can be cast solely in terms of average number operators, ⟨P⟩ q = g∆⟨N h −N c ⟩ q , which are reproduced exactly by the particle model [33].
For each model, we find the power noise, where we wrote our results in terms of equilibrium noise, E, and shot noise, S (S p ) [41].The equilibrium noise, is proportional to the response coefficient in power when a temperature bias is applied, in agreement with the fluctuationdissipation theorem [42].For simplicity, we present shot noise in the case General expressions for the noise can be found in [33].
The wave model reproduces the shot noise of the quantum model, but equilibrium fluctuations are reduced since the terms linear in nα are absent in Eq. (10b).These linear contributions stem from vacuum fluctuations and can be traced back to the quantum white noise auto-correlation function (5).This mismatch is shown in Fig. 2 (a) (green shade).While the wave model fails to capture the dominant contributions of noise at low temperatures, nh , nc ≲ 1, it reproduces the quantum noise at high temperatures nα ≫ 1.
In contrast to the wave model, the particle model captures the equilibrium noise but fails to reproduce the shot noise of the quantum model.Indeed, as for the wave model, the quantum noise is an upper bound for the classical one since (S p − S) ≥ 0. Note that the terms linear in nα , interpreted as vacuum fluctuations so far, are related to detailed balance in the particle model since (n α + 1) = e βΩ α nα , leading to the same equilibrium noise.In Fig. 2 (a) the red-shaded region illustrates the mismatch between particle and quantum model, with a maximum mismatch at g/κ = 1/2(1 + √ 3) 1/2 ≈ 2/3.We observe that for both limits g/κ → 0 and g/κ → ∞ the particle model captures the noise of the quantum model, i.e. ⟨P⟩ q = ⟨P⟩ p and ⟨⟨P 2 ⟩⟩ q = ⟨⟨P 2 ⟩⟩ p .For g/κ α → 0, intersystem transitions provide a bottleneck and transport exhibits bi-directional Poissonian statistics, fully characterized by the rates Γ αβ = Γ I nα (n β + 1) [24], For g/κ α → ∞ the system modes hybridize and effectively behave as a single oscillator in contact with two thermal baths [24,43] 16) To supplement the discussion of noise, we now analyze the Fano factor, F = ⟨⟨P 2 ⟩⟩/(⟨P⟩∆).This quantity is a measure of bunching and connected to intensity correlations [44,45]; e.g., thermal light is (super-)Poissonian, F ≥ 1 (bunched) while the flux of single photons is sub-Poissonian (anti-bunched) , F < 1.For the quantum and particle models, we find F q ≥ F p ≥ 1, while the wave model can attain sub-Poissonian statistics since F w can be both smaller and larger than one, as exemplified in Figs. 2 (b/c).To understand these results, we introduce δF qw(p) = F q − F w(p) , which represent the shaded regions in Fig. 2 (b/c).
For the wave model, Eqs.(10a, 10b) give δF qw = (n h + nc )/(n h − nc ).The divergence at equilibrium, nh = nc , is due to the absence of power, while equilibrium fluctuations are still present.We also note that in the high-temperature limit δF qw does not vanish, but F q(w) are dominated by the quadratic terms in nα such that F q ≈ F w , as we see in Fig. 2(c).At low values of g, the wave model results in anti-bunching as F w drops below one when nh nc < nh − nc .In this regime, vacuum fluctuations are crucial to capture the correct statistics, which never show anti-bunching in this quantum model.
Jointly analyzing quantum, wave, and particle models, we conclude that for low temperatures, nα ≲ 1, and g/κ ≈ 2/3 neither wave nor particle models capture the quantum, bosonic, noise encoded in power statistics.Therefore, even a minimal quantum heat engine, once operated in the quantum regime contains complementary equilibrium and nonequilibrium effects stemming from wave-like and particle-like behavior.As we have shown, this is however not conflicting with two possible classical pictures emerging in different parameter regimes.
Thermodynamic Uncertainty Relations (TURs).-Instochastic thermodynamics, the trade-off between power and noise has a well-established bound in terms of the entropy production rate, σ [46][47][48].In contrast to fermionic systems, where the effect of quantum coherence can decrease noise [15], the so-called TUR, ⟨⟨P 2 ⟩⟩/⟨P⟩ 2 ≥ 2/ σ, cannot be violated in our bosonic model [49] and this bound immediately applies to the rate equation (8).For the power of our heat engine, the TUR is equivalent to a bound on the Fano factor, solid gray in Fig. 2(c)].The wave model violates this bound and the violation coincides with the spurious anti-bunching.The reason for TUR violations in the wave model is the choice of the classical white noise in Eq. (7).In classical Langevin equations, the strength of the white noise is given by k B T α instead of nα , which only coincide for large temperatures.Since classical Langevin equations obey the TUR [50], a modified bound holds for the wave model, Alternative classical models.-Theparticle and wave models aforementioned are in principle not unique.In the particle model, for instance, we chose Γ I = 4g 2 /(κ h + κ c ).This choice is motivated because it is the only one that reproduces the quantum average power.For the wave model, a different value for the strength of the white noise [c.f.Eq. ( 7)] could be chosen.Indeed, setting ξ * α,in (t ′ )ξ β,in (t) w = (n α +C)δ αβ δ(t ′ −t), with the same C for α = h, c leaves the average power (9) unchanged.Tuning C we can attempt to account for vacuum fluctuations in the wave model; this has the consequence of modifying the equilibrium part of Eq. (10b) as n2 α → (n α +C) 2 , which should be compared to nα (n α +1) in the quantum model, see Eq. (10a).For C = 1/2, the modified wave model captures the linear terms in nα present in Eq. (10a).However, for any C 0, the modified wave model predicts noise at nh = nc = 0 and thus cannot correctly reproduce vacuum fluctuations for arbitrary temperatures.
Conclusions and outlook.-Wehave shown that the waveparticle duality (WPD) plays a fundamental role in the power statistics of quantum heat engines.We considered a minimal model where, despite the presence of quantum coherence and vacuum fluctuations, two classical descriptions based on either particles or waves reproduce the average power of the quantum model.Power fluctuations, however, contain contributions from vacuum fluctuations and coherence which cannot be reproduced by our wave and particle models respectively.Our work thus highlights the connection between power statistics and the WPD, a cornerstone of quantum theory.Thereby, we provide a novel perspective for understanding engines in the quantum regime.
We stress that our approach of comparing a quantum model to a wave and a particle model may readily be extended to different systems and thereby opens up a novel avenue for determining non-classical behavior.For instance, quantum few-level systems and qubits can be contrasted to either classical few-level systems (particle-like models) or to classical magnets with oscillating magnetization (wave-like models).Furthermore, by considering the full Josephson interaction in a circuit QED implementation of the heat engine considered here [26], the WPD can be exploited in a richer model, which contains squeezing and non-Gaussian effects in the power statistics.
Supplemental Material: The Wave-Particle Duality in a Quantum Heat Engine Marcelo Janovitch, Matteo Brunelli, Patrick P. Potts Department of Physics and Swiss Nanoscience Institute, University of Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland In this Supplemental Material, we provide detailed calculations for quantum, wave and particle models, and a short note on the connection between noise and work variance.Here, for simplicity, we drop the subscripts "q/w/p" used in the main text up to the final results and we also often work with the current, I, connected to power through a multiplicative constant, P = I × ∆.

References 5
Quantum This quantum system is Gaussian, thus the average power and noise can be computed from a set of four equations of motion.As discussed in the main text, there are two ways of finding such equations of motion, and each provides a different insight into our classical models.In this Section, we recast the dynamics of the quantum model in a time-independent fashion and derive equations of motion through a Lindblad master equation and through quantum Langevin equations.Next, the equations of motion are used to produce the average current and its zero-frequency noise.For the quantum model, we do not discuss how to define statistics of power, see [24] for a detailed account of full counting statistics in this model.Here, we start from Eq. ( 2) and calculate it from the quantum regression theorem [35].

Lindblad Master Equation (LME)
Equation ( 1) of the main text can be recast in a rotating frame with respect to where we used the fact that ∆ = Ω h − Ω c and the local master equation (3) retains the same form.We now work with the adjoint Liouvillian, L ; and defined through picture equivalence tr(AL(ρ)) = tr L(A)ρ .From the above expression, we obtain a closed set of equations of motion for the operators (a † h a h , a † c a c , a † h a c , a † c a h ) =: Θ, which can be cast in the form LΘ = XΘ + Y: Therefore, the equations of motion for the averages are d/dt⟨Θ⟩ = X⟨Θ⟩ + Y.
Quantum Langevin Equations for the quantum model (QLEs) We now discuss the connection between the above Lindblad master equation and the QLE (4), and derive the same equations of motion d/dt⟨Θ⟩ = X⟨Θ⟩ + Y.We start by recasting the QLE's of the main text (4) in a rotating frame generated by r , where b α (ω) are bosonic field operators of the baths.In this frame, we can eliminate the time-dependence and the QLEs write, With the above equation, we can again compute the equations of motion for Θ = (a † h a h , a † c a c , a † h a c , a † c a h ), using Leibniz's rule.After some algebra, we obtain Θ Unlike the LME, the QLE considers both environment and system degrees of freedom; the meaning of averaging here is thus inherently different from the one in the LME approach but must produce the same predictions.That is, it must only generate the same equations of motion of the reduced system dynamics at the level of expectation values.
We then concentrate on computing averages, d/dt ⟨Θ⟩; note that the term proportional to ⟨Θ⟩ is already X, the homogeneous term in Eq. (S3), and it remains to be shown that ⟨y⟩ = Y, corresponding to the non-homogeneous term in Eq. (S3).This is done by exploiting the quantum white noise relations (5).We thus need to eliminate a ( †)  α terms in favor of input modes.And for that, we formally solve the equations of motion for the ladder operators (S4) in Fourier space, a α (t) Solving for the system modes we obtain the following relations, We now consider, for concreteness, the first entry of ⟨y⟩, with each term written as a Fourier transform, where, in the last passage, we considered the quantum white noise relations in Fourier space, b † α,in (ν)b β,in (−ν ′ ) = nα δ αβ δ(ν−ν ′ ).The last integral can be evaluated by finding the poles of the integrand, We note that the poles are complex and we can evaluate the integral through a counter-clockwise contour in the upper-half complex plane.Thus, from the residue theorem, We are then left with, Similarly, The remaining entries of ⟨y⟩ are integrals of the form where we again use quantum white noise relations and the same complex-contour integration.Thus, we have shown that ⟨y⟩ = Y.Therefore, the equations of motion for the averages produced by the quantum Langevin approach (S5) are the same as the LME approach (S3).From this point on, we use the equations of motion to obtain steady-state power and noise.

Average current and zero-frequency noise
We start by recasting the equations of motion in a new basis vector of operators Θ → σ = (I, H, N h , N c ) where I = ig(a † h a c − a † c a h ) is the current operator in the new frame, H is Eq.(S1) and N α = a † α a α .The equations of motion are now in the form d/dtσ = Gσ + F with and F = (0, 0, nh κ h , nc κ c ).
We readily obtain steady-state averages by solving ⟨σ⟩ = −G −1 F. In particular, we obtain the average current, Setting κ h = κ c = κ we obtain Eq. ( 9) of the main text.We now turn to the evaluation of the zero-frequency cumulant given by Eq. ( 2).First, we introduce δσ = σ − ⟨σ⟩ which evolves according to the homogeneous equation, and consider the correlator ⟨δσδI⟩, whose first entry is ⟨δIδI⟩.The above homogeneous equation satisfies the requirements of the quantum regression theorem, which implies that ⟨δσ(t)δI(0)⟩ = e tG [δσ(0)]δI(0) .(S19) Integrating the above as in Eq.( 2), we have whose first entry is the desired cumulant.Given the propagator in Eq. (S16), we just need to compute the initial conditions, For instance, to compute I 2 we need four-point functions, due to Gaussianity, we can compute this correlator through Wick contractions, e.g.
where we also used a ( †) µ a ( †) ν = 0. To compute the desired initial conditions, all the following contractions are required: At the instance of Eq. (S22), we get by applying the above, The steady-state solutions for the covariance vector ⟨Θ⟩ fully determine ⟨I 2 ⟩.Similarly, for the other entries of the initial conditions, By inverting G and applying to ⟨δσδI⟩ we obtain a vector, whose first entry is the current noise, 3  .

WAVE MODEL
The calculations for the wave model follow very closely the calculations for the quantum model in the quantum Langevin equation formalism.However, from the beginning, we highlight that the central difference is that, since we no longer have operators, the role of commutation relations is lost in the computation of any correlation function.As we show now, this fact has no implication at the level of the equations of motion (S3), due to normal ordering.Technically, this is why we get the same average current.However, noise is computed using non-normal ordered operators -as can be seen from the initial conditions such as Eq.(S23).For clarity, we reproduce the explicit calculation, departing from the classical Langevin equations (already in the time-independent coordinate system, similar to Eqs.(S4)) From the above we obtain the equations of motion for and, in particular, for the averages Finally, we simplify the terms involving noise, again, by substituting the Fourier transforms and solving the same integrals in the complex plane.We get the equations of motion in the form d/dt⟨Θ⟩ where we have considered a generalized form of classical white noise ⟨ξ * α,in (t ′ )ξ β,in (t)⟩ = Φα δ αβ δ(t ′ − t) to account for both the first wave model introduced in the main text and the discussion of alternative classical models.As in the quantum model, we introduce the new variables σ = (I, H, N h , N c ), where H = g(A * h A c + A * c A c ), is the classical Hamiltonian, N α = A * α A α and the wave current, With this, d/dt⟨σ⟩ = G⟨σ⟩ + F ′ , where we emphasize that coefficient matrix G is the same as the quantum model's G, Eq.(S16) and F ′ = (0, 0, Φh κ h , Φc κ c ).We readily obtain the wave current, from ⟨σ⟩ = −GF ′ , which reduces to the current in the quantum model for Φα = nα + C. Turning to wave noise, we now combine the regression theorem, valid for any Markovian process satisfying a linear set of equations of motion [38] with the Isserlis'(Wick's probability) theorem [39] to compute the current fluctuations.These are similar to the quantum regression theorem and Wick's theorem used before.From the regression theorem, we have, where G is again the coefficient matrix of the equations of motion for σ.It remains to compute the initial conditions, Inspecting I 2 we realize that we need four-point functions, due to Gaussianity, we can compute these correlators through Isserlis contractions, e.g.
In general, all the following contractions will be needed In contrast to Eqs. (S24), above, all "1"-factors are missing; they emerge in the quantum model due to the bosonic commutation rules.Here, the A α are random variables and commute.In summary, in the wave model, we have the same propagator as for the quantum model but different initial conditions; this leads to Finally, using ⟨⟨P 2 ⟩⟩ = ∆ 2 ⟨⟨I⟩⟩ and writing the above in terms of equilibrium and shot noise, where χ = [(κ h + κ c )(4g 2 + κ h κ c )] −1 .Setting κ h = κ c = κ, Φα = nα and using ⟨⟨P 2 ⟩⟩ = ∆ 2 ⟨⟨I⟩⟩ we obtain Eq. (10b) of the main text.The discussion for alternative classical models can be supplemented by taking Φα = nα + C.

PARTICLE MODEL
We start by recalling Eq. ( 8) of the main text, For convenience, we can rewrite the above in matrix form; we start by introducing the vector |p) = n h ,n c p n h ,n c |n h , n c ) and the stochastic operator L , s.t.
As in the previous models, we introduce new variables N h =: (V + I )/(2Γ I ), N c =: (V − I )/(2Γ I ) and write the equations of motion in vector form σ = (I , V ) as where, and F = Γ I (κ h nh − κ c nc , κ h nh + κ c nc ).Again, we have a linear system for the steady-state averages G⟨σ⟩ + F = 0, and ⟨I ⟩ = ⟨I⟩, is given by Using Γ I = 4g 2 /(κ h + κ c ) and that ⟨P⟩ = ⟨I⟩∆, we obtain Eq. ( 9) of the main text.

Noise
We now concentrate in computing the zero-frequency noise.Using Eq. (S59) the first term in Eq. (S57) gives after some algebra, Above we have terms that can be readily computed from the equations of motion for the first moments ⟨N α ⟩, solved in the latter Subsection.However, we also have higher order terms, which demand addressing the equation of motion of ⟨N h N c ⟩. Before we compute the higher moments, we will show a regression theorem that also reduces the remaining term in Eq. (S57) to computing the steady-state averages of ⟨N 2 α ⟩, ⟨N h N c ⟩.We start by rewriting Eq. (S63) as an homogeneous equation, where δσ = σ − ⟨σ⟩.We now consider the LHS of the above, and that | ṗ) = L |p), We now use the above result.In particular, we are interested in the case that A = W 1 .Recalling the second term in Eq. (S57), where we used that the Drazin inverse can be written as L D = − ∞ 0 e tL (1 − |p * ) (1|), with |p * ) the eigenvector with zeroeigenvalue (steady-state) of L .Finally, we can use the regression theorem for ⟨δσe tL δW 1 ⟩ to evaluate the correlator, whose first entry is precisely Eq. (S72).Thus, we obtain the relation, Thus, in order to evaluate Eq. (S72), we need the initial conditions ⟨δσδW 1 ⟩, where ⟨W 1 ⟩ = ⟨I⟩ and ⟨V ⟩ have been computed already.Turning to the remaining terms, which we solve for the steady-state average by setting the LHS to zero.We readily see that these equations are coupled to each other and to the lower moments calculated previously.We can thus solve them and obtain all the steady-state second moments; together with the inverse of G, Eq. (S64), we have obtained all the ingredients to evaluate the particle current, Eq. (S74), and we finally obtain, , where we used the notation nT = nh + nc , δn = nh − nc and substituted Γ I = 4g 2 /(κ h + κ c ).Finally, considering ⟨⟨P 2 ⟩⟩ = ∆ 2 ⟨⟨I⟩⟩ and writing the above in terms of equilibrium and shot noise, we obtain (power) noise, P 2 p = E(n h (n h + 1) + nc (n c + 1)) − S p (n h − nc ) 2 , (S82) where χ = [(κ h + κ c )(4g 2 + κ h κ c )] −1 and S is given in Eq. (S30).Setting κ h = κ c = κ we obtain Eq. (10c) of the main text.
We also mention two important limits discussed in the main text.Expanding the second term in the RHS of Eq. (S84) around g = 0 and 1/g = 0, we have

FIG. 1 .
FIG. 1.(a) Schematic representation of the quantum heat engine.We consider two bosonic modes, described by the time-dependent Hamiltonian in Eq. (1), coupled to two thermal reservoirs.(b) Sketch of the wave model: Two classical waves interact by a χ (2) non-linear crystal, resulting in difference frequency generation.(c) Sketch of the particle model: Particles moving from a hot bath to a cold bath produce work by turning a gear, whereby they lose part of their energy.

FIG. 2 .
FIG.2.Noise and Fano factor.In all plots we set nc = 0.1, κ h = κ c = κ.(a) (linear-log) Power noise as a function of coupling between system modes, g/κ, and average power in the inset.We fixed nh = 2.The green (red)-shaded region indicates the mismatch between wave (particle) and quantum models; the particle model matches the quantum model for both small and large g/κ.(b) (linear-log) Fano factor as a function of g/κ, with nh = 2 (solid) and nh = 10 (dashed).For these temperatures, both classical models show a clear mismatch with the quantum one in the region where g is of the order of κ.The gray dashed line indicates F = 1.(c) (log-log) Fano factor as a function of nh , with g/κ = 2/3 (solid) and g/κ = 10 (dashed).We see that for either coupling regime, the wave model converges to the quantum model for large nh .The gray dotted line indicates F = 1, the solid gray line is the bound provided by the TUR, and the purple line is the bound from the modified TUR which holds for the wave model.We observe sub-Poissonian statistics for the wave model close for nh ≲ 1.