# The Wave-Particle Duality in a Quantum Heat Engine

Marcelo Janovitch<EMAIL_ADDRESS>Matteo Brunelli  Patrick P. Potts
<EMAIL_ADDRESS>Department of Physics and Swiss Nanoscience Institute,
University of Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland

###### Abstract

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.–The wave-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
[5, 6, 3], 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.

(a)(b)(c) Figure 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 $\chi^{(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.

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. 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.–We consider a quantum heat engine composed of two bosonic
modes [25, 26, 24], with frequencies $\Omega_{h/c}$, described by the
Hamiltonian,

$H(t)=\sum_{\alpha=h,c}\Omega_{\alpha}a^{\dagger}_{\alpha}a_{\alpha}+g\quantity(a_{h}^{\dagger}a_{c}e^{-i\Delta
t}+a^{\dagger}_{c}a_{h}e^{i\Delta t}),$ (1)

in the Schrödinger picture. We work with units $\hbar=k_{B}=1$ and
$\Delta=\Omega_{h}-\Omega_{c}$. Each system mode is connected to a bath with
different temperatures, and the heat flow leads to power output,
$P(t)=-\partial_{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 photon-assisted
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, $\langle
P\rangle_{q}$, and its zero-frequency noise,

$\langle\langle P^{2}\rangle\rangle_{q}=2\real\int_{0}^{\infty}\differential
t\expectationvalue{\delta P(t)\delta P(0)}_{q},$ (2)

with $\delta x:=x-\expectationvalue{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),

$\dot{\rho}=-i[H(t),\rho]+\sum_{\alpha=h,c}(\bar{\text{n}}_{\alpha}+1)\kappa_{\alpha}D[a_{\alpha}]\rho+\bar{\text{n}}_{\alpha}\kappa_{\alpha}D[a_{\alpha}^{\dagger}]\rho,$
(3)

with Bose-Einstein occupations
$\bar{\text{n}}_{\alpha}=(e^{\Omega_{\alpha}/T_{\alpha}}-1)^{-1}$,
$\bar{\text{n}}_{h}\geq\bar{\text{n}}_{c}$, and super-operators
$D[L]\rho=L\rho L^{\dagger}-\frac{1}{2}\\{L^{\dagger}L,\rho\\}$. 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]

$\displaystyle\dot{a}_{h}$
$\displaystyle=-\quantity(i\Omega_{h}+\frac{\kappa_{h}}{2})a_{h}-iga_{c}e^{-it\Delta}-\sqrt{\kappa_{h}}b_{h,\text{in}},$
(4a) $\displaystyle\dot{a}_{c}$
$\displaystyle=-\quantity(i\Omega_{c}+\frac{\kappa_{c}}{2})a_{c}-iga_{h}e^{+it\Delta}-\sqrt{\kappa_{c}}b_{c,\text{in}},$
(4b)

where the thermal baths are captured by input fields, $b_{\alpha,\text{in}}$,
and with quantum white noise auto-correlation function,

$\displaystyle\langle
b_{\alpha,\text{in}}^{\dagger}(t^{\prime})b_{\beta,\text{in}}(t)\rangle_{q}$
$\displaystyle=\bar{\text{n}}_{\alpha}\delta_{\alpha\beta}\delta(t^{\prime}-t),$
(5)

with
$[b_{\alpha,\text{in}}(t^{\prime}),b_{\beta,\text{in}}^{\dagger}(t)]=\delta_{\alpha\beta}\delta(t^{\prime}-t)$
and $\alpha,\beta=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
$\dot{(ab)}=\dot{a}b+a\dot{b}$.

In the long-time limit, the average power reduces to $\langle
P\rangle_{q}=g\Delta\langle N_{h}-N_{c}\rangle_{q}$, with
$N_{\alpha}=a^{\dagger}_{\alpha}a_{\alpha}$ and explicitly evaluated from a
closed set of equations of motion, $\langle\differential/\differential
t(a_{\alpha}^{\dagger}a_{\beta})\rangle_{q};~{}\alpha,\beta=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.– Our wave 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_{\alpha}\leftrightarrow A_{\alpha}$. The classical dynamics is given by
classical Langevin equations,

$\displaystyle\dot{A}_{h}$
$\displaystyle=-\quantity(i\Omega_{h}+\frac{\kappa_{h}}{2})A_{h}-igA_{c}e^{-it\Delta}-\sqrt{\kappa_{h}}\xi_{h,\text{in}},$
(6a) $\displaystyle\dot{A}_{c}$
$\displaystyle=-\quantity(i\Omega_{c}+\frac{\kappa_{c}}{2})A_{c}-igA_{h}e^{+it\Delta}-\sqrt{\kappa_{c}}\xi_{c,\text{in}}.$
(6b)

Above, the input fields encompass classical white noise, with

$\expectationvalue{\xi_{\alpha,\text{in}}^{*}(t^{\prime})\xi_{\beta,\text{in}}(t)}_{w}=\bar{\text{n}}_{\alpha}\delta_{\alpha\beta}\delta(t^{\prime}-t),$
(7)

where we indicate the averages of the wave model with “$w$”. Notably,
$\xi_{\alpha,\text{in}}$ are scalars and commute; thus, the classical fields,
$A_{\alpha}$, are functions of the random inputs, $\xi_{\alpha,\text{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 $\Omega_{\alpha}$ are supplied by thermal
fluctuations and a $\chi^{(2)}$ crystal amounts to difference frequency
generation, producing a power-output-field with frequency
$\Delta=\Omega_{h}-\Omega_{c}$ [36, 37].

In this case, the average power is given by $\langle
P\rangle_{w}=g\Delta\langle|A_{h}|^{2}-|A_{c}|^{2}\rangle_{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 input-
output formalism. From Eqs. (6) and the white noise relation (7), we compute
$\langle\differential/\differential t(A_{\alpha}^{*}A_{\beta})\rangle_{w}$ and
a (classical) regression theorem [38] combined with Wick-Isserlis’ theorem
[39] gives the noise. The procedure is carefully addressed in [33].

(a)(b)(c) Figure 2: Noise and Fano factor. In all plots we set
$\bar{\text{n}}_{c}=0.1$, $\kappa_{h}=\kappa_{c}=\kappa$. (a) (linear-log)
Power noise as a function of coupling between system modes, $g/\kappa$, and
average power in the inset. We fixed $\bar{\text{n}}_{h}=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/\kappa$. (b) (linear-log) Fano factor as a function of $g/\kappa$, with
$\bar{\text{n}}_{h}=2$ (solid) and $\bar{\text{n}}_{h}=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 $\kappa$. The gray dashed line
indicates $\mathcal{F}=1$. (c) (log-log) Fano factor as a function of
$\bar{\text{n}}_{h}$, with $g/\kappa=2/3$ (solid) and $g/\kappa=10$ (dashed).
We see that for either coupling regime, the wave model converges to the
quantum model for large $\bar{\text{n}}_{h}$. The gray dotted line indicates
$\mathcal{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
$\bar{\text{n}}_{h}\lesssim 1$.

Particle heat engine.– In our 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

$\displaystyle\dot{p}_{n_{h},n_{c}}$
$\displaystyle=\kappa_{h}(\bar{\text{n}}_{h}+1)(n_{h}+1)p_{n_{h}+1,n_{c}}+\kappa_{h}\bar{\text{n}}_{h}n_{h}p_{n_{h}-1,n_{c}}$
(8)
$\displaystyle+\kappa_{c}(\bar{\text{n}}_{c}+1)(n_{c}+1)p_{n_{h},n_{c}+1}+\kappa_{c}\bar{\text{n}}_{c}n_{c}p_{n_{h},n_{c}-1}$
$\displaystyle+\Gamma_{I}(n_{h}+1)n_{c}p_{n_{h}+1,n_{c}-1}+\Gamma_{I}(n_{c}+1)n_{h}p_{n_{h}-1,n_{c}+1}$
$\displaystyle-\Gamma^{0}_{n_{h},n_{c}}p_{n_{h},n_{c}},$

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
$\Gamma_{I}=4g^{2}/(\kappa_{h}+\kappa_{c})$, and the rate,
$\Gamma^{0}_{n_{h},n_{c}}$ is s.t.
$\sum_{n_{h},n_{c}}\dot{p}_{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 $\rho$ 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 $\Omega_{h}-\Omega_{c}$. An opposite and
less likely process is also allowed and would decrease the power output.

We find in the long-time limit, $\langle P\rangle_{p}=\Gamma_{I}\Delta\langle
n_{h}-n_{c}\rangle_{p}$, where, $\langle
x\rangle_{p}=\sum_{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.– We find that the quantum, as well as the wave and
particle models, lead to the same average power,

$\displaystyle\langle P\rangle_{q}=\langle P\rangle_{w}=\langle P\rangle_{p}$
$\displaystyle=\frac{4g^{2}\kappa_{h}\kappa_{c}\Delta(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})}{(4g^{2}+\kappa_{h}\kappa_{c})(\kappa_{h}+\kappa_{c})}$
(9)
$\displaystyle=\Delta(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})(\kappa_{h}^{-1}+\kappa_{c}^{-1}+\Gamma_{I}^{-1})^{-1},$

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, $\langle P\rangle_{q}=g\Delta\langle N_{h}-N_{c}\rangle_{q}$, which
are reproduced exactly by the particle model [33].

For each model, we find the power noise,

$\displaystyle\langle\langle P^{2}\rangle\rangle_{q}$
$\displaystyle=\mathcal{E}[\bar{\text{n}}_{h}(\bar{\text{n}}_{h}+1)+\bar{\text{n}}_{c}(\bar{\text{n}}_{c}+1)]-\mathcal{S}(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})^{2},$
(10a) $\displaystyle\langle\langle P^{2}\rangle\rangle_{w}$
$\displaystyle=\mathcal{E}(\bar{\text{n}}_{h}^{2}+\bar{\text{n}}_{c}^{2})-\mathcal{S}(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})^{2},$
(10b) $\displaystyle\langle\langle P^{2}\rangle\rangle_{p}$
$\displaystyle=\mathcal{E}[\bar{\text{n}}_{h}(\bar{\text{n}}_{h}+1)+\bar{\text{n}}_{c}(\bar{\text{n}}_{c}+1)]-\mathcal{S}_{p}(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})^{2},$
(10c)

where we wrote our results in terms of equilibrium noise, $\mathcal{E},$ and
shot noise, $\mathcal{S}$ ($\mathcal{S}_{p}$) [41]. The equilibrium noise,

$\displaystyle\mathcal{E}$ $\displaystyle=\frac{\langle
P\rangle\Delta}{\bar{\text{n}}_{h}-\bar{\text{n}}_{c}},$ (11)

is proportional to the response coefficient in power when a temperature bias
is applied, in agreement with the fluctuation-dissipation theorem [42]. For
simplicity, we present shot noise in the case $\kappa_{h}=\kappa_{c}=\kappa$,

$\displaystyle\mathcal{S}$
$\displaystyle=\mathcal{E}\quantity[1-2g^{2}\frac{(4g^{2}+5\kappa^{2})}{(4g^{2}+\kappa^{2})^{2}}],$
(12) $\displaystyle\mathcal{S}_{p}$
$\displaystyle=\mathcal{S}+\mathcal{E}\frac{24g^{4}\kappa^{2}}{(6g^{2}+\kappa^{2})(4g^{2}+\kappa^{2})^{2}}.$
(13)

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 $\bar{\text{n}}_{\alpha}$
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, $\bar{\text{n}}_{h},~{}\bar{\text{n}}_{c}\lesssim 1$, it
reproduces the quantum noise at high temperatures $\bar{\text{n}}_{\alpha}\gg
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 $(\mathcal{S}_{p}-\mathcal{S})\geq 0$. Note that the terms linear in
$\bar{n}_{\alpha}$, interpreted as vacuum fluctuations so far, are related to
detailed balance in the particle model since
$(\bar{n}_{\alpha}+1)=e^{\beta\Omega_{\alpha}}\bar{n}_{\alpha}$, 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/\kappa=1/2(1+\sqrt{3})^{1/2}\approx 2/3$. We observe that for both limits
$g/\kappa\to 0$ and $g/\kappa\to\infty$ the particle model captures the noise
of the quantum model, i.e. $\langle P\rangle_{q}=\langle P\rangle_{p}$ and
$\langle\langle P^{2}\rangle\rangle_{q}=\langle\langle
P^{2}\rangle\rangle_{p}$. For $g/\kappa_{\alpha}\to 0$, inter-system
transitions provide a bottleneck and transport exhibits bi-directional
Poissonian statistics, fully characterized by the rates
$\Gamma_{\alpha\beta}=\Gamma_{I}\bar{\text{n}}_{\alpha}(\bar{\text{n}}_{\beta}+1)$
[24],

$\langle P\rangle_{p}=\Delta(\Gamma_{hc}-\Gamma_{ch}),\hskip
14.22636pt\langle\langle
P^{2}\rangle\rangle_{p}=\Delta^{2}(\Gamma_{hc}+\Gamma_{ch}).$ (14)

For $g/\kappa_{\alpha}\to\infty$ the system modes hybridize and effectively
behave as a single oscillator in contact with two thermal baths [24, 43]

$\displaystyle\langle P\rangle_{p}$
$\displaystyle=\frac{\kappa_{c}\kappa_{h}\Delta}{\kappa_{c}+\kappa_{h}}(\bar{\text{n}}_{h}-\bar{\text{n}}_{c}),$
(15) $\displaystyle\langle\langle P^{2}\rangle\rangle_{p}$
$\displaystyle=\frac{\kappa_{c}\kappa_{h}\Delta}{(\kappa_{c}+\kappa_{h})^{3}}\big{\\{}(\kappa_{c}+\kappa_{h})^{2}[\bar{\text{n}}_{h}(\bar{\text{n}}_{h}+1)+\bar{\text{n}}_{c}(\bar{\text{n}}_{c}+1)]$
(16) $\displaystyle\hskip
56.9055pt-(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})^{2}(\kappa_{h}^{2}+\kappa_{c}^{2})\big{\\}}.$

To supplement the discussion of noise, we now analyze the Fano factor,
$\mathcal{F}=\langle\langle P^{2}\rangle\rangle/(\langle P\rangle\Delta)$.
This quantity is a measure of bunching and connected to intensity correlations
[44, 45]; e.g., thermal light is (super-)Poissonian, $\mathcal{F}\geq 1$
(bunched) while the flux of single photons is sub-Poissonian (anti-bunched) ,
$\mathcal{F}<1$. For the quantum and particle models, we find
$\mathcal{F}_{q}\geq\mathcal{F}_{p}\geq 1$, while the wave model can attain
sub-Poissonian statistics since $\mathcal{F}_{w}$ can be both smaller and
larger than one, as exemplified in Figs. 2 (b/c). To understand these results,
we introduce $\delta\mathcal{F}_{qw(p)}=\mathcal{F}_{q}-\mathcal{F}_{w(p)}$,
which represent the shaded regions in Fig. 2 (b/c).

For the wave model, Eqs. (10a, 10b) give
$\delta\mathcal{F}_{qw}=(\bar{\text{n}}_{h}+\bar{\text{n}}_{c})/(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})$.
The divergence at equilibrium, $\bar{\text{n}}_{h}=\bar{\text{n}}_{c}$, is due
to the absence of power, while equilibrium fluctuations are still present. We
also note that in the high-temperature limit $\delta\mathcal{F}_{qw}$ does not
vanish, but $\mathcal{F}_{q(w)}$ are dominated by the quadratic terms in
$\bar{\text{n}}_{\alpha}$ such that $\mathcal{F}_{q}\approx\mathcal{F}_{w}$,
as we see in Fig. 2(c). At low values of $g$, the wave model results in anti-
bunching as $\mathcal{F}_{w}$ drops below one when
$\bar{n}_{h}\bar{n}_{c}<\bar{n}_{h}-\bar{n}_{c}$. In this regime, vacuum
fluctuations are crucial to capture the correct statistics, which never show
anti-bunching in this quantum model.

Equations (10a, 10c) readily give
$\delta\mathcal{F}_{qp}=(\mathcal{S}_{p}-\mathcal{S})/(\langle
P\rangle\Delta)\geq 0$, implying reduced bunching for the particle model. This
can be seen in the red shade in Fig. 2 (b/c): For $g/\kappa\approx 2/3$,
$\delta\mathcal{F}_{qp}$ has a maximum, but the particle model still has
$\mathcal{F}_{p}\geq 1$. In contrast, in the limits $g/\kappa\to
0~{}(\infty)$, $\delta\mathcal{F}_{qp}=0$, independently of the temperatures.

Jointly analyzing quantum, wave, and particle models, we conclude that for low
temperatures, $\bar{\text{n}}_{\alpha}\lesssim 1$, and $g/\kappa\approx 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 non-
equilibrium 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).– In stochastic thermodynamics, the
trade-off between power and noise has a well-established bound in terms of the
entropy production rate, $\dot{\sigma}$ [46, 47, 48]. In contrast to fermionic
systems, where the effect of quantum coherence can decrease noise [15], the
so-called TUR, $\langle\langle P^{2}\rangle\rangle/\langle P\rangle^{2}\geq
2/\dot{\sigma}$, 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, $\mathcal{F}\geq
2(\Omega_{c}/T_{c}-\Omega_{h}/T_{h})^{-1}$ [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_{\alpha}$ instead of
$\bar{n}_{\alpha}$, which only coincide for large temperatures. Since
classical Langevin equations obey the TUR [50], a modified bound holds for the
wave model, $\mathcal{F}_{w}\geq
2(\bar{\text{n}}_{h}^{-1}-\bar{\text{n}}_{c}^{-1})^{-1}$ [solid purple in Fig.
2(c)].

Alternative classical models.– The particle and wave models aforementioned are
in principle not unique. In the particle model, for instance, we chose
$\Gamma_{I}=4g^{2}/(\kappa_{h}+\kappa_{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
$\expectationvalue{\xi_{\alpha,\text{in}}^{*}(t^{\prime})\xi_{\beta,\text{in}}(t)}_{w}=(\bar{\text{n}}_{\alpha}+C)\delta_{\alpha\beta}\delta(t^{\prime}-t)$,
with the same $C$ for $\alpha=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 $\bar{\text{n}}_{\alpha}^{2}\to(\bar{\text{n}}_{\alpha}+C)^{2}$, which
should be compared to $\bar{\text{n}}_{\alpha}(\bar{\text{n}}_{\alpha}+1)$ in
the quantum model, see Eq. (10a). For $C=1/2$, the modified wave model
captures the linear terms in $\bar{n}_{\alpha}$ present in Eq. (10a). However,
for any $C\neq 0$, the modified wave model predicts noise at
$\bar{\text{n}}_{h}=\bar{\text{n}}_{c}=0$ and thus cannot correctly reproduce
vacuum fluctuations for arbitrary temperatures.

Conclusions and outlook.–We have shown that the wave-particle 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.

_Acknowledgements._ We acknowledge fruitful discussions with G. Landi and J.
Rabin, and A. Tettamanti for carefully reading the manuscript. This work was
supported by the Swiss National Science Foundation (Eccellenza Professorial
Fellowship PCEFP2_194268).

## References

  * Einstein [1905] A. Einstein, On a heuristic point of view concerning the production and transformation of light, Ann. Phys. 17, 132 (1905).
  * Einstein [1909] A. Einstein, On the present status of the radiation problem, Phys. Zeit. 10, 185 (1909).
  * Weinberg [2015] S. Weinberg,  _Lectures on Quantum Mechanics_, 2nd ed. (Cambridge University Press, 2015).
  * Compton [1923] A. H. Compton, A quantum theory of the scattering of x-rays by light elements, Phys. Rev. 21, 483 (1923).
  * Davisson and Germer [1927] C. Davisson and L. H. Germer, Diffraction of electrons by a crystal of nickel, Phys. Rev. 30, 705 (1927).
  * Oliver  _et al._ [1999] W. D. Oliver, J. Kim, R. C. Liu, and Y. Yamamoto, Hanbury brown and twiss-type experiment with electrons, Science 284, 299 (1999).
  * Wootters and Zurek [1979] W. K. Wootters and W. H. Zurek, Complementarity in the double-slit experiment: Quantum nonseparability and a quantitative statement of bohr’s principle, Phys. Rev. D 19, 473 (1979).
  * Englert [1996] B.-G. Englert, Fringe visibility and which-way information: An inequality, Phys. Rev. Lett. 77, 2154 (1996).
  * Dittel  _et al._ [2021] C. Dittel, G. Dufour, G. Weihs, and A. Buchleitner, Wave-particle duality of many-body quantum states, Phys. Rev. X 11, 031041 (2021).
  * Wang  _et al._ [2019] K. Wang, Q. Xu, S. Zhu, and X. song Ma, Quantum wave–particle superposition in a delayed-choice experiment, Nat. Photonics 13, 872 (2019).
  * Mitchison  _et al._ [2015] M. T. Mitchison, M. P. Woods, J. Prior, and M. Huber, Coherence-assisted single-shot cooling by quantum absorption refrigerators, New J. Phys. 17, 115013 (2015).
  * Nimmrichter  _et al._ [2017] S. Nimmrichter, J. Dai, A. Roulet, and V. Scarani, Quantum and classical dynamics of a three-mode absorption refrigerator, Quantum 1, 37 (2017).
  * Kalaee  _et al._ [2021] A. A. S. Kalaee, A. Wacker, and P. P. Potts, Violating the thermodynamic uncertainty relation in the three-level maser, Phys. Rev. E 104, L012103 (2021).
  * Korzekwa and Lostaglio [2021] K. Korzekwa and M. Lostaglio, Quantum advantage in simulating stochastic processes, Phys. Rev. X 11, 021019 (2021).
  * Prech  _et al._ [2022] K. Prech, P. Johansson, E. Nyholm, G. T. Landi, C. Verdozzi, P. Samuelsson, and P. P. Potts, Entanglement and thermo-kinetic uncertainty relations in coherent mesoscopic transport, arXiv:2212.03835 (2022).
  * Unruh and Zurek [1989] W. G. Unruh and W. H. Zurek, Reduction of a wave packet in quantum brownian motion, Phys. Rev. D 40, 1071 (1989).
  * Zurek  _et al._ [1993] W. H. Zurek, S. Habib, and J. P. Paz, Coherent states via decoherence, Phys. Rev. Lett. 70, 1187 (1993).
  * Zurek [2003] W. H. Zurek, Decoherence, einselection, and the quantum origins of the classical, Rev. Mod. Phys. 75, 715 (2003).
  * Spekkens [2007] R. W. Spekkens, Evidence for the epistemic view of quantum states: A toy theory, Phys. Rev. A 75, 032110 (2007).
  * Spekkens [2016] R. W. Spekkens, Quasi-quantization: Classical statistical theories with an epistemic restriction, in  _Quantum Theory: Informational Foundations and Foils_, edited by G. Chiribella and R. W. Spekkens (Springer Netherlands, Dordrecht, 2016) pp. 83–135.
  * Harrow and Montanaro [2017] A. W. Harrow and A. Montanaro, Quantum computational supremacy, Nature 549, 203 (2017).
  * Mandel [1986] L. Mandel, Non-classical states of the electromagnetic field, Physica Scripta 1986, 34 (1986).
  * Ghosh  _et al._ [2018] A. Ghosh, W. Niedenzu, V. Mukherjee, and G. Kurizki, Thermodynamic principles and implementations of quantum machines, in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, edited by F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso (Springer International Publishing, Cham, 2018) pp. 37–66.
  * Kerremans  _et al._ [2022] T. Kerremans, P. Samuelsson, and P. Potts, Probabilistically violating the first law of thermodynamics in a quantum heat engine, SciPost Physics 12, 168 (2022).
  * Kosloff [1984] R. Kosloff, A quantum mechanical open system as a model of a heat engine, J. Chem. Phys. 80, 1625 (1984).
  * Hofer  _et al._ [2016] P. P. Hofer, J.-R. Souquet, and A. A. Clerk, Quantum heat engine based on photon-assisted cooper pair tunneling, Phys. Rev. B 93, 041418 (2016).
  * Alicki and Kosloff [2018] R. Alicki and R. Kosloff, Introduction to quantum thermodynamics: History and prospects, in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, edited by F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso (Springer International Publishing, Cham, 2018) pp. 1–33.
  * Westig  _et al._ [2017] M. Westig, B. Kubala, O. Parlavecchio, Y. Mukharsky, C. Altimiras, P. Joyez, D. Vion, P. Roche, D. Esteve, M. Hofheinz, M. Trif, P. Simon, J. Ankerhold, and F. Portier, Emission of nonclassical radiation by inelastic cooper pair tunneling, Phys. Rev. Lett. 119, 137001 (2017).
  * Peugeot  _et al._ [2021] A. Peugeot, G. Ménard, S. Dambach, M. Westig, B. Kubala, Y. Mukharsky, C. Altimiras, P. Joyez, D. Vion, P. Roche, D. Esteve, P. Milman, J. Leppäkangas, G. Johansson, M. Hofheinz, J. Ankerhold, and F. Portier, Generating two continuous entangled microwave beams using a dc-biased josephson junction, Phys. Rev. X 11, 031008 (2021).
  * Menard  _et al._ [2022] G. C. Menard, A. Peugeot, C. Padurariu, C. Rolland, B. Kubala, Y. Mukharsky, Z. Iftikhar, C. Altimiras, P. Roche, H. L. Sueur, P. Joyez, D. Vion, D. Esteve, J. Ankerhold, and F. Portier, Emission of photon multiplets by a dc-biased superconducting circuit, Phys. Rev. X 12, 021006 (2022).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Rev. Mod. Phys. 86, 1391 (2014).
  * Monsel  _et al._ [2021] J. Monsel, N. Dashti, S. K. Manjeshwar, J. Eriksson, H. Ernbrink, E. Olsson, E. Torneus, W. Wieczorek, and J. Splettstoesser, Optomechanical cooling with coherent and squeezed light: The thermodynamic cost of opening the heat valve, Phys. Rev. A 103, 063519 (2021).
  * [33] In the supplemental material we provide detailed calculations for the average current and noise in each model and a short note on the connection between power and work variance.
  * Gardiner and Collett [1985] C. W. Gardiner and M. J. Collett, Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation, Phys. Rev. A 31, 3761 (1985).
  * [35] C. W. Gardiner and P. Zoller,  _Quantum Noise A Handbook of Markovian and Non-Markovian Quantum Stochastic Mehtods With Applications to Quantum Optics_.
  * Franken  _et al._ [1961] P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich, Generation of optical harmonics, Phys. Rev. Lett. 7, 118 (1961).
  * Murti and Vijayan [2021] Y. V. G. S. Murti and C. Vijayan, A phenomenological view of nonlinear optics, in  _Physics of Nonlinear Optics_ (Springer International Publishing, Cham, 2021) pp. 9–25.
  * [38] C. W. Gardiner,  _Handbook of Stochastics Methods for Physics, Chemistry and Natural Sciences_.
  * Isserlis [1918] L. Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika 12, 134 (1918).
  * Flindt  _et al._ [2010] C. Flindt, T. Novotny, A. Braggio, and A.-P. Jauho, Counting statistics of transport through coulomb blockade nanostructures: High-order cumulants and non-markovian effects, Phys. Rev. B 82, 155407 (2010).
  * Blanter and Büttiker [2000] Y. Blanter and M. Büttiker, Shot noise in mesoscopic conductors, Phys. Rep. 336, 1 (2000).
  * Clerk  _et al._ [2010] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Introduction to quantum noise, measurement, and amplification, Rev. Mod. Phys. 82, 1155 (2010).
  * Brange  _et al._ [2019] F. Brange, P. Menczel, and C. Flindt, Photon counting statistics of a microwave cavity, Phys. Rev. B 99, 085418 (2019).
  * Basano  _et al._ [2005] L. Basano, P. Ottonello, and B. Torre, Bunching, antibunching, and the poisson limit of bose-einstein processes at low-degeneracy parameters, J. Opt. Soc. Am. B 22, 1314 (2005).
  * Kronwald  _et al._ [2013] A. Kronwald, M. Ludwig, and F. Marquardt, Full photon statistics of a light beam transmitted through an optomechanical system, Phys. Rev. A 87, 013847 (2013).
  * Barato and Seifert [2015] A. C. Barato and U. Seifert, Thermodynamic uncertainty relation for biomolecular processes, Phys. Rev. Lett. 114, 158101 (2015).
  * Gingrich  _et al._ [2016] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, Dissipation bounds all steady-state current fluctuations, Phys. Rev. Lett. 116, 120601 (2016).
  * Horowitz and Gingrich [2020] J. M. Horowitz and T. R. Gingrich, Thermodynamic uncertainty relations constrain non-equilibrium fluctuations, Nat. Phys. 16, 15 (2020).
  * Saryal  _et al._ [2019] S. Saryal, H. M. Friedman, D. Segal, and B. K. Agarwalla, Thermodynamic uncertainty relation in thermal transport, Phys. Rev. E 100, 042101 (2019).
  * Gingrich  _et al._ [2017] T. R. Gingrich, G. M. Rotskoff, and J. M. Horowitz, Inferring dissipation from current fluctuations, J. Phys. A: Math. Theor. 50, 184004 (2017).

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\times\Delta$.

###### Contents

  1. I Quantum Model
    1. I.1 Lindblad Master Equation (LME)
    2. I.2 Quantum Langevin Equations for the quantum model (QLEs)
    3. I.3 Average current and zero-frequency noise
  2. II Wave Model
  3. III Particle Model
    1. III.1 Full counting statistics (FCS)
    2. III.2 Average power
    3. III.3 Noise
    4. III.4 Limitting regimes $g/\kappa\to 0,~{}g/\kappa\to\infty$.
  4. IV Connection between power noise and work variance

##  I Quantum Model

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].

###  I.1 Lindblad Master Equation (LME)

Equation (1) of the main text can be recast in a rotating frame with respect
to $H_{0}=\sum_{\alpha}\Omega_{\alpha}a^{\dagger}_{\alpha}a_{\alpha}$ as

$H=g(a^{\dagger}_{h}a_{c}+a_{c}^{\dagger}a_{h}),$ (S1)

where we used the fact that $\Delta=\Omega_{h}-\Omega_{c}$ and the local
master equation (3) retains the same form. We now work with the adjoint
Liouvillian,

$\bar{\mathcal{L}}A=+i[H,A]+\kappa_{\alpha}(\bar{\text{n}}_{\alpha}+1)\bar{D}[a_{\alpha}]A+\kappa_{\alpha}\bar{\text{n}}_{\alpha}\bar{D}[a_{\alpha}^{\dagger}]A,$
(S2)

with $\bar{D}[L]A=1/2\quantity(L^{\dagger}[A,L]+[L^{\dagger},A]L)$; and
defined through picture equivalence
$\tr(A\mathcal{L}(\rho))=\tr(\bar{\mathcal{L}}(A)\rho)$. From the above
expression, we obtain a closed set of equations of motion for the operators
$(a^{\dagger}_{h}a_{h},a_{c}^{\dagger}a_{c},a^{\dagger}_{h}a_{c},a^{\dagger}_{c}a_{h})=:\bm{\Theta}$,
which can be cast in the form $\bar{\mathcal{L}}\bm{\Theta}=X\bm{\Theta}+Y$:

$\bar{\mathcal{L}}\bm{\Theta}=\underbrace{\begin{pmatrix}-\kappa_{h}&0&-ig&ig\\\
0&-\kappa_{c}&ig&-ig\\\ ig&-ig&-\frac{\kappa_{h}+\kappa_{c}}{2}&0\\\
-ig&ig&0&-\frac{\kappa_{h}+\kappa_{c}}{2}\end{pmatrix}}_{X}\bm{\Theta}+\underbrace{\begin{pmatrix}\kappa_{h}\bar{\text{n}}_{h}\\\
\kappa_{c}\bar{\text{n}}_{c}\\\ 0\\\ 0\end{pmatrix}}_{Y}.$ (S3)

Therefore, the equations of motion for the averages are
$\differential/\differential
t\langle\bm{\Theta}\rangle=X\langle\bm{\Theta}\rangle+Y$.

###  I.2 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
$\differential/\differential
t\langle\bm{\Theta}\rangle=X\langle\bm{\Theta}\rangle+Y$. We start by
recasting the QLE’s of the main text (4) in a rotating frame generated by
$\hat{r}=\sum_{\alpha}\Omega_{\alpha}\quantity[a^{\dagger}_{\alpha}a_{\alpha}+\int_{0}^{\infty}\differential\omega
b_{\alpha}^{\dagger}(\omega)b_{\alpha}(\omega)],$ where $b_{\alpha}(\omega)$
are bosonic field operators of the baths. In this frame, we can eliminate the
time-dependence and the QLEs write,

$\displaystyle\dot{a}_{h}$
$\displaystyle=-\frac{\kappa_{h}}{2}a_{h}-iga_{c}-\sqrt{\kappa_{h}}b_{h,\text{in}},$
(S4a) $\displaystyle\dot{a}_{c}$
$\displaystyle=-\frac{\kappa_{c}}{2}a_{h}-iga_{h}-\sqrt{\kappa_{c}}b_{c,\text{in}}.$
(S4b)

With the above equation, we can again compute the equations of motion for
$\bm{\Theta}=(a^{\dagger}_{h}a_{h},a_{c}^{\dagger}a_{c},a^{\dagger}_{h}a_{c},a^{\dagger}_{c}a_{h})$,
using Leibniz’s rule. After some algebra, we obtain $\bm{\Theta}$

$\displaystyle\derivative{\Theta}{t}=\underbrace{\begin{pmatrix}-\kappa_{h}&0&-ig&ig\\\
0&-\kappa_{c}&ig&-ig\\\ ig&-ig&-\frac{\kappa_{h}+\kappa_{c}}{2}&0\\\
-ig&ig&0&-\frac{\kappa_{h}+\kappa_{c}}{2}\end{pmatrix}}_{X}\Theta-\underbrace{\begin{pmatrix}\sqrt{\kappa_{h}}a^{\dagger}_{h}b_{h,\text{in}}\\\
\sqrt{\kappa_{c}}a^{\dagger}_{c}b_{c,\text{in}}\\\
\sqrt{\kappa_{h}}b_{h,\text{in}}^{\dagger}a_{c}+\sqrt{\kappa_{c}}a_{h}^{\dagger}b_{c,\text{in}}\\\
\sqrt{\kappa_{c}}a_{h}b_{c,\text{in}}^{\dagger}+\sqrt{\kappa_{c}}a_{c}^{\dagger}b_{h,\text{in}}\end{pmatrix}}_{y}.$
(S5)

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, $\differential/\differential
t\expectationvalue{\bm{\Theta}}$; note that the term proportional to
$\expectationvalue{\bm{\Theta}}$ is already $X$, the homogeneous term in Eq.
(S3), and it remains to be shown that $\expectationvalue{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_{\alpha}^{(\dagger)}$ 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_{\alpha}(t)\xrightarrow{\mathcal{F}_{\nu}}a_{\alpha}(\nu)$

$\displaystyle a_{h}(\nu)$
$\displaystyle=\frac{-iga_{c}(\nu)-\sqrt{\kappa_{h}}b_{h,\text{in}}(\nu)}{i\nu+\frac{\kappa_{h}}{2}},$
(S6) $\displaystyle a_{c}(\nu)$
$\displaystyle=\frac{-iga_{h}(\nu)-\sqrt{\kappa_{c}}b_{c,\text{in}}(\nu)}{i\nu+\frac{\kappa_{c}}{2}}.$
(S7)

Solving for the system modes we obtain the following relations,

$\displaystyle a_{h}(\nu)$
$\displaystyle=\frac{-2\sqrt{\kappa_{h}}b_{h,\text{in}}(\nu)(\kappa_{c}+2i\nu)+4ig\sqrt{\kappa_{c}}b_{c,\text{in}}(\nu)}{4g^{2}+(\kappa_{h}+2i\nu)(\kappa_{c}+2i\nu)},$
(S8a) $\displaystyle a_{c}(\nu)$
$\displaystyle=\frac{-2\sqrt{\kappa_{c}}b_{c,\text{in}}(\nu)(\kappa_{h}+2i\nu)+4ig\sqrt{\kappa_{h}}b_{h,\text{in}}(\nu)}{4g^{2}+(\kappa_{c}+2i\nu)(\kappa_{h}+2i\nu)}.$
(S8b)

We now consider, for concreteness, the first entry of $\expectationvalue{y}$,
with each term written as a Fourier transform,

$\displaystyle\sqrt{\kappa_{h}}\expectationvalue{a_{h}^{\dagger}(t)b_{h,\text{in}}}$
$\displaystyle=\sqrt{\kappa_{h}}\int\frac{\differential\nu\differential\nu^{\prime}}{2\pi}e^{i(\nu+\nu^{\prime})t}\frac{-2\sqrt{\kappa_{h}}\expectationvalue{b_{h,\text{in}}^{\dagger}(\nu)b_{h,\text{in}}(\nu^{\prime})}(\kappa_{c}+2i\nu)-4ig\sqrt{\kappa_{c}}\expectationvalue{b_{c,\text{in}}^{\dagger}(\nu)b_{h,\text{in}}(\nu^{\prime})}}{4g^{2}+(\kappa_{h}+2i\nu)(\kappa_{c}+2i\nu)}$
(S9)
$\displaystyle=-\frac{\kappa_{h}\bar{\text{n}}_{h}}{\pi}\int_{-\infty}^{+\infty}\differential\nu~{}\frac{(\kappa_{c}+2i\nu)}{4g^{2}+(\kappa_{h}+2i\nu)(\kappa_{c}+2i\nu)},$
(S10)

where, in the last passage, we considered the quantum white noise relations in
Fourier space,
$\expectationvalue{b_{\alpha,\text{in}}^{\dagger}(\nu)b_{\beta,\text{in}}(-\nu^{\prime})}=\bar{\text{n}}_{\alpha}\delta_{\alpha\beta}\delta(\nu-\nu^{\prime}).$
The last integral can be evaluated by finding the poles of the integrand,

$\displaystyle\nu_{1}=\frac{1}{4}\quantity(i(\kappa_{h}+\kappa_{c})-\sqrt{16g^{2}-(\kappa_{h}-\kappa_{c})^{2}}),$
(S11a)
$\displaystyle\nu_{2}=\frac{1}{4}\quantity(i(\kappa_{h}+\kappa_{c})+\sqrt{16g^{2}-(\kappa_{h}-\kappa_{c})^{2}}).$
(S11b)

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,

$\displaystyle\int_{-\infty}^{\infty}\differential\nu\frac{\kappa_{c}+2i\nu}{4g^{2}+(\kappa_{h}-2i\nu)(\kappa_{c}-2i\nu)}=:\int_{-\infty}^{\infty}\differential\nu
f(\nu)=\int_{C}\differential zf(z)=2\pi
i\sum_{i=1,2}\text{Res}(f(z),\nu_{i})=\pi.$ (S12)

We are then left with,

$\displaystyle\sqrt{\kappa_{h}}\expectationvalue{a^{\dagger}_{h}(t)b_{h,\text{in}}(t)}=-\kappa_{h}\bar{\text{n}}_{h}.$
(S13)

Similarly,

$\displaystyle\sqrt{\kappa_{c}}\expectationvalue{a^{\dagger}_{c}(t)b_{c,\text{in}}(t)}=-\kappa_{c}\bar{\text{n}}_{c}.$
(S14)

The remaining entries of $\langle y\rangle$ are integrals of the form

$\displaystyle\sqrt{\kappa_{c}}\expectationvalue{a_{h}^{\dagger}b_{c,\text{in}}}=4ig\bar{\text{n}}_{h}\kappa_{c}\int_{-\infty}^{\infty}\differential\nu\frac{1}{4g^{2}+(\kappa_{c}+2i\nu)(\kappa_{h}+2i\nu)}=0,$
(S15)

where we again use quantum white noise relations and the same complex-contour
integration. Thus, we have shown that $\langle y\rangle=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.

###  I.3 Average current and zero-frequency noise

We start by recasting the equations of motion in a new basis vector of
operators $\bm{\Theta}\to\bm{\sigma}=(I,H,N_{h},N_{c})$ where
$I=ig(a^{\dagger}_{h}a_{c}-a_{c}^{\dagger}a_{h})$ is the current operator in
the new frame, $H$ is Eq. (S1) and
$N_{\alpha}=a^{\dagger}_{\alpha}a_{\alpha}$. The equations of motion are now
in the form $\differential/\differential t\bm{\sigma}=G\bm{\sigma}+F$ with

$G=\begin{pmatrix}-\frac{\kappa_{h}+\kappa_{c}}{2}&0&2g^{2}&-2g^{2}\\\
0&-\frac{\kappa_{h}+\kappa_{c}}{2}&0&0\\\ -1&0&-\kappa_{h}&0\\\
1&0&0&-\kappa_{c}\end{pmatrix},$ (S16)

and $F=(0,0,\bar{\text{n}}_{h}\kappa_{h},\bar{\text{n}}_{c}\kappa_{c})$.

We readily obtain steady-state averages by solving
$\expectationvalue{\bm{\sigma}}=-G^{-1}F$. In particular, we obtain the
average current,

$\expectationvalue{I}_{q}=\frac{4g^{2}\kappa_{h}\kappa_{c}(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})}{(4g^{2}+\kappa^{2}_{h}\kappa_{c})(\kappa_{h}+\kappa_{c})}.$
(S17)

Setting $\kappa_{h}=\kappa_{c}=\kappa$ 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
$\delta\bm{\sigma}=\bm{\sigma}-\expectationvalue{\bm{\sigma}}$ which evolves
according to the homogeneous equation,

$\bar{\mathcal{L}}\delta\bm{\sigma}=G\delta\bm{\sigma},$ (S18)

and consider the correlator $\expectationvalue{\delta\bm{\sigma}\delta I}$,
whose first entry is $\expectationvalue{\delta I\delta I}$. The above
homogeneous equation satisfies the requirements of the quantum regression
theorem, which implies that

$\expectationvalue{\delta\bm{\sigma}(t)\delta
I(0)}=\expectationvalue{e^{tG}\quantity[\delta\bm{\sigma}(0)]\delta I(0)}.$
(S19)

Integrating the above as in Eq.(2), we have

$\displaystyle\expectationvalue{\expectationvalue{\bm{\sigma}I}}=2\real\int_{0}^{\infty}\differential
t\expectationvalue{\delta\bm{\sigma}(t)\delta
I(0)}=-2G^{-1}\expectationvalue{\delta\bm{\sigma}(0)\delta I(0)},$ (S20)

whose first entry is the desired cumulant. Given the propagator in Eq. (S16),
we just need to compute the initial conditions,

$\expectationvalue{\delta\bm{\sigma}(0)\delta I(0)}=\begin{pmatrix}\langle
I^{2}\rangle-\langle I\rangle^{2}\\\ \langle HI\rangle-\langle I\rangle\langle
H\rangle\\\ \langle N_{h}I\rangle-\langle I\rangle\langle N_{h}\rangle\\\
\langle N_{c}I\rangle-\langle I\rangle\langle N_{c}\rangle\end{pmatrix}.$
(S21)

For instance, to compute $\expectationvalue{I^{2}}$ we need four-point
functions,

$\displaystyle\expectationvalue{I^{2}}=-g^{2}\expectationvalue{a^{\dagger}_{h}a_{c}a^{\dagger}_{h}a_{c}+a^{\dagger}_{c}a_{h}a^{\dagger}_{c}a_{h}-a_{h}^{\dagger}a_{c}a^{\dagger}_{c}a_{h}-a^{\dagger}_{c}a_{h}a^{\dagger}_{h}a_{c}},$
(S22)

due to Gaussianity, we can compute this correlator through Wick contractions,
e.g.

$\displaystyle\expectationvalue{a^{\dagger}_{\mu}a_{\nu}a^{\dagger}_{\gamma}a_{\sigma}}=\expectationvalue{a^{\dagger}_{\mu}a_{\nu}}\expectationvalue{a^{\dagger}_{\gamma}a_{\sigma}}+\expectationvalue{a^{\dagger}_{\mu}a_{\sigma}}\expectationvalue{a^{\dagger}_{\gamma}a_{\nu}},$
(S23)

where we also used
$\expectationvalue{a_{\mu}^{(\dagger)}a_{\nu}^{(\dagger)}}=0$. To compute the
desired initial conditions, all the following contractions are required:

$\displaystyle\langle a^{\dagger}_{h}a_{c}a^{\dagger}_{h}a_{c}\rangle$
$\displaystyle=2\langle a^{\dagger}_{h}a_{c}\rangle^{2},$ (S24a)
$\displaystyle\langle a^{\dagger}_{c}a_{h}a^{\dagger}_{c}a_{h}\rangle$
$\displaystyle=2\langle a^{\dagger}_{c}a_{h}\rangle^{2},$ (S24b)
$\displaystyle\langle a^{\dagger}_{h}a_{c}a^{\dagger}_{c}a_{h}\rangle$
$\displaystyle=\quantity|\langle a^{\dagger}_{h}a_{c}\rangle|^{2}+\langle
a^{\dagger}_{h}a_{h}\rangle\langle a^{\dagger}_{c}a_{c}+1\rangle,$ (S24c)
$\displaystyle\langle a^{\dagger}_{c}a_{h}a^{\dagger}_{h}a_{c}\rangle$
$\displaystyle=\quantity|\langle a^{\dagger}_{c}a_{h}\rangle|^{2}+\langle
a^{\dagger}_{c}a_{c}\rangle\langle a^{\dagger}_{h}a_{h}+1\rangle,$ (S24d)
$\displaystyle\langle a^{\dagger}_{h}a_{h}a^{\dagger}_{h}a_{c}\rangle$
$\displaystyle=\langle
a^{\dagger}_{h}a_{c}\rangle\expectationvalue{2a^{\dagger}_{h}a_{h}+1},$ (S24e)
$\displaystyle\langle a^{\dagger}_{c}a_{c}a^{\dagger}_{c}a_{h}\rangle$
$\displaystyle=\langle a^{\dagger}_{c}a_{h}\rangle\langle
2a^{\dagger}_{c}a_{c}+1\rangle,$ (S24f) $\displaystyle\langle
a^{\dagger}_{h}a_{h}a^{\dagger}_{c}a_{h}\rangle$ $\displaystyle=2\langle
a^{\dagger}_{h}a_{h}\rangle\langle a^{\dagger}_{c}a_{h}\rangle,$ (S24g)
$\displaystyle\langle a^{\dagger}_{c}a_{c}a^{\dagger}_{h}a_{c}\rangle$
$\displaystyle=2\langle a^{\dagger}_{c}a_{c}\rangle\langle
a^{\dagger}_{h}a_{c}\rangle.$ (S24h)

At the instance of Eq. (S22), we get by applying the above,

$\displaystyle\langle I^{2}\rangle=-g^{2}\quantity[2\langle
a^{\dagger}_{h}a_{c}\rangle^{2}+2\langle
a^{\dagger}_{c}a_{h}\rangle^{2}-2\quantity|\langle
a^{\dagger}_{h}a_{c}\rangle|^{2}-\langle a^{\dagger}_{h}a_{h}\rangle\langle
a^{\dagger}_{c}a_{c}+1\rangle-\langle a^{\dagger}_{c}a_{c}\rangle\langle
a^{\dagger}_{h}a_{h}+1\rangle].$ (S25)

The steady-state solutions for the covariance vector
$\expectationvalue{\bm{\Theta}}$ fully determine $\langle I^{2}\rangle$.
Similarly, for the other entries of the initial conditions,

$\displaystyle\langle HI\rangle$ $\displaystyle=ig^{2}\quantity[2\langle
a^{\dagger}_{h}a_{c}\rangle^{2}-2\langle
a^{\dagger}_{c}a_{h}\rangle^{2}+\langle a^{\dagger}_{c}a_{c}\rangle\langle
a^{\dagger}_{h}a_{h}+1\rangle-\langle a^{\dagger}_{h}a_{h}\rangle\langle
a^{\dagger}_{c}a_{c}+1\rangle],$ (S26) $\displaystyle\langle N_{h}I\rangle$
$\displaystyle=ig\quantity[\langle a^{\dagger}_{h}a_{c}\rangle\langle
2a^{\dagger}_{h}a_{h}+1\rangle-2\langle a^{\dagger}_{c}a_{h}\rangle\langle
a^{\dagger}_{h}a_{h}\rangle],$ (S27) $\displaystyle\langle N_{c}I\rangle$
$\displaystyle=-ig\quantity[\langle a^{\dagger}_{c}a_{h}\rangle\langle
2a^{\dagger}_{c}a_{c}+1\rangle-2\langle a^{\dagger}_{h}a_{c}\rangle\langle
a^{\dagger}_{c}a_{c}\rangle].$ (S28)

By inverting $G$ and applying to $\langle\delta\bm{\sigma}\delta I\rangle$ we
obtain a vector, whose first entry is the current noise,

$\expectationvalue{\expectationvalue{I^{2}}}_{q}=\frac{2g^{2}\kappa_{h}\kappa_{c}\quantity((\bar{\text{n}}_{h}+\bar{\text{n}}_{c})^{2}+2(\bar{\text{n}}_{h}+\bar{\text{n}}_{c}))}{(4g^{2}+\kappa_{h}\kappa_{c})(\kappa_{h}+\kappa_{c})}+\frac{2g^{2}\kappa_{h}\kappa_{c}(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})^{2}\quantity(16g^{4}(\kappa_{h}-\kappa_{c})^{2}+\kappa_{h}^{2}\kappa_{c}^{2}(\kappa_{h}+\kappa_{c})^{2}-8g^{2}\kappa_{h}\kappa_{c}(\kappa_{h}^{2}+\kappa_{c}^{2}+4\kappa_{h}\kappa_{c}))}{(4g^{2}+\kappa_{h}\kappa_{c})^{3}(\kappa_{h}+\kappa_{c})^{3}}.$
(S29)

Finally, considering $\langle\langle P^{2}\rangle\rangle=\langle\langle
I^{2}\rangle\rangle\Delta^{2}$, we write the above in terms of equilibrium and
shot noise,

$\displaystyle\expectationvalue{\expectationvalue{P^{2}}}_{q}$
$\displaystyle=\mathcal{E}[\bar{\text{n}}_{h}(\bar{\text{n}}_{h}+1)+\bar{\text{n}}_{c}(\bar{\text{n}}_{c}+1)]-\mathcal{S}(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})^{2},$
(S30) $\displaystyle\mathcal{E}$
$\displaystyle=4g^{2}\kappa_{h}\kappa_{c}\Delta^{2}\chi,$ (S31)
$\displaystyle\mathcal{S}$
$\displaystyle=4g^{2}\kappa_{h}\kappa_{c}(-8g^{2}\kappa_{h}^{2}\kappa_{c}^{2}+\kappa_{h}^{2}\kappa_{c}^{2}(\kappa_{h}+\kappa_{c})^{2}+16g^{4}(\kappa_{h}^{2}+\kappa_{c}^{2}))\Delta^{2}\chi^{3},$
(S32)

where $\chi=[(\kappa_{h}+\kappa_{c})(4g^{2}+\kappa_{h}\kappa_{c})]^{-1}$.
Setting $\kappa_{h}=\kappa_{c}=\kappa$ we obtain Eq. (10a) of the main text.

##  II 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))

$\displaystyle\dot{A}_{h}=-\frac{\kappa_{h}}{2}A_{h}-igA_{c}-\sqrt{\kappa_{h}}\xi_{h,\text{in}},$
(S33a)
$\displaystyle\dot{A}_{c}=-\frac{\kappa_{c}}{2}A_{h}-igA_{h}-\sqrt{\kappa_{c}}\xi_{c,\text{in}}.$
(S33b)

From the above we obtain the equations of motion for
$\Theta=\quantity(A^{*}_{h}A_{h},A^{*}_{c}A_{c},A^{*}_{h}A_{c},A^{*}_{c}A_{h})$,
and, in particular, for the averages

$\displaystyle\derivative{\Theta}{t}=\begin{pmatrix}-\kappa_{h}&0&-ig&ig\\\
0&-\kappa_{c}&ig&-ig\\\ ig&-ig&-\frac{\kappa_{h}+\kappa_{c}}{2}&0\\\
-ig&ig&0&-\frac{\kappa_{h}+\kappa_{c}}{2}\end{pmatrix}\Theta-\begin{pmatrix}\sqrt{\kappa_{h}}\expectationvalue{A^{*}_{h}\xi_{h,\text{in}}}\\\
\sqrt{\kappa_{c}}\expectationvalue{A^{*}_{c}\xi_{c,\text{in}}}\\\
\sqrt{\kappa_{h}}\expectationvalue{A_{c}\xi_{h,\text{in}}^{*}}+\sqrt{\kappa_{c}}\expectationvalue{A_{h}^{*}\xi_{c,\text{in}}}\\\
\sqrt{\kappa_{c}}\expectationvalue{A_{h}\xi_{c,\text{in}}^{*}}+\sqrt{\kappa_{c}}\expectationvalue{A_{c}^{*}\xi_{h,\text{in}}}\end{pmatrix}.$
(S34)

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 $\differential/\differential
t\langle\bm{\Theta}\rangle=X\langle\bm{\Theta}\rangle+Y^{\prime}$

$\displaystyle\derivative{\expectationvalue{\Theta}}{t}=\underbrace{\begin{pmatrix}-\kappa_{h}&0&-ig&ig\\\
0&-\kappa_{c}&ig&-ig\\\ ig&-ig&-\frac{\kappa_{h}+\kappa_{c}}{2}&0\\\
-ig&ig&0&-\frac{\kappa_{h}+\kappa_{c}}{2}\end{pmatrix}}_{X}\expectationvalue{\Theta}+\underbrace{\begin{pmatrix}\kappa_{h}\bar{\Phi}_{h}\\\
\kappa_{c}\bar{\Phi}_{c}\\\ 0\\\ 0\end{pmatrix}}_{Y^{\prime}},$ (S35)

where we have considered a generalized form of classical white noise
$\langle\xi_{\alpha,\text{in}}^{*}(t^{\prime})\xi_{\beta,\text{in}}(t)\rangle=\bar{\Phi}_{\alpha}\delta_{\alpha\beta}\delta(t^{\prime}-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 $\bm{\sigma}=(I,\mathcal{H},N_{h},N_{c})$, where
$\mathcal{H}=g(A^{*}_{h}A_{c}+A^{*}_{c}A_{c})$, is the classical Hamiltonian,
$N_{\alpha}=A^{*}_{\alpha}A_{\alpha}$ and the wave current,

$\displaystyle I=ig(A^{*}_{h}A_{c}-A^{*}_{c}A_{c}).$ (S36)

With this, $\differential/\differential
t\langle\bm{\sigma}\rangle=G\langle\bm{\sigma}\rangle+F^{\prime}$, where we
emphasize that coefficient matrix $G$ is the same as the quantum model’s $G$,
Eq.(S16) and
$F^{\prime}=(0,0,\bar{\Phi}_{h}\kappa_{h},\bar{\Phi}_{c}\kappa_{c})$. We
readily obtain the wave current, from
$\langle\bm{\sigma}\rangle=-GF^{\prime}$,

$\expectationvalue{I}_{w}=\frac{4g^{2}\kappa_{h}\kappa_{c}(\bar{\Phi}_{h}-\bar{\Phi}_{c})}{(4g^{2}+\kappa^{2}_{h}\kappa_{c})(\kappa_{h}+\kappa_{c})},$
(S37)

which reduces to the current in the quantum model for
$\bar{\Phi}_{\alpha}=\bar{\text{n}}_{\alpha}+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,

$\displaystyle\langle\langle\bm{\sigma}I\rangle\rangle=2\int_{0}^{\infty}\differential
te^{-tG}\langle\delta\bm{\sigma}(t)\delta
I(0)\rangle=-2G^{-1}\langle\delta\bm{\sigma}(0)\delta I(0)\rangle,$ (S38)

where $G$ is again the coefficient matrix of the equations of motion for
$\bm{\sigma}$. It remains to compute the initial conditions,

$\displaystyle\expectationvalue{\delta\bm{\sigma}\delta
I}=\begin{pmatrix}\expectationvalue{I^{2}}-\expectationvalue{I}^{2}\\\
\expectationvalue{\mathcal{H}I}-\expectationvalue{I}\expectationvalue{\mathcal{H}}\\\
\expectationvalue{N_{h}I}-\expectationvalue{I}\expectationvalue{N_{h}}\\\
\expectationvalue{N_{c}I}-\expectationvalue{I}\expectationvalue{N_{c}}\end{pmatrix}.$
(S39)

Inspecting $\expectationvalue{I^{2}}$ we realize that we need four-point
functions,

$\displaystyle\expectationvalue{I^{2}}=-g^{2}\expectationvalue{A^{*}_{h}A_{c}A^{*}_{h}A_{c}+A^{*}_{c}A_{h}A^{*}_{c}A_{h}-A_{h}^{*}A_{c}A^{*}_{c}A_{h}-A^{*}_{c}A_{h}A^{*}_{h}A_{c}},$
(S40)

due to Gaussianity, we can compute these correlators through Isserlis
contractions, e.g.

$\displaystyle\expectationvalue{A^{*}_{\mu}A_{\nu}A^{*}_{\gamma}A_{\sigma}}=\expectationvalue{A^{*}_{\mu}A_{\nu}}\expectationvalue{A^{*}_{\gamma}A_{\sigma}}+\expectationvalue{A^{*}_{\mu}A_{\sigma}}\expectationvalue{A^{*}_{\gamma}A_{\nu}}.$
(S41)

In general, all the following contractions will be needed

$\displaystyle\langle A^{*}_{h}A_{c}A^{*}_{h}A_{c}\rangle$
$\displaystyle=2\langle A^{*}_{h}A_{c}\rangle^{2},$ (S42a)
$\displaystyle\langle A^{*}_{c}A_{h}A^{*}_{c}A_{h}\rangle$
$\displaystyle=2\langle A^{*}_{c}A_{h}\rangle^{2},$ (S42b)
$\displaystyle\langle A^{*}_{h}A_{c}A^{*}_{c}A_{h}\rangle$
$\displaystyle=\quantity|\langle A^{*}_{h}A_{c}\rangle|^{2}+\langle
A^{*}_{h}A_{c}\rangle\langle A^{*}_{c}A_{c}\rangle,$ (S42c)
$\displaystyle\langle A^{*}_{c}A_{h}A^{*}_{h}A_{c}\rangle$
$\displaystyle=\quantity|\langle A^{*}_{c}A_{h}\rangle|^{2}+\langle
A^{*}_{c}A_{h}\rangle\langle A^{*}_{h}A_{h}\rangle,$ (S42d)
$\displaystyle\langle A^{*}_{h}A_{h}A^{*}_{h}A_{c}\rangle$
$\displaystyle=2\langle A^{*}_{h}A_{c}\rangle\langle A^{*}_{h}A_{h}\rangle,$
(S42e) $\displaystyle\langle A^{*}_{c}A_{c}A^{*}_{c}A_{h}\rangle$
$\displaystyle=2\langle A^{*}_{c}A_{h}\rangle\langle A^{*}_{c}A_{c}\rangle,$
(S42f) $\displaystyle\langle A^{*}_{h}A_{h}A^{*}_{c}A_{h}\rangle$
$\displaystyle=2\langle A^{*}_{h}A_{h}\rangle\langle A^{*}_{c}A_{h}\rangle,$
(S42g) $\displaystyle\langle A^{*}_{c}A_{c}A^{*}_{h}A_{c}\rangle$
$\displaystyle=2\langle A^{*}_{c}A_{c}\rangle\langle A^{*}_{h}A_{c}\rangle.$
(S42h)

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_{\alpha}$
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

$\expectationvalue{\expectationvalue{I^{2}}}_{w}=\frac{2g^{2}\kappa_{h}\kappa_{c}\quantity(\bar{\Phi}_{h}+\bar{\Phi}_{c})^{2}}{(4g^{2}+\kappa_{h}\kappa_{c})(\kappa_{h}+\kappa_{c})}+\frac{2g^{2}\kappa_{h}\kappa_{c}(\bar{\Phi}_{h}-\bar{\Phi}_{c})^{2}\quantity(16g^{2}(\kappa_{h}-\kappa_{c})^{2}+\kappa_{h}^{2}\kappa_{c}^{2}(\kappa_{h}+\kappa_{c})^{2}-8g^{2}\kappa_{h}\kappa_{c}(\kappa_{h}^{2}+\kappa_{c}^{2}+4\kappa_{h}\kappa_{c}))}{(4g^{2}+\kappa_{h}\kappa_{c})^{3}(\kappa_{h}+\kappa_{c})^{3}}.$
(S43)

Finally, using $\langle\langle P^{2}\rangle\rangle=\Delta^{2}\langle\langle
I\rangle\rangle$ and writing the above in terms of equilibrium and shot noise,

$\displaystyle\expectationvalue{\expectationvalue{P^{2}}}_{w}$
$\displaystyle=\mathcal{E}[\bar{\Phi}_{h}^{2}+\bar{\Phi}_{c}^{2}]-\mathcal{S}(\bar{\Phi}_{h}-\bar{\Phi}_{c})^{2},$
(S44) $\displaystyle\mathcal{E}$
$\displaystyle=4g^{2}\kappa_{h}\kappa_{c}\Delta^{2}\chi,$ (S45)
$\displaystyle\mathcal{S}$
$\displaystyle=4g^{2}\kappa_{h}\kappa_{c}(-8g^{2}\kappa_{h}^{2}\kappa_{c}^{2}+\kappa_{h}^{2}\kappa_{c}^{2}(\kappa_{h}+\kappa_{c})^{2}+16g^{4}(\kappa_{h}^{2}+\kappa_{c}^{2}))\Delta^{2}\chi^{3},$
(S46)

where $\chi=[(\kappa_{h}+\kappa_{c})(4g^{2}+\kappa_{h}\kappa_{c})]^{-1}$.
Setting $\kappa_{h}=\kappa_{c}=\kappa$,
$\bar{\Phi}_{\alpha}=\bar{\text{n}}_{\alpha}$ and using $\langle\langle
P^{2}\rangle\rangle=\Delta^{2}\langle\langle I\rangle\rangle$ we obtain Eq.
(10b) of the main text. The discussion for alternative classical models can be
supplemented by taking $\bar{\Phi}_{\alpha}=\bar{\text{n}}_{\alpha}+C$.

##  III Particle Model

We start by recalling Eq. (8) of the main text,

$\displaystyle\dot{p}_{n_{h},n_{c}}$
$\displaystyle=\kappa_{h}(\bar{\text{n}}_{h}+1)(n_{h}+1)p_{n_{h}+1,n_{c}}+\kappa_{h}\bar{\text{n}}_{h}n_{h}p_{n_{h}-1,n_{c}}+\kappa_{c}(\bar{\text{n}}_{c}+1)(n_{c}+1)p_{n_{h},n_{c}+1}+\kappa_{c}\bar{\text{n}}_{c}n_{c}p_{n_{h},n_{c}-1}$
$\displaystyle+\Gamma_{I}(n_{h}+1)n_{c}p_{n_{h}+1,n_{c}-1}+\Gamma_{I}(n_{c}+1)n_{h}p_{n_{h}-1,n_{c}+1}-\Gamma^{0}_{n_{h},n_{c}}p_{n_{h},n_{c}},$

with

$\Gamma_{n_{h},n_{c}}^{0}=\sum_{\alpha}\quantity(\Gamma^{\downarrow}_{\alpha}n_{\alpha}+\Gamma^{\uparrow}_{\alpha}(n_{\alpha}+1))+\Gamma_{I}\quantity((n_{h}+1)n_{c}+(n_{c}+1)n_{h}).$
(S47)

For convenience, we can rewrite the above in matrix form; we start by
introducing the vector
$\left|p\right)=\sum_{n_{h},n_{c}}p_{n_{h},n_{c}}\left|n_{h},n_{c}\right)$ and
the stochastic operator $\mathscr{L}$, s.t.

$\displaystyle\derivative{t}\left|p\right)$
$\displaystyle=\mathscr{L}\left|p\right),$ (S48) $\displaystyle\mathscr{L}$
$\displaystyle=\sum_{\alpha}\Gamma^{\uparrow}_{\alpha}\mathscr{J}_{\alpha}^{\uparrow}+\Gamma^{\downarrow}_{\alpha}\mathscr{J}_{\alpha}^{\downarrow}-\Gamma^{\uparrow}_{\alpha}\mathscr{N}_{\alpha}-\Gamma^{\downarrow}_{\alpha}(\mathscr{N}_{\alpha}+1)$
(S49)
$\displaystyle~{}~{}~{}~{}~{}+\Gamma_{I}\quantity(\mathscr{V}^{+}+\mathscr{V}^{-}-(\mathscr{N}_{h}+1)\mathscr{N}_{c}-(\mathscr{N}_{c}+1)\mathscr{N}_{h}).$

Above, we introduced the number operators

$\displaystyle\mathscr{N}_{h}$
$\displaystyle=\sum_{n_{h}}n_{h}\left|n_{h}\right)\left(n_{h}\right|,$
$\displaystyle\mathscr{N}_{c}$
$\displaystyle=\sum_{n_{c}}n_{c}\left|n_{c}\right)\left(n_{c}\right|,$ (S50)

and the jump operators,

$\displaystyle\mathscr{J}_{\alpha}^{\uparrow}$
$\displaystyle=\sum_{n_{\alpha}}n_{\alpha}\left|n_{\alpha}-1\right)\left(n_{\alpha}\right|,$
$\displaystyle\mathscr{J}_{\alpha}^{\downarrow}$
$\displaystyle=\sum_{n_{\alpha}}(n_{\alpha}+1)\left|n_{\alpha}+1\right)\left(n_{\alpha}\right|,$
(S51) $\displaystyle\mathscr{V}^{+}$
$\displaystyle=\sum_{n_{h},n_{c}}n(m+1)\left|n_{h}-1,n_{c}+1\right)\left(n_{h},n_{c}\right|,$
$\displaystyle\mathscr{V}^{-}$
$\displaystyle=\sum_{n_{h},n_{c}}m(n+1)\left|n_{h}+1,n_{c}-1\right)\left(n_{h},n_{c}\right|.$
(S52)

We note that in this vector space formalism, $\quantity(1\big{|}p)=1$ where
$\left(1\right|$ is a vector whose entries are all 1, represents the
normalization of a probability distribution and $\langle
X\rangle=\left(1\right|X\left|p\right)$ gives the expected value relative to
the probability $p$.

###  III.1 Full counting statistics (FCS)

Equation (S48) is our starting point to establish current fluctuations through
full counting statistics [40]. We want to monitor the intra-system particle
flux; in our particle model, these jumps are captured by the operators
$\mathscr{V}^{\pm}$ and, thus, we dress them with the counting field, $\chi$,
through $\mathscr{V}^{\pm}\to\mathscr{V}^{\pm}e^{\pm i\chi}$. We then obtain a
generalized master equation for the counting-field dressed probability
distribution

$\displaystyle\derivative{t}\left|p_{\chi}\right)=\mathscr{L}(\chi)\left|p_{\chi}\right)=(\mathscr{L}+\Delta\mathscr{L}(\chi))\left|p_{\chi}\right),$
(S53)

where, we have that

$\displaystyle\Delta\mathscr{L}(\chi)=\Gamma_{I}(e^{i\chi}-1)\mathscr{V}^{+}+\Gamma_{I}(e^{-i\chi}-1)\mathscr{V}^{-}.$
(S54)

The cumulant generating function of jump counts is connected to the
normalization of $\left|p_{\chi}\right)$ through
$S(\chi)=\ln\quantity(1\big{|}p_{\chi})$, where $\left(1\right|$ is a vector
whose entries are all 1. $S(\chi)$ is directly connected to the cumulants of
particle current,

$\langle\langle
I^{n}\rangle\rangle=\derivative{t}\frac{\differential^{n}}{\differential\chi^{n}}S(\chi)\Big{|}_{\chi=0}.$
(S55)

The first and second cumulants of the current are given by [40]

$\displaystyle\langle\langle I\rangle\rangle$
$\displaystyle=\langle\mathscr{W}_{1}\rangle,$ (S56)
$\displaystyle\langle\langle I^{2}\rangle\rangle$
$\displaystyle=\langle\mathscr{W}_{2}\rangle-2\langle\mathscr{W}_{1}\mathscr{L}^{D}\mathscr{W}_{1}\rangle,$
(S57)

where
$\mathscr{W}_{n}=\differential^{n}/\differential(i\chi)^{n}\Delta\mathscr{L}(\chi)|_{\chi=0}$
and $\mathscr{L}^{D}$ is the Drazin inverse of $\mathscr{L}$. In particular,
we obtain from Eq. (S54),

$\displaystyle\mathscr{W}_{1}$
$\displaystyle=\Gamma_{I}(\mathscr{V}^{+}-\mathscr{V}^{-}),$ (S58)
$\displaystyle\mathscr{W}_{2}$
$\displaystyle=\Gamma_{I}(\mathscr{V}^{+}+\mathscr{V}^{-}).$ (S59)

###  III.2 Average power

Taking the average of Eq. (S58), we obtain after some algebra the expression
for the current in the particle model,

$\displaystyle\expectationvalue{I}=\langle\mathscr{W}_{1}\rangle=\Gamma_{I}\expectationvalue{\mathscr{N}_{h}-\mathscr{N}_{c}},$
(S60)

also mentioned inline in the main text. Thus, to obtain the steady-state
current, we need the steady-state averages
$\langle\mathscr{N}_{\alpha}\rangle$. Using Eq. (S48), we obtain the
differential equations,

$\displaystyle\derivative{t}\langle\mathscr{N}_{h}\rangle$
$\displaystyle=\kappa_{h}\quantity(\bar{\text{n}}_{h}-\langle\mathscr{N}_{h}\rangle)-\Gamma_{I}\langle\mathscr{N}_{h}-\mathscr{N}_{c}\rangle,$
(S61) $\displaystyle\derivative{t}\langle\mathscr{N}_{c}\rangle$
$\displaystyle=\kappa_{c}\quantity(\bar{\text{n}}_{c}-\langle\mathscr{N}_{c}\rangle)+\Gamma_{I}\langle\mathscr{N}_{h}-\mathscr{N}_{c}\rangle,$
(S62)

As in the previous models, we introduce new variables
$\mathscr{N}_{h}=:(\mathscr{V}+\mathscr{I})/(2\Gamma_{I}),\mathscr{N}_{c}=:(\mathscr{V}-\mathscr{I})/(2\Gamma_{I})$
and write the equations of motion in vector form
$\bm{\sigma}=(\mathscr{I},\mathscr{V})$ as

$\displaystyle\derivative{t}\langle\bm{\sigma}\rangle=G\langle\bm{\sigma}\rangle+F,$
(S63)

where,

$\displaystyle
G=\begin{pmatrix}-\Gamma_{I}-\frac{\kappa_{h}+\kappa_{c}}{2}&-\frac{\kappa_{h}-\kappa_{c}}{2}\\\
-\frac{\kappa_{h}-\kappa_{c}}{2}&-\frac{\kappa_{h}+\kappa_{c}}{2}\end{pmatrix},$
(S64)

and
$F=\Gamma_{I}(\kappa_{h}\bar{\text{n}}_{h}-\kappa_{c}\bar{\text{n}}_{c},\kappa_{h}\bar{\text{n}}_{h}+\kappa_{c}\bar{\text{n}}_{c})$.
Again, we have a linear system for the steady-state averages
$G\langle\bm{\sigma}\rangle+F=0$, and $\langle\mathscr{I}\rangle=\langle
I\rangle$, is given by

$\displaystyle\langle
I\rangle_{p}=\frac{\Gamma_{I}\kappa_{h}\kappa_{c}(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})}{\kappa_{h}\kappa_{c}+\Gamma_{I}(\kappa_{h}+\kappa_{c})}.$
(S65)

Using $\Gamma_{I}=4g^{2}/(\kappa_{h}+\kappa_{c})$ and that $\langle
P\rangle=\langle I\rangle\Delta$, we obtain Eq. (9) of the main text.

###  III.3 Noise

We now concentrate in computing the zero-frequency noise. Using Eq. (S59) the
first term in Eq. (S57) gives after some algebra,

$\displaystyle\langle\mathscr{W}_{2}\rangle$
$\displaystyle=\Gamma_{I}\langle\mathscr{N}_{h}(\mathscr{N}_{c}+1)+\mathscr{N}_{c}(\mathscr{N}_{h}+1)\rangle,$
(S66)
$\displaystyle=\Gamma_{I}\langle\mathscr{N}_{h}+\mathscr{N}_{c}+2\mathscr{N}_{h}\mathscr{N}_{c}\rangle.$
(S67)

Above we have terms that can be readily computed from the equations of motion
for the first moments $\langle\mathscr{N}_{\alpha}\rangle$, solved in the
latter Subsection. However, we also have higher order terms, which demand
addressing the equation of motion of
$\langle\mathscr{N}_{h}\mathscr{N}_{c}\rangle$. 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
$\langle\mathscr{N}_{\alpha}^{2}\rangle,\langle\mathscr{N}_{h}\mathscr{N}_{c}\rangle$.

We start by rewriting Eq. (S63) as an homogeneous equation,

$\displaystyle\derivative{t}\langle\delta\bm{\sigma}\rangle=G\langle\delta\bm{\sigma}\rangle,$
(S68)

where $\delta\bm{\sigma}=\bm{\sigma}-\langle\bm{\sigma}\rangle$. We now
consider the LHS of the above, and that
$\left|\dot{p}\right)=\mathscr{L}\left|p\right)$,

$\displaystyle\derivative{t}\langle\delta\bm{\sigma}\rangle=\left(1\right|\delta\bm{\sigma}\mathscr{L}\left|p\right),~{}\forall\left|p\right).$
(S69)

This is true for all $\left|p\right)$ and, further, not only to normalized
probabilities but to any vector $\left|s\right)$ s.t. $\sum_{i}s_{i}<\infty$.
If we then consider the RHS of Eq. (S68), we conclude that
$\left(1\right|\delta\bm{\sigma}\mathscr{L}=\left(1\right|G\delta\bm{\sigma}$
and similarly,

$\displaystyle\left(1\right|\delta\bm{\sigma}e^{t\mathscr{L}}=\left(1\right|e^{tG}\delta\bm{\sigma},$
(S70)

multiplying both sides by $\left|s\right):=\delta\mathscr{A}\ket{p}$, with
$\mathscr{A}$ generic operator, we obtain a regression theorem,

$\displaystyle\langle\delta\bm{\sigma}e^{t\mathscr{L}}\delta\mathscr{A}\rangle=\langle
e^{tG}\delta\bm{\sigma}\delta\mathscr{A}\rangle.$ (S71)

We now use the above result. In particular, we are interested in the case that
$\mathscr{A}=\mathscr{W}_{1}$. Recalling the second term in Eq. (S57),

$\displaystyle\langle\mathscr{W}_{1}\mathscr{L}^{D}\mathscr{W}_{1}\rangle=-\int_{0}^{\infty}dt\quantity[\expectationvalue{\mathscr{W}_{1}e^{t\mathscr{L}}\mathscr{W}_{1}}-\expectationvalue{\mathscr{W}_{1}}^{2}]=-\int_{0}^{\infty}dt\expectationvalue{\delta\mathscr{W}_{1}e^{t\mathscr{L}}\delta\mathscr{W}_{1}},$
(S72)

where we used that the Drazin inverse can be written as
$\mathscr{L}^{D}=-\int_{0}^{\infty}e^{t\mathscr{L}}(1-\left|p^{*}\right)\left(1\right|)$,
with $\left|p^{*}\right)$ the eigenvector with zero-eigenvalue (steady-state)
of $\mathscr{L}$. Finally, we can use the regression theorem for
$\langle\delta\bm{\sigma}e^{t\mathscr{L}}\delta\mathscr{W}_{1}\rangle$ to
evaluate the correlator,

$\displaystyle\langle\langle\delta\bm{\sigma}\delta\mathscr{W}_{1}\rangle\rangle=\int_{0}^{\infty}dt\expectationvalue{\delta\bm{\sigma}e^{t\mathscr{L}}\delta\mathscr{W}_{1}}=G^{-1}\expectationvalue{\delta\bm{\sigma}\delta\mathscr{W}_{1}},$
(S73)

whose first entry is precisely Eq. (S72). Thus, we obtain the relation,

$\displaystyle\langle\langle
I^{2}\rangle\rangle=\langle\mathscr{W}_{2}\rangle-[G^{-1}\langle\delta\bm{\sigma}\delta\mathscr{W}_{1}\rangle]_{i=1}.$
(S74)

Thus, in order to evaluate Eq. (S72), we need the initial conditions
$\expectationvalue{\delta\bm{\sigma}\delta\mathscr{W}_{1}}$,

$\displaystyle\expectationvalue{\delta\mathscr{V}\delta\mathscr{W}_{1}}$
$\displaystyle=\expectationvalue{\mathscr{V}\mathscr{W}_{1}}-\expectationvalue{\mathscr{V}}\expectationvalue{\mathscr{W}_{1}},$
(S75)

where $\expectationvalue{\mathscr{W}_{1}}=\expectationvalue{I}$ and
$\expectationvalue{\mathscr{V}}$ have been computed already. Turning to the
remaining terms,

$\displaystyle\expectationvalue{\mathscr{W}_{1}\mathscr{W}_{1}}=\Gamma_{I}^{2}\expectationvalue{(\mathscr{V}^{+}-\mathscr{V}^{-})^{2}}=\Gamma_{I}^{2}\expectationvalue{(\mathscr{N}_{h}-\mathscr{N}_{c})^{2}}-2\Gamma_{I}^{2}\expectationvalue{\mathscr{N}_{h}+\mathscr{N}_{c}+2\mathscr{N}_{h}\mathscr{N}_{c}},$
(S76)
$\displaystyle\expectationvalue{\mathscr{V}\mathscr{W}_{1}}=\Gamma_{I}^{2}\expectationvalue{(\mathscr{V}^{+}+\mathscr{V}^{-})(\mathscr{V}^{+}-\mathscr{V}^{-})}=\Gamma_{I}^{2}\expectationvalue{\mathscr{N}^{2}_{h}-\mathscr{N}^{2}_{c}}.$
(S77)

we conclude that the computing the zero-frequency noise reduces to computing
the steady-state averages of
$\langle\mathscr{N}_{\alpha}^{2}\rangle,\langle\mathscr{N}_{h}\mathscr{N}_{c}\rangle$.

Now, we use the rate equation (S48) to obtain the equations of motion,

$\displaystyle\derivative{t}\expectationvalue{\mathscr{N}_{h}^{2}}=-2(\kappa_{h}+\Gamma_{I})\expectationvalue{\mathscr{N}_{h}^{2}}+4\Gamma_{I}\expectationvalue{\mathscr{N}_{h}\mathscr{N}_{c}}+(4\kappa_{h}\bar{\text{n}}_{h}+\kappa_{h}+\Gamma_{I})\expectationvalue{\mathscr{N}_{h}}+\Gamma_{I}\mathscr{N}_{c}+\kappa_{h}\bar{\text{n}}_{h},$
(S78)
$\displaystyle\derivative{t}\expectationvalue{\mathscr{N}_{c}^{2}}=-2(\kappa_{c}+\Gamma_{I})\expectationvalue{\mathscr{N}_{c}^{2}}+4\Gamma_{I}\expectationvalue{\mathscr{N}_{h}\mathscr{N}_{c}}+(4\kappa_{c}\bar{\text{n}}_{c}+\kappa_{c}+\Gamma_{I})\expectationvalue{\mathscr{N}_{c}}+\Gamma_{I}{\mathscr{N}_{h}}+\kappa_{c}\bar{\text{n}}_{c},$
(S79)
$\displaystyle\derivative{t}\expectationvalue{\mathscr{N}_{h}\mathscr{N}_{c}}=\Gamma_{I}\expectationvalue{\mathscr{N}_{h}^{2}}+\Gamma_{I}\expectationvalue{\mathscr{N}_{c}^{2}}-(4\Gamma_{I}+\kappa_{h}+\kappa_{c})\expectationvalue{\mathscr{N}_{h}\mathscr{N}_{c}}+(\bar{\text{n}}_{c}\kappa_{c}-\Gamma_{I})\expectationvalue{\mathscr{N}_{h}}+(\bar{\text{n}}_{h}\kappa_{h}-\Gamma_{I})\expectationvalue{\mathscr{N}_{c}},$
(S80)

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,

$\displaystyle\langle\langle
I^{2}\rangle\rangle_{p}=\frac{2(\bar{\text{n}}_{T}^{2}+2\bar{\text{n}}_{T})g^{2}\bar{\text{n}}_{T}\kappa_{h}\kappa_{c}}{(4g^{2}+\kappa_{h}\kappa_{c})(\kappa_{h}+\kappa_{c})}-8\delta\bar{\text{n}}^{2}g^{2}\kappa_{h}\kappa_{c}$
$\displaystyle\Bigg{(}\frac{192g^{8}(\kappa_{h}-\kappa_{c})^{2}+64g^{6}(\kappa_{h}^{4}+\kappa_{c}^{4}-10\kappa_{h}^{2}\kappa_{c}^{2}-2\kappa_{c}^{3}\kappa_{h}-2\kappa_{h}^{3}\kappa_{c})}{(4g^{2}+\kappa_{h}\kappa_{c})^{3}(\kappa_{h}+\kappa_{c}^{3})(48g^{4}+\kappa_{h}\kappa_{c}(\kappa_{h}+\kappa_{c})^{2})+4g^{2}(\kappa_{h}^{2}+\kappa_{c}^{2}+6\kappa_{h}\kappa_{c})}$
(S81) $\displaystyle-$
$\displaystyle\frac{16g^{4}\kappa_{h}\kappa_{c}(\kappa_{h}+\kappa_{c})^{2}(\kappa_{h}^{2}+\kappa_{c}^{2}+7\kappa_{h}\kappa_{c})-4g^{2}\kappa_{h}^{2}\kappa_{c}^{2}(\kappa_{h}+\kappa_{c})^{4}+\kappa_{h}^{3}\kappa_{c}^{3}(\kappa_{h}+\kappa_{c})^{4}}{(4g^{2}+\kappa_{h}\kappa_{c})^{3}(\kappa_{h}+\kappa_{c}^{3})(48g^{4}+\kappa_{h}\kappa_{c}(\kappa_{h}+\kappa_{c})^{2})+4g^{2}(\kappa_{h}^{2}+\kappa_{c}^{2}+6\kappa_{h}\kappa_{c})}\Bigg{)},$

where we used the notation
$\bar{\text{n}}_{T}=\bar{\text{n}}_{h}+\bar{\text{n}}_{c},\delta\bar{\text{n}}=\bar{\text{n}}_{h}-\bar{\text{n}}_{c}$
and substituted $\Gamma_{I}=4g^{2}/(\kappa_{h}+\kappa_{c})$. Finally,
considering $\langle\langle P^{2}\rangle\rangle=\Delta^{2}\langle\langle
I\rangle\rangle$ and writing the above in terms of equilibrium and shot noise,
we obtain (power) noise,

$\displaystyle\expectationvalue{\expectationvalue{P^{2}}}_{p}$
$\displaystyle=\mathcal{E}(\bar{\text{n}}_{h}(\bar{\text{n}}_{h}+1)+\bar{\text{n}}_{c}(\bar{\text{n}}_{c}+1))-\mathcal{S}_{p}(\bar{\text{n}}_{h}-\bar{\text{n}}_{c})^{2},$
(S82) $\displaystyle\mathcal{E}$
$\displaystyle=4g^{2}\kappa_{h}\kappa_{c}\Delta^{2}\chi,$ (S83)
$\displaystyle\mathcal{S}_{p}$
$\displaystyle=\mathcal{S}+\frac{4(2g)^{6}\kappa_{h}^{3}\kappa_{c}^{3}(12g^{2}+(\kappa_{h}+\kappa_{c})^{2})\Delta^{2}\chi^{3}}{48g^{4}+\kappa_{h}\kappa_{c}(\kappa_{h}+\kappa_{c})^{2}+4g^{2}(\kappa_{h}^{2}+\kappa_{c}^{2}+6\kappa_{h}\kappa_{c})},$
(S84)

where $\chi=[(\kappa_{h}+\kappa_{c})(4g^{2}+\kappa_{h}\kappa_{c})]^{-1}$ and
$\mathcal{S}$ is given in Eq. (S30). Setting $\kappa_{h}=\kappa_{c}=\kappa$ we
obtain Eq. (10c) of the main text.

###  III.4 Limitting regimes $g/\kappa\to 0,~{}g/\kappa\to\infty$.

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

$\displaystyle\mathcal{S}_{p}-\mathcal{S}$
$\displaystyle=\frac{256g^{6}}{\kappa_{h}\kappa_{c}(\kappa_{h}+\kappa_{c})^{2}}+\mathcal{O}(g^{7}),$
(S85) $\displaystyle\mathcal{S}_{p}-\mathcal{S}$
$\displaystyle=\frac{\kappa_{h}^{3}\kappa_{c}^{3}}{g^{2}(\kappa_{h}+\kappa_{c})^{3}}+\mathcal{O}\quantity(\frac{1}{g^{4}}).$
(S86)

##  IV Connection between power noise and work variance

Consider the operator $W(t)=\int_{0}^{t}P(t^{\prime})\differential
t^{\prime}$. Now, we shall show that

$\displaystyle\partial_{t}\quantity(\langle W^{2}(t)\rangle-\langle
W(t)\rangle^{2})=2\real\int_{0}^{t}\differential\tau\langle\delta
P(\tau)\delta P(0)\rangle,$ (S87)

and, thus, zero-frequency noise, Eq. (2), corresponds to the steady-state of
the above, $t\to\infty$.

First,

$\displaystyle\partial_{t}\langle W(t)^{2}\rangle=2\int_{0}^{t}\differential
t^{\prime}\langle P(t^{\prime})\rangle\langle
P(t)\rangle=2\real\int_{0}^{t}\differential t^{\prime}\langle
P(t^{\prime})\rangle\langle P(t)\rangle,$ (S88)

where the $\real$ comes for free, since $P(t)$ is Hermitian. Second,

$\displaystyle\partial_{t}\langle
W^{2}(t)\rangle=\partial_{t}\int_{0}^{t}\int_{0}^{t}\differential
t^{\prime}\differential t^{\prime\prime}\langle
P(t^{\prime})P(t^{\prime\prime})\rangle=\int_{0}^{t}\differential
t^{\prime}\langle\\{P(t),P(t^{\prime})\\}\rangle=2\real\int_{0}^{t}\differential
t^{\prime}\langle P(t)P(t^{\prime})\rangle.$ (S89)

Combining Eqs.(S88, S89),

$\displaystyle\partial_{t}\quantity(\langle W^{2}(t)\rangle-\langle
W(t)\rangle^{2})$ $\displaystyle=2\real\int_{0}^{t}\differential
t^{\prime}\langle\delta P(t)\delta
P(t^{\prime})\rangle=2\real\int_{0}^{t}\differential t^{\prime}\langle\delta
P(t-t^{\prime})\delta P(0)\rangle,$ (S90)
$\displaystyle=2\real\int_{0}^{t}\differential\tau\langle\delta P(\tau)\delta
P(0)\rangle,$ (S91)

where we can safely assume that the correlation function only depends on
$\tau:=t^{\prime}-t$ if we are interested in steady-state averages. Finally,

$\displaystyle\partial_{t}\quantity(\langle W(t)^{2}\rangle-\langle
W(t)\rangle^{2})\Big{|}_{t\to\infty}=2\real\int_{0}^{\infty}\differential\tau\langle\delta
P(\tau)\delta P(0)\rangle,$ (S92)

which corresponds to Eq. (2) of the main text.