# Spontaneous superradiant photon current

Lei Qiao<EMAIL_ADDRESS>Department of Physics, National University of
Singapore, Singapore 117551, Singapore  Jiangbin Gong<EMAIL_ADDRESS>Department of Physics, National University of Singapore, Singapore 117551,
Singapore  Center for Quantum Technologies, National University of Singapore,
Singapore 117543, Singapore

(August 31, 2024)

###### Abstract

This work reports the spontaneous emergence of a photon current in a class of
spin-cavity systems, where an assemble of quantum emitters interact with
distinct photon modes confined in tunneling-coupled cavities. Specifically,
with necessary symmetry breaking, photons in a superradiant phase afforded by
coherent photon-emitter interaction spontaneously flow from a cavity with a
lower resonance frequency to a different cavity with a higher resonance
frequency. Theoretical analysis reveals that cavity dissipation is the key to
alter spin-cavity coherence, which then makes it possible to extract photons
from, and later return photons to the vaccum through the cavities. The
interplay between photon loss and emitter coherence hence sustains a counter-
intuitive steady current of photons between cavities without an external
pumping field.

Quantum phase transitions (QPTs) are a type of phase transitions that occur at
zero temperature and are accessed by varying a nonthermal parameter between
two different quantum phases. As a fundamental phenomenon in physics, QPTs
have been extensively investigated in various fields [1, 2, 3]. For the
second-order QPTs, which are also called continuous phase transitions because
of the associated gapless excitations, one important feature is that the
quantum fluctuations become increasingly important and may cause dramatic
changes in the system’s behavior. Because of the experimental progress for the
scalability and integrability of cavity-quantum-electrodynamics systems with
excellent controllability [4, 5, 6], there has been a growing interest in the
studies of photonic QPTs in coupled-cavity platforms such as the Jaynes-
Cummings Hubbard lattice systems [7, 8, 9, 10, 11, 12] and the Rabi lattice
models [13, 14, 15, 16, 17]. As a typical QPT with light-matter interactions,
superradiant phase transitions in cavity Rabi systems describe the change of
system’s ground states from the normal phase (with cavity-mode in the vacuum
state and atoms in their ground states) to the symmetry broken phase
characterized by two degenerate ground states with atomic and photonic
excitations [18, 19, 20, 21, 22, 23, 24]. More related to this work, in the
presence of dissipation for realistic experimental setups, the characteristics
of superradiant phase transitions can be significantly affected [25, 26, 27,
28, 29, 30, 31, 32], such as anomalous multicritical phenomena [28] and self-
oscillating pump [32] in atom-cavity systems. Recently, the superradiant phase
transition with a tricritical point was proposed by introducing a pump
squeezed field [33, 34, 35]. In addition, the finite-component
multicriticality at the superradiant QPT was observed in the qubit-boson
system via engineered qubit biases [36].

This work aims to show one striking consequence of superradiant phase
transition in a spin-cavity system with photon loss. Specifically, the quantum
coherence arising from the interplay of spin-cavity coupling and dissipation
induces a sustained photon current between different tunneling-coupled
cavities, without the use of any external pump field. In essence, in the
superradiant phase, dissipation is found to act as a pump that breaks the
time-reversal symmetry and generates a steady-state current in a class of
spin-cavity systems that involve quantum emitters coupled to multiple
cavities. As seen below, to generate the said spontaneous photon current, some
other symmetries of the system are also to be broken, either spontaneously or
by design. Interestingly, photons flow in a sustained way, from a cavity with
lower resonance frequency to a different cavity with a higher resonance
frequency. Such spontaneous photon current is intriguing because it challenges
our naive understanding of what dissipation can really do. Our results
indicate that photon loss applied to a cavity can actually assist the cavity
to first acquire photons from vacuum and then pass them to other cavities with
ultimately giving them back to its surrounding.

The system we consider here consists of a collection of $N$ identical two-
level quantum emitters with resonant frequency $\Omega$, which interact with
the bosonic modes confined with multiple tunneling-coupled cavities, with Fig.
1 showcasing an example with three cavities in a triangular configuration.
Assuming $a_{n}^{{\dagger}}$ and $a_{n}$ are the creation and annihilation
operators of the $n$th cavity mode respectively, the tight-binding Hamiltonian
of such spin-cavity system is given by

$\displaystyle H$
$\displaystyle=\sum_{n=1}^{N_{c}}\hbar\omega_{n}a_{n}^{{\dagger}}a_{n}+\sum_{n=1}^{N_{c}}\hbar
J\left(a_{n+1}^{{\dagger}}a_{n}+a_{n}^{{\dagger}}a_{n+1}\right)$
$\displaystyle+\hbar\Omega S_{z}+\sum_{n=1}^{N_{c}}\frac{2\hbar
g}{\sqrt{N}}\left(a_{n}+a_{n}^{{\dagger}}\right)S_{x}\text{,}$ (1)

where $\omega_{n}$ is the resonance frequency of photon mode in the $n$th
cavity with a periodic boundary condition $a_{N_{c}+1}=a_{1}$ ($N_{c}=3$ as an
example to elaborate on our theoretical calculations below), and $J$ is the
photon hopping strength. The collective spin operators
$S_{\alpha}=\sum_{j=1}^{N}\sigma_{\alpha}^{j}/2$ ($\alpha=x,y,z$) are defined
by using the individual Pauli spin operators $\sigma_{\alpha}^{j}$, where $j$
denotes the $j$th two-level emitter. At the $n$th cavity site, the local spin-
cavity coupling strength is characterized by $g$. Due to the presence of the
counter-rotating-wave (CRW) part in the interaction between emitters and
photons, the continuous $U\left(1\right)$ symmetry of the system is broken.
However, the Hamiltonian $H$ commutes with the parity operator
$\Pi=\exp(i\pi\hat{P})$ where
$\hat{P}=S_{z}+\sum_{n=1}^{N_{c}}a_{n}^{{\dagger}}a_{n}$, and thus the system
possesses a global $Z_{2}$ symmetry, the breaking of which signals the optical
switching from a normal phase to two degenerate superradiant phases. It is
important to note that the system can directly extract photons from the vacuum
due to the CRW terms $a_{n}^{{\dagger}}S_{+}$ in the Hamiltonian $H$, where
$S_{\pm}=S_{x}\pm iS_{y}$ are the collective spin ladder operators. Because
the spin system here is coupled to the tunneling-coupled cavities at multiple
locations, such a configuration may be regarded as cavities plus “giant atoms”
[37, 38, 39, 40], which has been realized at different experimental platforms
such as ferromagnetic spin ensemble and superconducting artificial atoms [41,
42, 43].

Figure 1: Schematic of a spin-cavity setup. A loop of three cavities with
nearest-neighbor tunneling rate $J$ is also coupled to the same set of two-
level emitters.

To account for the impact of dissipation, the dynamical evolution of the spin-
cavity system can be described by a density matrix $\hat{\rho}\left(t\right)$
satisfying the following Lindblad master equation

$\frac{\partial\hat{\rho}}{\partial
t}=-\frac{i}{\hbar}\left[H\left(t\right),\hat{\rho}\right]+\kappa\sum_{n}\mathcal{D}\left[a_{n}\right]\hat{\rho}+\frac{\gamma}{N}\mathcal{D}\left[S_{-}\right]\hat{\rho}\text{,}$
(2)

where the Lindblad dissipators are
$\mathcal{D}[a_{n}]\hat{\rho}=2a_{n}\hat{\rho}a_{n}^{{\dagger}}-\\{a_{n}^{{\dagger}}a_{n},\hat{\rho}\\}$
and
$\mathcal{D}[S_{-}]\hat{\rho}=2S_{-}\hat{\rho}S_{+}-\\{S_{+}S_{-},\hat{\rho}\\}$.
Clearly, the first term on the right-hand side of Eq. (2) represents the
standard Hamiltonian evolution while the other two terms describe the
dissipation of the cavities and the spin system respectively. Given that some
ferromagnetic spin ensembles or artificial atom systems like superconducting
qubits can be controlled to work with a dissipation rate in the kilohertz
regime, much lower in comparison with other physical quantities in megahertz
and even gigahertz regime [44, 45, 46], here we neglect the dissipation rate
of two-level emitters and hence set $\gamma=0$ in our following analysis. This
also helps us to focus on the impact of cavity dissipation at the rate
$\kappa$. It should be noted that only for free cavity dynamics, such cavity
dissipation rate means an exponential photon loss at rate $\kappa$ on the
average. When there is coherent interaction between emitters and cavities, it
is far from obvious what this cavity dissipation rate can lead to.

The complex order parameters capturing the superradiant phase transition are
$\left\langle a_{n}\right\rangle=\sqrt{N}\alpha_{n}$ where
$\alpha_{n}=\alpha_{n,\operatorname{Re}}+i\alpha_{n,\operatorname{Im}}$,
namely, its real and imaginary parts. To study the phase transition, we apply
a standard mean-field treatment to Eq. (2), yielding

$i\frac{d\alpha_{n}}{dt}=\omega_{n}\alpha_{n}+J\left(\alpha_{n-1}+\alpha_{n+1}\right)+2gX-i\kappa\alpha_{n}\text{,}$
(3) $\frac{dX}{dt}=-\Omega Y\text{,}$ (4) $\frac{dY}{dt}=\Omega
X-Z\sum_{n}2g\left(\alpha_{n}+\alpha_{n}^{\ast}\right)\text{,}$ (5)
$\frac{dZ}{dt}=Y\sum_{n}2g\left(\alpha_{n}+\alpha_{n}^{\ast}\right)\text{,}$
(6)

where $\left\langle S_{x}\right\rangle=NX$, $\left\langle
S_{y}\right\rangle=NY$, $\left\langle S_{z}\right\rangle=NZ$. Since
$[S_{\alpha},\mathbf{S}^{2}]=0$ for $\alpha=x$, $y$, $z$, one observes that
the Hamiltonian in Eq. (1) commutes with $\mathbf{S}^{2}$. As such the total
length of the collective spin is conserved. This relation guarantees that the
equations of motion for the collective spin properties satisfies
$\partial\mathbf{S}^{2}/\partial t=0$, even with cavity dissipation. One
therefore always has $X^{2}+Y^{2}+Z^{2}=1/4$ when analysing the dynamics.

Figure 2: Order parameters $\alpha_{n}$ as functions of the cavity dissipation
rate $\kappa/\Omega$ and the coupling strength $g$, with (a) the real part
$\alpha_{n,\operatorname{Re}}$ and (b) the imaginary part
$\alpha_{n,\operatorname{Im}}$. Note that $\alpha_{n}$ for $n=1,2,3$ are equal
to each other due to mode symmetry here. Other parameters are $J/\Omega=0.1$
and $\omega_{c}/\Omega=0.5$.

To gain some insights into the superradiant transition, we first consider the
system under mode symmetry, namely,
$\omega_{1}=\omega_{2}=\omega_{3}=\omega_{c}$. Under the steady-state
condition, the dynamical equations in Eq. (3)-(6) can be solved analytically.
For the trivial phase, all the two-level emitters are in their ground states
and there is no photon in the cavities at the steady state. This corresponds
to $Z=-1/2$ and $\alpha_{n,\operatorname{Re}}=\alpha_{n,\operatorname{Im}}=0$.
It is found that if the spin-cavity coupling strength $g$ exceeds the
following critical value

$g_{c}=\frac{1}{2}\sqrt{\frac{\Omega}{3}\left(\omega_{c}+2J+\frac{\kappa^{2}}{\omega_{c}+2J}\right)}\text{,}$
(7)

then the nontrivial solutions with nonzero $\alpha_{n,\operatorname{Re}}$ and
$\alpha_{n,\operatorname{Im}}$ can be obtained as
$\alpha_{1,\operatorname{Re}}=\alpha_{2,\operatorname{Re}}=\alpha_{3,\operatorname{Re}}=A/3$
and
$\alpha_{1,\operatorname{Im}}=\alpha_{2,\operatorname{Im}}=\alpha_{3,\operatorname{Im}}=B/3$,
where $A$ and $B$ are given by

$A=\pm\frac{\Omega}{4g}\sqrt{\frac{1}{4Z^{2}}-1}\text{, }B=\frac{\kappa
A}{\omega_{c}+2J}\text{.}$ (8)

The corresponding steady-state values of $X$, $Y$ and $Z$ are given by
$X=\pm\sqrt{1/4-Z^{2}}$, $Y=0$ and

$Z=-\frac{\Omega}{24g^{2}}\left(\omega_{c}+2J+\frac{\kappa^{2}}{\omega_{c}+2J}\right)\text{.}$
(9)

The stability of these nontrivial solutions can be checked by expanding the
order parameters to the first order with small fluctuations (see Supplementary
Material [47]). The positivity condition for the real parts of the eigenvalues
of the dynamical matrix satisfied by these fluctuations guarantees that the
order parameters shall evolve to its steady-state values when they are
perturbed. The above two superradiant solutions for $\alpha_{n}$ and $X$ have
opposite values, thus choosing any one of them spontaneously breaks the
$Z_{2}$ symmetry. The critical coupling strength $g_{c}$ defines a second-
order phase transition.

From the above superradiant phase solution, the cavity coupling $J$ induces an
effective shift $2J$ of the cavity frequency $\omega_{c}$ because
$\omega_{c}+2J$ always appear together. For small $\kappa$, the critical value
$g_{c}$ will increase as $J$ becomes stronger. Furthermore, compared to the
case of one cavity with a collection of identical two-level emitters [49],
$g_{c}$ here becomes smaller as it is scaled by a factor $1/\sqrt{3}$ if both
$J$ and $\kappa$ are much less than $\omega_{c}$. This indicates that coupling
the emitters with multiple cavities has an interesting impact: it makes the
superradiant phase transition more likely to occur. Indeed, one can obtain
that as the number $N_{c}$ of cavities increases, $g_{c}$ will decrease with
the scaling $1/\sqrt{N_{c}}$. This fact provides us a way to reduce the
critical coupling strength $g_{c}$ for the superradiant phase transition. In
Fig. 2, we plot $\alpha_{n,\operatorname{Re}}$ and
$\alpha_{n,\operatorname{Im}}$ as functions of $g$ and the cavity dissipation
rate $\kappa$ in panel (a) and panel (b), respectively. The white areas
represent the normal phase. The areas with $\alpha_{n,\operatorname{Re}}\neq
0$ and $\alpha_{n,\operatorname{Im}}\neq 0$ indicate the superradiant phase.
The boundary between these two phases corresponds to $Z=-1/2$. It can also be
seen from Fig. 2 that the superradiant phase can be gradually restored back to
the normal phase if we continue to increase $\kappa$ while fixing other system
parameters. Among all these results discussed here, of central importance to
our key results below is that the cavity dissipation can induce imaginary
values of the system’s order parameter $\alpha_{n}$. As shown in Eq. (8), the
superradiant phase will necessarily have nonzero
$\alpha_{n,\operatorname{Im}}$ if $\kappa\neq 0$.

Figure 3: (a) Order parameters $\left|\alpha_{n}\right|$ as a function of the
coupling strength $g$ with $\Delta/\Omega=0.5$. (b) Photon number fluctuations
as a function of $g$ with $\Delta/\Omega=0.5$. (c) Photon current
$\langle\hat{I}\rangle/N$ as functions of the difference frequency $\Delta$
and the coupling strength $g$. Other parameters are $J/\Omega=0.2$,
$\omega_{c}/\Omega=0.1$ and $\kappa/\Omega=0.3$.

We are now ready to break the symmetry between the cavity modes in terms of
their different resonance frequencies. In the case of using three cavities as
an example, one can still fortunately find analytical solutions corresponding
to steady-state superradiant phases, characterized by
$Z=\Omega/(4g\tilde{A})$, where
$\tilde{A}=\sum_{n}g(\tilde{\alpha}_{n}+\tilde{\alpha}_{n}^{\ast})/2$ with

$\tilde{\alpha}_{n}=-\frac{2(J+i\kappa-\omega_{n-1})\left(J+i\kappa-\omega_{n+1}\right)}{2J^{3}+J^{2}\left(3i\kappa-\omega_{\mathrm{tot}}\right)+i\prod_{n^{\prime}}\left(\kappa+i\omega_{n^{\prime}}\right)}\text{.}$
(10)

Here $\omega_{\mathrm{tot}}$ is defined as
$\omega_{\mathrm{tot}}=\sum_{n^{\prime}}\omega_{n^{\prime}}$. The collective
spin along the other two directions are found to be $X=\pm\sqrt{1/4-Z^{2}}$,
$Y=0$. The order parameters for cavity modes are
$\alpha_{1}=gX\tilde{\alpha}_{1}$, $\alpha_{2}=gX\tilde{\alpha}_{2}$,
$\alpha_{3}=gX\tilde{\alpha}_{3}$. Again, the presence of cavity dissipation
leads to complex $\alpha_{n}$, but now in different manners due to mode
symmetry breaking. The critical strength of $g$ for the superradiant phase
occurs at
$g_{c}=\\{-\Omega/[\sum_{n}(\tilde{\alpha}_{n}+\tilde{\alpha}_{n}^{\ast})]\\}^{1/2}$.
For the case of one cavity coupled by a collection of two-level emitters, the
critical coupling strength $g_{c}$ is proportional to the square root of the
resonant frequency of the cavity [49]. Curiously, despite having different
resonance frequencies in general, all the three cavity modes here share this
common critical coupling strength $g_{c}$. In Fig. 3(a), for a fixed value of
$\Delta$ we present how the nonzero values of the order parameters for the
three cavity modes arise. It is seen that they make a sharp phase transition
at the same critical value of $g_{c}$. To further show the smoking-gun
evidence of second-order phase transitions underlying the photon current to be
highlighted below, we also examine the quantum fluctuations of operators
$a_{n}$ and $a_{n}^{{\dagger}}$ in each cavity on top of the mean-field steady
states. By defining $a_{n}=\alpha_{n}\sqrt{N}+c_{n}$ with $c_{n}$ being the
cavity fluctuation operator [47], the photon number fluctuations $\langle
c_{n}^{{\dagger}}c_{n}\rangle$ is presented in Fig. 3(b) after the system has
reached its steady state. Clearly then, as the system makes transitions from
the normal phase to the superradiant phase with cavity dissipation, the photon
number fluctuations continuously diverge when approaching either side of the
phase boundary, thus confirming the nature of a continuous second-order
transition.

The broken time-reversal symmetry due to dissipation yields different complex
order parameters for the different cavities and this allows us to investigate
if a steady current at the superadiant phase can be established. To that end
we define the photon current operator as
$\hat{I}=iJ[(a_{1}^{{\dagger}}a_{2}+a_{2}^{{\dagger}}a_{3}+a_{3}^{{\dagger}}a_{1})-h.c.]$
as a discretized form of the continuous current operator [47, 50, 51, 52]. The
expression of $\hat{I}$ hints that we need complex order parameters to
possibly have a net nonzero current. To quantitatively characterize asymmetry
between the cavities, we assume $\omega_{1}=\omega_{c}$,
$\omega_{2}=\omega_{c}+\Delta$, $\omega_{3}=\omega_{c}+2\Delta$, with $\Delta$
being the asymmetry parameter.

Fig. 3(c) shows the phase diagram for $\langle\hat{I}\rangle$ in the
$\Delta-g$ plane with dissipation rate $\kappa=0.3\Omega$. First of all, there
is no photon current under $\Delta\rightarrow 0$, regardless of the spin-
cavity strength $g$. That is, even when the system enters into the
superradiant phase like the previous case with mode symmetry, there is no
photon current. Next we increase the asymmetry parameter $\Delta$. It is
observed that when the system is in the superradiant phase,
$\langle\hat{I}\rangle$ at a fixed $g$ emerges, increases first with $\Delta$.
In this regime, the increasing $\Delta$ as the difference in cavity
frequencies boosts the inequality of order parameters $\alpha_{n}$ between two
neighbouring cavities, leading to an enhancement in the photon current.
However, as $\Delta$ continues to increase, the current begins to decline
after it reaches its maximum value. The critical point $\Delta_{c}$ where the
current vanishes is found to be exactly at the boundary between the normal
phase and the superradiant phase, which can be determined by the equation
$\Omega+\sum_{n}g^{2}\left(\tilde{\alpha}_{n}+\tilde{\alpha}_{n}^{\ast}\right)=0$.
In this second regime where the photon current decreases as a consequence of
an increasing $\Delta$, another factor becomes dominant: if the resonance
frequencies $\omega_{2}$ and $\omega_{3}$ are much higher than $\omega_{1}$,
there is a severe energy mismatch for the photons to tunnel. In the end, as
$\Delta$ reaches $\Delta_{c}$, it causes a fundamental change to the steady
state, i.e., back to the normal phase that can no longer accomodate the photon
current.

Figure 4: (a) Current $\langle\hat{I}_{n,n+1}\rangle$ as a function of
$\Delta$ with $\kappa/\Omega=0.3$. (b) Current $\langle\hat{I}_{n,n+1}\rangle$
as a function of $\kappa$ with $\Delta/\Omega=0.2$. (c) Current
$\langle\hat{I}_{s,n}\rangle$ as a function of $\Delta$ with
$\kappa/\Omega=0.3$. (d) Time evolution of $\left|\alpha_{n}\right|$ for
$\kappa/\Omega=0.3$ and $\Delta/\Omega=0.5$, with the initial values
$\alpha_{n}(0)=0$, $Y\left(0\right)=0$ and $Z\left(0\right)=-0.499$. The time
is in units of $1/\Omega$. Other parameters are $J/\Omega=0.2$,
$\omega_{c}/\Omega=0.1$ and $g/\Omega=0.35$.

To digest the generated photon current in more details, we now examine the
local photon current operator
$\hat{I}_{n,n+1}=iJ(a_{n}^{{\dagger}}a_{n+1}-h.c)$ from cavity $n$ to cavity
$n+1$. In Fig. 4(a) and (b), the local currents $\langle\hat{I}_{1,2}\rangle$,
$\langle\hat{I}_{2,3}\rangle$ and $\langle\hat{I}_{3,1}\rangle$ are plotted as
a function of $\Delta$ and $\kappa$ respectively. It is seen that
$\langle\hat{I}_{1,2}\rangle$ and $\langle\hat{I}_{2,3}\rangle$ are positive
while $\langle\hat{I}_{3,1}\rangle$ is negative. This indicates that the
nonzero net current $\langle
iJ[(a_{1}^{{\dagger}}a_{2}+a_{2}^{{\dagger}}a_{3}+a_{3}^{{\dagger}}a_{1})-h.c.]\rangle$
analyzed above is not because of a net chiral photon flow, but because of an
imbalance between inter-cavity currents. That is, there are now two routes of
photon flow: via the first route, photons flow from cavity 1 to cavity 3; via
the second route, photons flow from cavity 1 to cavity 2, and then from cavity
2 to cavity 3. Interestingly, one trend is clear: overall photons tend to flow
from a cavity with lower resonance frequency to another cavity with higher
resonance frequency. Furthermore, a simple analysis based on the expression of
the inter-cavity currents suffices to show that the non-chiral nature of the
currents is the rule.

The above detailed characteristics of the photon currents call for a more
careful analysis in terms of local current conservation. Due to the presence
of the rotating-wave (RW) terms and the CRW terms in the system’s Hamiltonian,
the spin-cavity interaction can either produce a photon while deexciting one
emitter (plus the opposite process) or produce a photon while also exciting
one emitter (plus the opposite process). These two processes lead to the
“photon currents” between the emitters and a local cavity, namely,
$\langle\hat{I}_{s,n}\rangle=ig\langle
a_{n}S_{+}-a_{n}^{{\dagger}}S_{-}\rangle/\sqrt{N}$ for the RW terms and
$\langle\hat{I}_{s,n}^{c}\rangle=ig\langle
a_{n}S_{-}-a_{n}^{{\dagger}}S_{+}\rangle/\sqrt{N}$ for the CRW terms [47].
Because the spin component $Y$ in the $y$-axis direction is always zero,
$\langle\hat{I}_{s,n}\rangle$ is found to be equal to
$\langle\hat{I}_{s,n}^{c}\rangle$ for the steady states at the mean-field
level. Fig. 4(c) shows the current $\langle\hat{I}_{s,n}\rangle$ between the
emitters and cavity $n$. One sees that even for $\Delta=0$,
$\langle\hat{I}_{s,n}\rangle$ is nonzero (at the same value) for all the
cavities. As $\Delta$ increases, the currents $\langle\hat{I}_{s,3}\rangle$
and $\langle\hat{I}_{s,2}\rangle$ sequentially change the flow directions.
Given all these details, one wonders if there is a conservation law connecting
all the currents. This is indeed true because the Kirchhoff’s law for the
photon currents can be satisfied at the mean-field level, which assumes the
form
$\langle\hat{I}_{n-1,n}\rangle-\langle\hat{I}_{n,n+1}\rangle+\langle\hat{I}_{s,n}\rangle+\langle\hat{I}_{s,n}^{c}\rangle-\langle
I_{n,d}\rangle=0$ for each cavity as a node. Here $\langle
I_{n,d}\rangle=2\kappa\langle a_{n}^{{\dagger}}a_{n}\rangle$ is the
dissipation photon current [47]. For the spin system as a node, Kirchhoff’s
law is also manifested by
$\langle\hat{I}_{s,n}\rangle=\langle\hat{I}_{s,n}^{c}\rangle$. It follows that
the specific behavior of the asymmetric inter-cavity photon currents should be
always digested together with the photon flows between the emitters and local
cavities as well as the photon dissipation. In particular, the inter-cavity
photon currents can be said to start with extracting photons from the vaccum
and eventually return the photons to the environment through dissipation.

Finally, to demonstrate that the above-analyzed photon currents are truly
spontaneous, we specifically present a computational example below to show why
a spontaneous photon current can be expected as the system evolves by itself
from an arbitrary initial state. To that end we show the time evolution of
$\left|\alpha_{n}(t)\right|$ in the regime of superradiant phase but with the
initial values of $\alpha_{n}(0)=0$, $Y\left(0\right)=0$ and
$Z\left(0\right)=-0.499$ in Fig. 4(d). As time evolves, one sees that
$\alpha_{n}(t)$ indeed spontaneously approaches its steady-state values, as
predicted by our mean-field treatment and also matching our analytical results
[the dots in Fig. 4(d), obtained by use of $\alpha_{n}=gX\tilde{\alpha}_{n}$
with Eq. (10)]. Because nonzero imaginary values of
$\left|\alpha_{n}(t)\right|$ are a necessary condition for inter-cavity photon
currents, the results here further indicate the importance of the spontaneous
symmetry-breaking superradiant phase in generating the spontaneous photon
current.

In summary, we have discovered a spontaneous superradiant photon current
between multiple cavities coupled by a common set of two-level emitters. By
introducing cavity loss and hence breaks the time reversal symmetry, steady
photon current between different cavities can be generated in the superradiant
phase when we break the cavity mode symmetry. Dissipation hence acts on a pump
to generate the photon currents going from a cavity with lower resonance
frequency to another with higher resonance frequency. As to experimental
realizations, the platform consisting of superconducting coplanar cavities and
a yttrium iron garnet sphere is a highly promising system. For tunneling-
coupled cavities, the hopping energy $J/(2\pi)$ ranges from $20$ MHz to $730$
MHz in the existing experimental setups [53, 54, 55, 56, 57, 58]. The
superconducting cavities have frequencies between $4$ and $15$ GHz [45, 59].
The spin-cavity coupling strength $g=g_{0}\sqrt{N}$ can been enhanced by the
number of emitters $N$, the largest reported coupling strength $g/(2\pi)$ for
yttrium iron garnet system can reach $3.06$ GHz [60]. The ferromagnetic spin
ensemble of yttrium iron garnet is an easily tunable system with the magnon
mode frequency ranging from a few hundred MHz to several tens of GHz. Thus,
spontaneous superradiant photon current proposed here is within the reach of
these achieved technologies [61, 62, 63, 64].

Acknowledgments.— J.G. was supported by the National Research Foundation,
Singapore and A*STAR under its CQT Bridging Grant.

## References

  * [1] S. Sachdev, Quantum Phase Transitions, 2nd Edition (Cambridge University Press, 2011).
  * [2] D. Rossini and E. Vicari, Coherent and dissipative dynamics at quantum phase transitions, Phys. Rep. 936, 1 (2021).
  * [3] S. L. Sondhi, S. M. Girvin, J. P. Carini, and D. Shahar, Continuous quantum phase transitions, Rev. Mod. Phys. 69, 315 (1997).
  * [4] M. J. Hartmann, F. G. S. L. Brandao, and M. B. Plenio, Strongly interacting polaritons in coupled arrays of cavities, Nat. Phys. 2, 849 (2006).
  * [5] M. Notomi, E. Kuramochi, and T. Tanabe, Large-scale arrays of ultrahigh-Q coupled nanocavities, Nat. Photon. 2, 741 (2008).
  * [6] A. A. Houck, H. E. Tureci, and J. Koch, On-chip quantum simulation with superconducting circuits, Nat. Phys. 8, 292 (2012).
  * [7] M. J. Hartmann, F. G. S. L. Brandao, and M. B. Plenio, Quantum many-body phenomena in coupled cavity arrays, Laser Photonics Rev. 2, 527 (2008).
  * [8] A. D. Greentree, C. Tahan, J. H. Cole, and L. C. Hollenberg, Quantum phase transitions of light, Nat. Phys. 2, 856 (2006).
  * [9] J. Koch, A. A. Houck, K. L. Hur, and S. M. Girvin, Time-Reversal-Symmetry Breaking in Circuit-QED-Based Photon Lattices, Phys. Rev. A 82, 043811 (2010).
  * [10] C. Noh and D. G Angelakis, Quantum simulations and many-body physics with light, Rep. Prog. Phys. 80, 016401 (2016).
  * [11] M. J. Hwang and M. B. Plenio, Quantum Phase Transition in the Finite Jaynes-Cummings Lattice Systems, Phys. Rev. Lett. 117, 123602 (2016).
  * [12] H. Zheng and Y. Takada, Importance of counter-rotating coupling in the superfluid-to-Mott-insulator quantum phase transition of light in the Jaynes-Cummings lattice, Phys. Rev. A 84, 043819 (2011).
  * [13] M. Schiró, M. Bordyuh, B. Öztop, and H. E. Türeci, Phase transition of light in cavity QED lattices, Phys. Rev. Lett. 109, 053601 (2012).
  * [14] M. J. Hwang, R. Puebla, and M. B. Plenio, Quantum Phase Transition and Universal Dynamics in the Rabi Model, Phys. Rev. Lett. 115, 180404 (2015).
  * [15] D. Lv, S. An, Z. Liu, J.-N. Zhang, J. S. Pedernales, L. Lamata, E. Solano, and K. Kim, Quantum Simulation of the Quantum Rabi Model in a Trapped Ion, Phys. Rev. X 8, 021027 (2018).
  * [16] Q. X. Mei, B. W. Li, Y. K. Wu, M. L. Cai, Y. Wang, L. Yao, Z. C. Zhou, and L. M. Duan, Experimental Realization of the Rabi-Hubbard Model with Trapped Ions, Phys. Rev. Lett. 128, 160504 (2022).
  * [17] J. Zhao and M. J. Hwang, Frustrated Superradiant Phase Transition, Phys. Rev. Lett. 128, 163601 (2022).
  * [18] K. Hepp and E. H. Lieb, On the superradiant phase transition for molecules in a quantized radiation field: The Dicke Maser model, Ann. Phys. (N. Y) 76, 360 (1973).
  * [19] Y. K. Wang and F. T. Hioe, Phase transition in the Dicke model of superradiance, Phys. Rev. A 7, 831 (1973).
  * [20] A. Baksic and C. Ciuti, Controlling discrete and continuous symmetries in “superradiant” phase transitions with circuit QED systems, Phys. Rev. Lett. 112, 173601 (2014).
  * [21] A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Ultrastrong Coupling between Light and Matter, Nat. Rev. Phys. 1, 19 (2019).
  * [22] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Ultrastrong Coupling Regimes of Light-Matter Interaction, Rev. Mod. Phys. 91, 025005 (2019).
  * [23] M. Feng, Y. P. Zhong, T. Liu, L. L. Yan, W. L. Yang, J. Twamley, and H. Wang, Exploring the quantum critical behaviour in a driven Tavis-Cummings circuit, Nat. Commun. 6, 7111 (2015).
  * [24] P. Nataf, T. Champel, G. Blatter, and D. M. Basko, Rashba cavity QED: a route towards the superradiant quantum phase transition, Phys. Rev. Lett. 123, 207402 (2019).
  * [25] J. Keeling, M. J. Bhaseen, and B. D. Simons, Collective Dynamics of Bose-Einstein Condensates in Optical Cavities, Phys. Rev. Lett. 105, 043001 (2010).
  * [26] F. Piazza and H. Ritsch, Self-Ordered Limit Cycles, Chaos, and Phase Slippage with a Superfluid inside an Optical Resonator, Phys. Rev. Lett. 115, 163601 (2015).
  * [27] H. J. Carmichael, Breakdown of Photon Blockade: A Dissipative Quantum Phase Transition in Zero Dimensions, Phys. Rev. X 5, 031028 (2015).
  * [28] M. Soriente, T. Donner, R. Chitra, and O. Zilberberg, Dissipation-Induced Anomalous Multicritical Phenomena, Phys. Rev. Lett. 120, 183603 (2018).
  * [29] M. J. Hwang, P. Rabl, and M. B. Plenio, Dissipative phase transition in the open quantum Rabi model, Phys. Rev. A 97, 013825 (2018).
  * [30] R. Gutiérrez-Jáuregui and H. J. Carmichael, Dissipative quantum phase transitions of light in a generalized Jaynes-Cummings-Rabi model, Phys. Rev. A 98, 023804 (2018).
  * [31] N. Dogra, M. Landini, K. Kroeger, L. Hruby, T. Donner, and T. Esslinger, Dissipation induced structural instability and chiral dynamics in a quantum gas, Science 366, 1496 (2019).
  * [32] D. Dreon, A. Baumgärtner, X. Li, S. Hertlein, T. Esslinger, and T. Donner, Self-oscillating pump in a topological dissipative atom–cavity system, Nature (London) 608, 494 (2022).
  * [33] C. J. Zhu, L. L. Ping, Y. P. Yang, and G. S. Agarwal, Squeezed Light Induced Symmetry Breaking Superradiant Phase Transition, Phys. Rev. Lett. 124, 073602 (2020).
  * [34] G. Q. Zhang, Z. Chen, W. Xiong, C. H. Lam, and J. Q. You, Parity-symmetry-breaking quantum phase transition via parametric drive in a cavity magnonic system, Phys. Rev. B 104, 064423 (2021).
  * [35] Y. Qin, S. C. Li, K. Li, and J. J. Song, Controllable quantum phase transition in a double-cavity magnonic system, Phys. Rev. B 106, 054419 (2022).
  * [36] H. J. Zhu, K. Xu, G. F. Zhang, and W. M. Liu, Finite-Component Multicriticality at the Superradiant Quantum Phase Transition, Phys. Rev. Lett. 125, 050402 (2020).
  * [37] A. F. Kockum, G. Johansson, and F. Nori, Decoherencefree interaction between giant atoms in waveguide quantum electrodynamics, Phys. Rev. Lett. 120, 140404 (2018).
  * [38] L. Guo, A. F. Kockum, F. Marquardt, and G. Johansson, Oscillating bound states for a giant atom, Phys. Rev. Research 2, 043014 (2020).
  * [39] Wei Zhao and Zhihai Wang, Single-photon scattering and bound states in an atom-waveguide system with two or multiple coupling points, Phys. Rev. A 101, 053855 (2020).
  * [40] A. F. Kockum, P. Delsing, and G. Johansson, Designing frequency-dependent relaxation rates and lamb shifts for a giant artificial atom, Phys. Rev. A 90, 013837 (2014).
  * [41] Z.-Q. Wang, Y.-P. Wang, J. Yao, R.-C. Shen, W.-J. Wu, J. Qian, J. Li, S.-Y. Zhu, and J. You, Giant spin ensembles in waveguide magnonics, Nature Communications 13, 7580 (2022).
  * [42] B. Kannan, M. Ruckriegel, D. Campbell, A. F. Kockum, J. Braumüller, D. Kim, M. Kjaergaard, P. Krantz, A. Melville, B. M. Niedzielski, A. Vepsäläinen, R. Winik, J. Yoder, F. Nori, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Waveguide quantum electrodynamics with superconducting artificial giant atoms, Nature (London) 583, 775 (2020).
  * [43] G. Andersson, B. Suri, L. Guo, T. Aref, and P. Delsing, Non-exponential decay of a giant artificial atom, Nat. Phys. 15, 1123 (2019).
  * [44] X. Gu, A. F. Kockum, A. Miranowicz, Y. X. Liu, and F. Nori, Microwave photonics with superconducting quantum circuits, Phys. Rep. 718–719, 1 (2017).
  * [45] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Circuit quantum electrodynamics, Rev. Mod. Phys. 93, 025005 (2021).
  * [46] N. Meher and S. Sivakumar, A review on quantum information processing in cavities, Eur. Phys. J. Plus 137, 985 (2022).
  * [47] In this Supplemental Material, in which the Refs. [48] is cited, we provide details for the stability analysis, the derivation of current operators and the calculation for fluctuation Hamiltonian.
  * [48] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series and Products (Academic, Orlando, 1980).
  * [49] C. Emary and T. Brandes, Chaos and the Quantum Phase Transition in the Dicke Model, Phys. Rev. E 67, 066203 (2003).
  * [50] Y. Y. Zhang, Z. X. Hu, L. B. Fu, H. G. Luo, H. Pu, and X. F. Zhang, Quantum Phases in a Quantum Rabi Triangle, Phys. Rev. Lett. 127, 063602 (2021).
  * [51] K. A. Jung and J. Liu, Stationary current between hybrid vibrational polaritonic systems, Phys. Rev. B 108, L081404 (2023).
  * [52] A. Messiah, Quantum Mechanics (Dover Publications, New York, 1995).
  * [53] N. M. Sundaresan, R. Lundgren, G. Zhu, A. V. Gorshkov, and A. A. Houck, Interacting qubit-photon bound states with superconducting circuits, Phys. Rev. X 9, 011021 (2019).
  * [54] V. S. Ferreira, J. Banker, A. Sipahigil, M. H. Matheny, A. J. Keller, E. Kim, M. Mirhosseini, and O. Painter, Collapse and revival of an artificial atom coupled to a structured photonic reservoir, Phys. Rev. X 11, 041043 (2021).
  * [55] M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch, and A. A. Houck, Observation of a dissipative phase transition in a one-dimensional circuit QED lattice, Phys. Rev. X 7, 011016 (2017).
  * [56] D. C. McKay, R. Naik, P. Reinhold, L. S. Bishop, and D. I. Schuster, High-Contrast Qubit Interactions Using Multimode Cavity QED , Phys. Rev. Lett. 114, 080501 (2015).
  * [57] M. Scigliuzzo, G. Calajò, F. Ciccarello, D. P. Lozano, A. Bengtsson, P. Scarlino, A. Wallraff, D. Chang, P. Delsing, and S. Gasparinetti. Controlling Atom-Photon Bound States in an Array of Josephson-Junction Resonators, Phys. Rev. X 12, 031036 (2022).
  * [58] E. Kim, X. Zhang, V. S. Ferreira, J. Banker, J. K. Iverson, A. Sipahigil, M. Bello, A. González-Tudela, M. Mirhosseini, and O. Painter, Quantum Electrodynamics in a Topological Waveguide, Phys. Rev. X 11, 011015 (2021).
  * [59] A. Blais, S. M. Girvin, and W. D. Oliver, Quantum information processing and quantum optics with circuit quantum electrodynamics, Nat. Phys. 16, 247 (2020).
  * [60] J. Bourhill, N. Kostylev, M. Goryachev, D. L. Creedon, and M. E. Tobar, Ultrahigh cooperativity interactions between magnons and resonant photons in a YIG sphere, Phys. Rev. B 93, 144420 (2016).
  * [61] H. Huebl, C.W. Zollitsch, J. Lotze, F. Hocke, M. Greifenstein, A. Marx, R. Gross, and S.T. Goennenwein, High Cooperativity in Coupled Microwave Resonator Ferrimagnetic Insulator Hybrids, Phys. Rev. Lett. 111, 127003 (2013).
  * [62] X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Strongly Coupled Magnons and Cavity Microwave Photons, Phys. Rev. Lett. 113, 156401 (2014).
  * [63] Y. Tabuchi, S. Ishino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Hybridizing Ferromagnetic Magnons and Microwave Photons in the Quantum Limit, Phys. Rev. Lett. 113, 083603 (2014).
  * [64] P. Pirro, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, Advances in coherent magnonics, Nat. Rev. Mater. 6, 1114 (2021).

Supplementary Material for:
Spontaneous superradiant photon current

##  SM 1 Stability analysis

The solutions obtained by solving the mean-field equations (Eq. (3)-(6) in the
main text) can be stable or unstable. So it is necessary to check the
stability of these solutions. For a group of linear differential equations
which can be written as the matrix form $\partial V\left(t\right)/\partial
t=MV\left(t\right)$, if the real parts of all eigenvalues of the stability
matrix $M$ are not positive, then the solutions acquired by the steady-state
condition $\partial V\left(t\right)/\partial t=0$ are stable [1].
Specifically, one can rewrite the the mean-field solutions in terms of steady-
state solutions plus small fluctuations, namely, $\left\langle
S_{x}\right\rangle=N(X+\delta X)$, $\left\langle S_{y}\right\rangle=N(Y+\delta
Y)$, $\left\langle S_{z}\right\rangle=N(Z+\delta Z)$, $\left\langle
a_{n}\right\rangle=\sqrt{N}(\alpha_{n}+\delta\alpha_{n})$. Replacing the
previous order parameters in the mean-field equations by these ansatz, one can
obtain the following dynamical equations describing the fluctuations as

$\frac{d\delta\alpha_{n}}{dt}=-i\omega_{n}\delta\alpha_{n}-iJ\left(\delta\alpha_{n-1}+\delta\alpha_{n+1}\right)-i2g\delta
X-\kappa\delta\alpha_{n}\text{,}$ (SM1) $\frac{d\delta X}{dt}=-\Omega\delta
Y\text{,}$ (SM2) $\frac{d\delta Y}{dt}=\Omega\delta X-2g\delta
Z\sum_{n=1}^{N}\left(\alpha_{n}+\alpha_{n}^{\ast}\right)-2gZ\sum_{n=1}^{N}\left(\delta\alpha_{n}+\delta\alpha_{n}^{\ast}\right)\text{,}$
(SM3) $\frac{d\delta Z}{dt}=2g\delta
Y\sum_{n=1}^{N}\left(\alpha_{n}+\alpha_{n}^{\ast}\right)+2gY\sum_{n=1}^{N}\left(\delta\alpha_{n}+\delta\alpha_{n}^{\ast}\right)\text{.}$
(SM4)

Here, quantities $\alpha_{n}$, $\alpha_{n}^{\ast}$, $X$, $Y$, $Z$ take the
steady-state solutions which have been obtained in the main text. Note that
the parameters $X$, $Y$, $Z$ are not independent with each other. They are
restricted by the spin constraint $X^{2}+Y^{2}+Z^{2}=1/4$. Thus, the
fluctuation $\delta Z$ can be eliminated by

$\delta Z=\frac{X\delta X+Y\delta Y}{Z}\text{.}$ (SM5)

To proceed further we rewrite the quadratures $\delta
q_{n}=(\delta\alpha_{n}+\delta\alpha_{n}^{\ast})/\sqrt{2}$, $\delta
p_{n}=i(\delta\alpha_{n}^{\ast}-\delta\alpha_{n})/\sqrt{2}$. Then Eqs.
(SM1)-(SM4) can be rewritten in the matrix form

$\frac{\partial V\left(t\right)}{\partial t}=MV\left(t\right)\text{,}$ (SM6)

where $V\left(t\right)=[\delta q_{1},\delta p_{1},\delta q_{2},\delta
p_{2},\delta q_{3},\delta p_{3},\delta X,\delta Y]^{T}$ is the column vector
of the fluctuations. By numerically calculating the eigenvalues of matrix $M$
and checking the positivity of the real parts of these eigenvalues, one can
analyse the stability of the analytical solutions obtained in the main text.

##  SM 2 The derivation of current operators

By use of Heisenberg’s equation of motion about $a_{n}$ with cavity
dissipation rate $\kappa$, namely,

$\frac{da_{n}}{dt}=-\frac{i}{\hbar}\left[a_{n},H\right]-\kappa a_{n}\text{,}$
(SM7)

one can obtain the dynamical equations for the photon number operator
$\hat{n}_{n}=a_{n}^{{\dagger}}a_{n}$. That is,

$\displaystyle\frac{d\hat{n}_{n}}{dt}$
$\displaystyle=iJ\left(a_{n}a_{n-1}^{{\dagger}}-a_{n}^{{\dagger}}a_{n-1}\right)+iJ\left(a_{n}a_{n+1}^{{\dagger}}-a_{n}^{{\dagger}}a_{n+1}\right)$
$\displaystyle\quad+i\frac{g}{\sqrt{N}}\left[a_{n}S_{+}-a_{n}^{{\dagger}}S_{-}+a_{n}S_{-}-a_{n}^{{\dagger}}S_{+}\right]-2\kappa
a_{n}^{{\dagger}}a_{n}$
$\displaystyle\equiv\hat{I}_{n-1,n}+\hat{I}_{n+1,n}+\hat{I}_{s,n}+\hat{I}_{s,n}^{c}-I_{n,d}$
(SM8)

where current operators $\hat{I}_{n\pm 1,n}$, $\hat{I}_{s,n}$,
$\hat{I}_{s,n}^{c}$ and $I_{n,d}$ are defined by

$\hat{I}_{n\pm 1,n}=iJ\left(a_{n}a_{n\pm
1}^{{\dagger}}-a_{n}^{{\dagger}}a_{n\pm 1}\right)\text{,}$ (SM9)
$\hat{I}_{s,n}=i\frac{g}{\sqrt{N}}\left(a_{n}S_{+}-a_{n}^{{\dagger}}S_{-}\right)\text{,}$
(SM10)
$\hat{I}_{s,n}^{c}=i\frac{g}{\sqrt{N}}\left(a_{n}S_{-}-a_{n}^{{\dagger}}S_{+}\right)\text{,}$
(SM11) $I_{n,d}=2\kappa a_{n}^{{\dagger}}a_{n}\text{.}$ (SM12)

Here, $\hat{I}_{n\pm 1,n}$ is the current operator from cavity $n\pm 1$ to
cavity $n$ and $\hat{I}_{s,n}$ is the current operator from the spin system to
cavity $n$. The operator $\hat{I}_{s,n}^{c}$ arises from the counter-rotating-
wave terms in the Hamiltonian $H$. This existence of $\hat{I}_{s,n}^{c}$
clearly indicates that the number of photons and the number of emitters’
excitations can gain or lose simultaneously and hence violating the
conservation in the total number of photons and excitations. $I_{n,d}$
represents the operator associated with the current from cavity $n$ to the
environment by dissipation.

##  SM 3 Holstein-Primakoff representation and fluctuation equations of motion

With the Holstein-Primakoff transformation,
$S_{+}=b^{{\dagger}}\sqrt{N-b^{{\dagger}}b}$ and $S_{z}=-N/2+b^{{\dagger}}b$,
the Hamiltonian in Eq. (1) in the main text can be rewritten as

$\displaystyle H$
$\displaystyle=\sum_{n=1}^{N_{c}}\hbar\omega_{n}a_{n}^{{\dagger}}a_{n}+\sum_{n=1}^{N_{c}}\hbar
J\left(a_{n+1}^{{\dagger}}a_{n}+a_{n}^{{\dagger}}a_{n+1}\right)+\hbar\Omega\left(-\frac{N}{2}+b^{{\dagger}}b\right)$
$\displaystyle\quad+\sum_{n=1}^{N_{c}}\frac{\hbar
g_{n}}{\sqrt{N}}\left(b^{{\dagger}}\sqrt{N-b^{{\dagger}}b}+\sqrt{N-b^{{\dagger}}b}b\right)\left(a_{n}+a_{n}^{{\dagger}}\right)\text{.}$
(SM13)

We introduce the operators in terms of a mean-field part plus quantum
fluctuation operators with $a_{n}=\alpha_{n}\sqrt{N}+c_{n}$ and
$b=\beta\sqrt{N}+d$. Here, $\beta=(X-iY)/\sqrt{1/2-Z}$. $c_{n}$ is bosonic
annihilation operator of fluctuations for $n$th cavity and $d$ is the
emitter’s fluctuation operator. In the thermodynamics limit with a large-$N$
expansion, one can obtain the Hamiltonian for the fluctuations as

$\displaystyle H_{\text{fl}}$
$\displaystyle=\sum_{n=1}^{N_{c}}\hbar\left[\omega_{n}c_{n}^{{\dagger}}c_{n}+J\left(c_{n+1}^{{\dagger}}c_{n}+c_{n}^{{\dagger}}c_{n+1}\right)\right]+\hbar\tilde{\Omega}d^{{\dagger}}d$
$\displaystyle\quad+\sum_{n}\hbar\left[\Gamma_{1}\left(c_{n}^{{\dagger}}d+c_{n}d\right)+\Gamma_{3,n}d^{{\dagger}}d^{{\dagger}}+H.c\right]\text{,}$
(SM14)

where the coefficients $\tilde{\Omega}$, $\Gamma_{1}$ and $\Gamma_{3,n}$ are

$\tilde{\Omega}=\Omega-g\left[\sum_{n}\operatorname{Re}\left(\alpha_{n}\right)\right]\operatorname{Re}\left(\beta\right)\left[\frac{\left(4-3\left|\beta\right|^{2}\right)}{B_{1}^{3}}\right]\text{,}$
(SM15)
$\Gamma_{1}=g\frac{1-\left|\beta\right|^{2}-\operatorname{Re}\left(\beta\right)\beta^{\ast}}{B_{1}}\text{,}$
(SM16)
$\Gamma_{3,n}=-g\operatorname{Re}\left(\alpha_{n}\right)\frac{\left[\operatorname{Re}\left(\beta\right)\beta^{2}+2\beta\left(1-\left|\beta\right|^{2}\right)\right]}{2B_{1}^{3}}\text{.}$
(SM17)

Here, $B_{1}$ is defined as $B_{1}\equiv\sqrt{1-\left|\beta\right|^{2}}$.
$\operatorname{Re}\left(\alpha_{n}\right)$ and
$\operatorname{Re}\left(\beta\right)$ represent the real parts of $\alpha_{n}$
and $\beta$ respectively. By using the Heisenberg’s equation of motion with
the cavity dissipation, one can derive a set of closed differential equations
for fluctuation correlators, $\langle c_{m}^{{\dagger}}c_{n}\rangle$, $\langle
c_{m}c_{n}\rangle$, $\langle c_{m}d\rangle$, $\langle
c_{m}d^{{\dagger}}\rangle$, $\langle dd\rangle$ and $\langle
d^{{\dagger}}d\rangle$. Assuming that the system has reached its steady state,
one can easily obtain the photon number fluctuations $\langle
c_{1}^{{\dagger}}c_{1}\rangle$, $\langle c_{2}^{{\dagger}}c_{2}\rangle$ and
$\langle c_{3}^{{\dagger}}c_{3}\rangle$ by numerical calculations.

## References

  * [1] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series and Products (Academic, Orlando, 1980).


$\displaystyle v^{\top}\bar{A}_{i}v<0\quad\text{for all}\quad
v\in\mathbb{R}^{K_{i}}.$

Thus, the matrix $\bar{A}_{i}$ is Hurwitz, which further implies that $A_{i}$
is also Hurwitz. Following the steps similar to those in eq. 56 \- (59), we
can conclude that the matrix $G$ is also Hurwitz, which completes the proof.

We then verify Theorem 2. We proceed as follows,

$\displaystyle\left\|{\rm\Phi}\theta^{*}-G_{\pi}\right\|_{\mu_{\pi},\alpha}$
$\displaystyle\quad=\left\|{\rm\Gamma}_{{\rm\Phi},\mu_{\pi}}{\mathcal{T}}_{\pi}{\rm\Phi}\theta^{*}-{\mathcal{T}}_{\pi}G_{\pi}\right\|_{\mu_{\pi},\alpha}$
$\displaystyle\quad\leq\left\|{\rm\Gamma}_{{\rm\Phi},\mu_{\pi}}{\mathcal{T}}_{\pi}{\rm\Phi}\theta^{*}-{\rm\Gamma}_{{\rm\Phi},\mu_{\pi}}{\mathcal{T}}_{\pi}G_{\pi}\right\|_{\mu_{\pi},\alpha}+\left\|{\rm\Gamma}_{{\rm\Phi},\mu_{\pi}}\mathcal{T}_{\pi}G_{\pi}-{\mathcal{T}}_{\pi}G_{\pi}\right\|_{U_{\pi}}$
$\displaystyle\quad\leq\left\|{\rm\Gamma}_{{\rm\Phi},\mu_{\pi}}\hat{M}_{\pi}({\rm\Phi}\theta^{*}-G_{\pi})\right\|_{U_{\pi}}+\left\|{\rm\Gamma}_{{\rm\Phi},\mu_{\pi}}\mathcal{T}_{\pi}G_{\pi}-\mathcal{T}_{\pi}G_{\pi}\right\|_{\mu_{\pi},\alpha}$
$\displaystyle\quad\leq\left\|\hat{M}_{\pi}({\rm\Phi}\theta^{*}-G_{\pi})\right\|_{U_{\pi}}+\left\|{\rm\Gamma}_{{\rm\Phi},\mu_{\pi}}\mathcal{T}_{\pi}G_{\pi}-\mathcal{T}_{\pi}G_{\pi}\right\|_{\mu_{\pi},\alpha}$
$\displaystyle\quad\leq
C_{\zeta}\left\|{\rm\Phi}\theta^{*}-G_{\pi}\right\|_{U_{\pi}}+\left\|{\rm\Gamma}_{{\rm\Phi},\mu_{\pi}}\mathcal{T}_{\pi}G_{\pi}-\mathcal{T}_{\pi}G_{\pi}\right\|_{\mu_{\pi},\alpha},$
(109)

where the last inequality in eq. 109 can be obtained as follows. Using [13,
Theorem 1.3.22], we have

$\displaystyle\text{eig}(U^{-1}_{\pi,i}[\mathsf{P}_{\pi,i}^{\top}\otimes\mathsf{I}_{d_{i}}]U_{\pi,i})=\text{eig}([\mathsf{P}_{\pi,i}^{\top}\otimes\mathsf{I}_{d_{i}}]U_{\pi,i}U^{-1}_{\pi,i})=\text{eig}([\mathsf{P}_{\pi}^{\top}\otimes\mathsf{I}_{d_{i}}])=\text{eig}(\mathsf{P}_{\pi}^{\top})=\text{eig}(\mathsf{P}_{\pi}).$

Following the steps similar to those in eq. 56-(59), we can conclude that
$\text{eig}(\hat{M}_{\pi})=\text{eig}(\mathsf{P}_{\pi})$. Following the steps
similar to those for obtaining eq. 106. We have

$\displaystyle\left\|{\rm\Phi}\theta^{*}-G_{\pi}\right\|_{\mu_{\pi},\alpha}\leq\frac{1}{1-C_{\zeta}}\left\|{\rm\Gamma}_{{\rm\Phi},\mu_{\pi}}\mathcal{T}_{\pi}G_{\pi}-\mathcal{T}_{\pi}G_{\pi}\right\|_{\mu_{\pi},\alpha},$

where $0<C_{\zeta}<1$, which completes the proof.

##  Appendix H Proof of Theorem 3

We first define the matrix $B$ in the following way:

  * •

Forward GVF:
$B=\mathbb{E}_{\mathcal{D}\cdot\pi}[[\phi(s^{\prime},a^{\prime})\otimes\mathsf{I}_{d}]m(x)[\phi(s,a)\otimes\mathsf{I}_{d}]]$

  * •

Backward GVF:
$B=\mathbb{E}_{\mathcal{D}\cdot\pi}[[\phi(s,a)\otimes\mathsf{I}_{d}]m(x)[\phi(s^{\prime},a^{\prime})\otimes\mathsf{I}_{d}]]$.

We further define the following stochastic matrices in both the forward and
backward GVF evaluation settings. Recall that $(s_{t},a_{t})\sim D(\cdot)$,
$s^{\prime}_{t}\sim\mathsf{P}(\cdot|s_{t},a_{t})$ and
$a^{\prime}_{t}\sim\pi(\cdot|s^{\prime}_{t}$.

  * •

Forward GVF:

$\displaystyle A_{t}$
$\displaystyle=[\phi(s_{t},a_{t})\otimes\mathsf{I}_{d}](m(x_{t})[\phi(s^{\prime}_{t},a^{\prime}_{t})\otimes\mathsf{I}_{d}]^{\top}-[\phi(s_{t},a_{t})\otimes\mathsf{I}_{d}]^{\top}),$
$\displaystyle B_{t}$
$\displaystyle=[\phi(s^{\prime}_{t},a^{\prime}_{t})\otimes\mathsf{I}_{d}]m(x_{t})[\phi(s_{t},a_{t})\otimes\mathsf{I}_{d}],$
$\displaystyle C_{t}$
$\displaystyle=(\phi(s_{t},a_{t})\phi(s_{t},a_{t})^{\top})\otimes\mathsf{I}_{d},$
$\displaystyle b_{t}$
$\displaystyle=[\phi(s_{t},a_{t})\otimes\mathsf{I}_{d}]C(x_{t}).$ (110)

  * •

Backward GVF:

$\displaystyle A_{t}$
$\displaystyle=[\phi(s^{\prime}_{t},a^{\prime}_{t})\otimes\mathsf{I}_{d}](m(x_{t})[\phi(s_{t},a_{t})\otimes\mathsf{I}_{d}]^{\top}-[\phi(s^{\prime}_{t},a^{\prime}_{t})\otimes\mathsf{I}_{d}]^{\top}),$
$\displaystyle B_{t}$
$\displaystyle=[\phi(s_{t},a_{t})\otimes\mathsf{I}_{d}]m(x_{t})[\phi(s^{\prime}_{t},a^{\prime}_{t})\otimes\mathsf{I}_{d}],$
$\displaystyle C_{t}$
$\displaystyle=(\phi(s^{\prime}_{t},a^{\prime}_{t})\phi(s_{t},a_{t})^{\top})\otimes\mathsf{I}_{d},$
$\displaystyle b_{t}$
$\displaystyle=[\phi(s^{\prime}_{t},a^{\prime}_{t})\otimes\mathsf{I}_{d}]C(x_{t}).$
(111)

Recall the matrices $A$ and $C$ defined in Appendix C. For a constant $\xi>0$,
we define

$\displaystyle H_{t}$ $\displaystyle=\left[\begin{array}[]{cc}A_{t}&B_{t}\\\
\xi A_{t}&\xi C_{t}\end{array}\right],\qquad
h_{t}=\left[\begin{array}[]{c}b_{t}\\\ 0\end{array}\right].$ (116)

and

$\displaystyle H$ $\displaystyle=\left[\begin{array}[]{cc}A&B\\\ \xi A&\xi
C\end{array}\right],\qquad h=\left[\begin{array}[]{c}b\\\
0\end{array}\right],$ (121)

where $A=\mathbb{E}[A_{t}]$, $B=\mathbb{E}[B_{t}]$, $C=\mathbb{E}[C_{t}]$ and
$b=\mathbb{E}[b_{t}]$.

For the matrix $H_{t}$, we have the following holds

$\displaystyle\left\|H_{t}\right\|_{F}^{2}$
$\displaystyle=(1+\xi^{2})\left\|A_{t}\right\|_{F}^{2}+\left\|B_{t}\right\|_{F}^{2}+\xi^{2}\left\|C_{t}\right\|_{F}^{2}$
$\displaystyle\leq(1+\xi^{2})[d^{2}C^{2}_{\phi}(C_{m}+1)]^{2}+d^{2}C^{2}_{\phi}C^{2}_{m}+\xi^{2}C^{4}_{\phi}d^{2}.$
(122)

which implies that $\left\|H_{t}\right\|_{F}\leq C_{H}$, where

$\displaystyle
C_{H}=\sqrt{(1+\xi^{2})[d^{2}C^{2}_{\phi}(C_{m}+1)]^{2}+d^{2}C^{2}_{\phi}C^{2}_{m}+\xi^{2}C^{4}_{\phi}d^{2}}.$

For the vector $h_{t}$, we can obtain $\left\|h_{t}\right\|_{2}\leq
C_{h}=dC_{\phi}R_{C}$ by following the steps similar to those for obtaining
eq. 122.

The update in Algorithm 2 can be rewritten as

$\displaystyle
v_{t+1}={\rm\Gamma}_{R_{v}}\left(v_{t}+\alpha_{t}(H_{t}v_{t}+h_{t})\right),$
(123)

where $v_{t}=[\theta_{t}^{\top},w_{t}^{\top}]^{\top}$, and
$R_{v}=R_{\theta}\times\mathbb{R}^{Kd_{g}\times 1}$. Following the proof
similar to those in [22, Section 5.3.3, Theorem 3], we can show that the
matrix $H$ is Hurwitz under 1 and 2 with an appropriately chosen
$\xi>\max\\{0,-\text{eig}_{\min}(C^{-1}[(A+A^{\top})/2])\\}$.

We define the following optimal point
$v^{*}=[\bar{\theta}^{*\top},w^{*\top}]^{\top}$ for the linear SA defined in
eq. 123

$\displaystyle\langle\varphi(v^{*}),v-v^{*}\rangle\leq 0,\quad\forall v\in
R_{v},$

where $\varphi(v)=Hv+b$. We also define $C_{v}=\left\|v^{*}\right\|_{2}$. It
can be checked that there exist a positive constant $\lambda^{\prime}_{G}$
such that

$\displaystyle\langle\varphi(v^{*})-\varphi(v),v^{*}-v\rangle\leq-\lambda^{\prime}_{G}\left\|v-v^{*}\right\|_{2}^{2}.$
(124)

We further define $\varphi(x_{t},v)=H_{t}v+h_{t}$.

Following the steps similar to those for proving Lemma 7 and Lemma 8, we can
obtain the following two lemmas.

######  Lemma 11.

Given a sample $(s_{t},a_{t},B_{t},s^{\prime}_{t})\sim\mathcal{D}_{d}$ and
$a^{\prime}_{t}\sim\pi(\cdot|s^{\prime}_{t})$ and any $v\in R_{v}$, we have
the following holds

$\displaystyle\left\|\varphi(x_{t},v)-\varphi(v)\right\|_{2}^{2}\leq
16C^{2}_{H}\left\|\kappa-\kappa^{*}\right\|_{2}^{2}+16C^{2}hPC^{2}_{v}+8C^{2}_{h}.$

###### Proof.

Based on the definition of $\varphi(x_{t},v)$ and $\varphi(v)$, we can obtain
the following

$\displaystyle\left\|\varphi(x_{t},v)-\varphi(v)\right\|_{2}^{2}$
$\displaystyle\leq
2\left\|(H_{t}-H)v\right\|_{2}^{2}+2\left\|h_{t}-h\right\|_{2}^{2}$
$\displaystyle\leq
4\left\|(H_{t}-H)(v-v^{*})\right\|_{2}^{2}+4\left\|(H_{t}-H)v^{*}\right\|_{2}^{2}+2\left\|h_{t}-h\right\|_{2}^{2}$
$\displaystyle\leq
16C^{2}_{H}\left\|\kappa-\kappa^{*}\right\|_{2}^{2}+16C^{2}_{h}C^{2}_{v}+8C^{2}_{h}.$

∎

######  Lemma 12.

Consider the population GTD update $\varphi(v)=Hv+b$. We have

$\displaystyle\langle-\varphi(v),v-v^{*}\rangle\geq\lambda^{\prime}_{G}\left\|v-v^{*}\right\|_{2}^{2},\quad\forall
v\in R_{v}.$

We also have the following "three-point lemma" holds for the GTD update.

######  Lemma 13.

Consider the update of $w_{t}$ and $\theta_{t}$ in Algorithm 2. For all $v\in
R_{v}$, we have the following holds

$\displaystyle-\alpha_{t}\langle\varphi(x_{t},v_{t}),v_{t+1}-v\rangle+\frac{1}{2}\left\|v_{t+1}-v_{t}\right\|_{2}^{2}\leq\frac{1}{2}\left\|v_{t}-v\right\|_{2}^{2}-\frac{1}{2}\left\|v_{t+1}-v\right\|_{2}^{2}.$
(125)

Using Lemma 13 and following the steps similar to those from eq. 88 to eq. 91,
we can obtain

$\displaystyle\frac{1}{2}\left\|v_{t}-v\right\|_{2}^{2}-\frac{1}{2}\left\|v_{t+1}-v\right\|_{2}^{2}$
$\displaystyle\quad\geq-\alpha_{t}\langle\varphi(v_{t+1}),v_{t+1}-v\rangle-\alpha_{t}\langle\varphi(x_{t},v_{t})-\varphi(v_{t}),v_{t}-v\rangle-\alpha^{2}_{t}C^{2}_{H}\left\|v_{t+1}-v\right\|_{2}^{2}$
$\displaystyle\quad\quad-\alpha^{2}_{t}\left\|\varphi(x_{t},\kappa_{t})-\varphi(\kappa_{t})\right\|_{2}^{2}.$
(126)

Taking expectation on both sides of eq. 126, letting $v=v^{*}$, and using the
fact that
$-\langle\varphi(v_{t+1}),v_{t+1}-v^{*}\rangle\leq\lambda^{\prime}_{G}\left\|v_{t+1}-v^{*}\right\|_{2}$
yield

$\displaystyle(1+2\alpha_{t}\lambda^{\prime}_{G}-2\alpha^{2}_{t}C^{2}_{H})\mathbb{E}[\left\|v_{t+1}-v^{*}\right\|_{2}^{2}]$
$\displaystyle\leq\mathbb{E}[\left\|v_{t}-v^{*}\right\|_{2}^{2}]+2\alpha^{2}_{t}\mathbb{E}[\left\|\varphi(x_{t},\kappa_{t})-\varphi(\kappa_{t})\right\|_{2}^{2}]$
$\displaystyle\leq(1+32\alpha^{2}_{t}C^{2}_{H})\mathbb{E}[\left\|v_{t}-v^{*}\right\|_{2}^{2}]++32C^{2}_{h}C^{2}_{v}+16C^{2}_{h},$
(127)

where the second inequality follows from Lemma 12. Multiplying both sides of
eq. 127 by $o_{t}$ and summing over iterations $t=0,\cdots,T-1$ yield

$\displaystyle\sum_{t=0}^{T-1}a^{\prime\prime}_{t}\mathbb{E}[\left\|v_{t+1}-v^{*}\right\|_{2}^{2}]\leq\sum_{t=0}^{T-1}b^{\prime\prime}_{t}\mathbb{E}[\left\|v_{t}-v^{*}\right\|_{2}^{2}]+c^{\prime\prime},$
(128)

where

$\displaystyle a^{\prime\prime}_{t}$
$\displaystyle=(1+2\alpha_{t}\lambda^{\prime}_{G}-2\alpha^{2}_{t}C^{2}_{H})o_{t},$
$\displaystyle b^{\prime\prime}_{t}$
$\displaystyle=(1+32\alpha^{2}_{t}C^{2}_{H})o_{t},$ $\displaystyle
c^{\prime\prime}$
$\displaystyle=(32C^{2}hPC^{2}_{v}+16C^{2}_{h})\sum_{t=0}^{T-1}\alpha^{2}_{t}o_{t}.$

Now we let

$\displaystyle o_{t}$ $\displaystyle=(t+t_{2})(t+t_{2}+1),$
$\displaystyle\alpha_{t}$
$\displaystyle=\frac{4}{\lambda^{\prime}_{G}(t+t_{2}-1)},$ $\displaystyle
t_{2}$ $\displaystyle=\frac{34C^{2}_{H}}{\lambda^{2}_{J}}+1.$

Then, we can obtain the following

$\displaystyle a^{\prime\prime}_{t}-b^{\prime\prime}_{t+1}$
$\displaystyle=(1+2\alpha_{t}\lambda^{\prime}_{G}-2\alpha^{2}_{t}C^{2}_{H})s_{t}-(1+32\alpha_{t+1}^{2}C^{2}_{H})o_{t+1}$
$\displaystyle\geq(1+2\alpha_{t}\lambda^{\prime}_{G})o_{t}-(1+2\alpha^{2}_{t}C^{2}_{H}+32\alpha_{t+1}^{2}C^{2}_{H})o_{t+1}$
$\displaystyle\geq(1+2\alpha_{t}\lambda^{\prime}_{G})o_{t}-(1+34\alpha_{t}^{2}C^{2}_{H})o_{t+1}$
$\displaystyle\geq(1+2\alpha_{t}\lambda^{\prime}_{G})o_{t}-(1+\alpha_{t}\lambda^{\prime}_{G})o_{t+1}$
$\displaystyle\geq(t+t_{2}+1)\frac{(t+t_{2}+3)(t+t_{2})-(t+t_{2}+2)^{2}}{t+t_{2}-1}$
$\displaystyle\geq 0,$

where the second inequality follows from the fact that
$\alpha_{t}\leq\frac{\lambda^{\prime}_{G}}{34C^{2}_{H}}$.

Applying the above property to eq. 128 yields

$\displaystyle
a^{\prime\prime}_{T-1}t\mathbb{E}[\left\|v_{T}-v^{*}\right\|_{2}^{2}]\leq
b^{\prime\prime}_{0}\left\|v_{0}-v^{*}\right\|_{2}^{2}+c^{\prime\prime},$

which implies

$\displaystyle\mathbb{E}[\left\|v_{T}-v^{*}\right\|_{2}^{2}]$
$\displaystyle\leq\frac{b^{\prime\prime}_{0}\left\|v_{0}-v^{*}\right\|_{2}^{2}}{a^{\prime\prime}_{T-1}}+\frac{c^{\prime\prime}}{a^{\prime\prime}_{T-1}}$
$\displaystyle\leq\frac{(1+16\alpha^{2}_{0}C^{2}_{H})(t_{2}+1)^{2}\left\|v_{0}-v^{*}\right\|_{2}^{2}}{(T+t_{2}-1)(T+t_{2})}+\frac{128C^{2}_{h}C^{2}_{v}+64C^{2}_{h}}{(T+t_{2})\lambda^{2}_{J}}.$

Using the fact
$\left\|\theta_{T}-\bar{\theta}^{*}\right\|_{F}^{2}\leq\left\|v_{T}-v^{*}\right\|_{2}^{2}$,
we have

$\displaystyle\mathbb{E}[\left\|\theta_{T}-\bar{\theta}^{*}\right\|_{F}^{2}]\leq\frac{(1+16\alpha^{2}_{0}C^{2}_{H})(t_{2}+1)^{2}\left\|v_{0}-v^{*}\right\|_{2}^{2}}{(T+t_{2}-1)(T+t_{2})}+\frac{128C^{2}_{h}C^{2}_{v}+64C^{2}_{h}}{(T+t_{2})\lambda^{2}_{J}},$

which completes the proof.