#  Self-organized cavity bosons beyond the adiabatic elimination approximation

Giuliano Orso<EMAIL_ADDRESS>Université Paris Cité, Laboratoire
Matéeriaux et Phéenomènes Quantiques (MPQ), CNRS, F-75013, Paris, France
Jakub Zakrzewski<EMAIL_ADDRESS>Institute of Theoretical Physics,
Jagiellonian University in Krakow, ul. Łojasiewicza 11, 30-348 Kraków, Poland
Mark Kac Complex Systems Research Center, Jagiellonian University in Krakow,
Łojasiewicza 11, 30-348 Kraków, Poland  Piotr Deuar<EMAIL_ADDRESS>Institute of Physics, Polish Academy of Sciences, Aleja Lotników 32/46, 02-668
Warsaw, Poland

###### Abstract

The long-time behavior of weakly interacting bosons moving in a two-
dimensional optical lattice and coupled to a lossy cavity is investigated
numerically in the regime of high particle filling. The truncated Wigner
representation allows us to take into full account the dynamics of the cavity
mode, quantum fluctuations, and self-organization of individual runs. We
observe metastability at very long times and superfluid quasi-long range
order, in sharp contrast with the true long range order found in the ground
state of the approximate Bose-Hubbard model with extended interactions,
obtained by adiabatically eliminating the cavity field. As the strength of the
light-matter coupling increases, the system first becomes supersolid at the
Dicke superradiant transition and then turns into a charge-density wave via
the Berezinskii-Kosterlitz-Thouless mechanism. The two phase transitions are
characterized via an accurate finite-size scaling.

today

_Introduction._ Recent experiments explore the regime of strong coupling of
quantum matter and quantum light on a variety of platforms, ranging from
electronic systems coupled to THz cavities [1] to atomic Bose gases coupled to
a dissipative cavity [2, 3], where self-organized density patterns have been
observed [4, 5] by tuning the system across the Dicke superradiant transition.
Depending on the strength of the atom-cavity coupling, the crystal order can
either coexist with the superfluid (SF) phase in the celebrated supersolid
(SS) phase or replace the SF order by a charge density wave (CDW) insulator
(see Fig. 1).

Figure 1: Long-time density (top) and phase (middle) distributions of bosons
of a given TW trajectory on a $64\times 64$ lattice at $tJ=10^{4}$. The
columns correspond to three increasing values of the coupling Rabi frequency
$\Omega=1.5J$ (left), $4J$ (center), and $10J$ (right), where the system is in
the SF, SS, and CDW phases, respectively. The square symbols in the middle row
signal the presence of a vortex (black) or an antivortex (red) in the
plaquette. Other parameters are $\delta=2J,\kappa=0.5J$, $U=0.2J$ and
$n_{b}=5$. Bottom: distribution of the density wave order parameter $\Delta$
over an ensemble of 800 trajectories on $16\times 16$ lattices.

Most theoretical works have concentrated on the so-called “bad cavity” regime:
If the detuning of the cavity $\delta$ and the photon leak rate $\kappa$ are
large compared to the characteristic energy scales of the gas, namely the
tunneling rate $J$ and the interaction strength $U$, then the cavity mode can
be adiabatically eliminated resulting in an extended Bose-Hubbard (EBH) model,
characterized by infinite-range instantaneous interactions. The ground state
phase diagram of the EBH model obtained using Gutzwiller mean field theories
[6, 7, 8], quantum Monte Carlo [9] and slave-boson [10] techniques, predicts
that the SF and SS phases are characterized by a true long-range order and
that the SS-CDW phase transition is of second order [6, 7, 8, 9]. Additional
interesting phases are found for one dimensional (1D) EBH cavity model [11,
12].

How does the physics change when we go beyond this adiabatic approximation? A
natural question is whether the true long-time state is the same as the one
provided by the simplified bad cavity model. A pioneering step was taken in
[13, 14] using tensor network simulations for 1D system. Follow on work
suggested [15] that cavity bosons in the steady state can be described as an
equilibrium thermal state of the EBH model, whose effective temperature is
obtained self-consistently, by including perturbatively the effect of quantum
fluctuations around the mean field solution. However, this requires computing
the retarded dynamical susceptibility in the EBH model as a function of the
frequency, which is challenging for the 2D systems relevant for ongoing
experiments [4, 5].

Importantly, the atom-cavity system differs markedly from photon-polariton
systems [16] in that here only the single cavity mode is dissipative, while
the overwhelming majority of modes are conservative. Therefore the atoms
experience dissipative effects only indirectly, which might lead to a very
slow relaxation towards the steady state, and make it out-of-reach for 1D
quasi-exact simulations [13, 14].

In this Letter we study large 2D systems through the use of the truncated
Wigner (TW) approach, which keeps full account of the cavity degree of freedom
with no a priori assumption on the structure of the combined mixed state, and
includes quantum fluctuations as well as nonlinear effects seeded by them
[17]. It also allows us to reach long times and compute with a good accuracy
correlation functions and distribution functions that are crucial to unveil
the order of phase transitions and to identify the associated critical points.

_Model and Method._ The Hamiltonian describing the system can be written as
$\widehat{H}=\widehat{H}_{b}+\widehat{H}_{c}+\widehat{H}_{bc}$, where

$\widehat{H}_{b}=-J\sum_{\langle
i,j\rangle}\widehat{b}^{\dagger}_{i}\widehat{b}_{j}+\frac{U}{2}\sum_{j}\widehat{n}_{j}(\widehat{n}_{j}-1)$
(1)

is the Bose-Hubbard model describing the matter component, with
$\widehat{b}^{\dagger}_{j},\widehat{b}_{j}$ the local field operators and
$\widehat{n}_{j}=\widehat{b}^{\dagger}_{j}\widehat{b}_{j}$ the local density
operator at each lattice site $j$. We assume that the cavity contains only one
relevant mode of energy $\omega_{c}$, so that its bare Hamiltonian in the
rotating wave approximation is
$\widehat{H}_{c}=\delta\,\widehat{a}^{\dagger}\widehat{a}$ with
$\widehat{a}^{\dagger},\widehat{a}$ being the photon field operators and
$\delta=\omega_{c}-\omega_{p}$ the detuning between the cavity and the pump
frequencies (hereafter we use $\hbar=1$ units). The coupling between the
bosons and the cavity is described by the Hamiltonian

$\widehat{H}_{bc}=-\frac{\Omega}{\sqrt{N}}\left(\widehat{a}+\widehat{a}^{\dagger}\right)\widehat{\Delta},$
(2)

where $\Omega$ is the (re-scaled) Rabi frequency,
$N=\langle\sum_{j}\widehat{n}_{j}\rangle$ is the conserved total number of
bosons, and $\widehat{\Delta}=\sum_{j}f_{j}\widehat{n}_{j}$ is a collective
operator in which $f_{j}$ is the overlap of the cavity mode with the
wavefunction amplitude in the $j$th lattice site taking, typically, a
staggered “checkerboard” form [4]; we use $f_{j}=(-1)^{j_{x}+j_{y}}$ here. The
single shot values of $\widehat{\Delta}$ tend to $\pm N$ when atoms are
completely arranged in even/odd lattice sites, so it can be considered as a
collective magnetization order parameter for even/odd arrangement. The
dynamics of the hybrid system is described by the master equation for the
density matrix $\widehat{\rho}$ in the Markovian limit with photon leak rate
$\kappa$:

$\frac{d\widehat{\rho}}{dt}=-i\left[\widehat{H},\widehat{\rho}\right]+\kappa\left[2\widehat{a}\widehat{\rho}\widehat{a}^{\dagger}-\widehat{a}^{\dagger}\widehat{a}\widehat{\rho}-\widehat{\rho}\widehat{a}^{\dagger}\widehat{a}\right].$
(3)

We address the evolution and long-time behavior of the system via the TW
approach [18, 19, 17, 20, 21, 22, 23]. The field operators
$\widehat{a},\widehat{b}_{j}$ become represented by an ensemble of complex
variables $\alpha,\beta_{j}$ and the evolution equation (3) by stochastic
differential equations [24]

$\displaystyle id\alpha$ $\displaystyle\\!\\!=$
$\displaystyle\\!\\!\Big{[}(\delta-i\kappa)\alpha-\frac{\Omega}{\sqrt{N}}\sum_{j}f_{j}\left(|\beta_{j}|^{2}-\tfrac{1}{2}\right)\Big{]}dt+i\sqrt{\kappa}dW_{c}$
$\displaystyle id\beta_{j}$ $\displaystyle=$
$\displaystyle\Big{[}U\left(|\beta_{j}|^{2}-1\right)\beta_{j}-\frac{\Omega}{\sqrt{N}}(\alpha+\alpha^{*})f_{j}\beta_{j}$
(4) $\displaystyle-J\sum_{r\in X(j)}\beta_{r}\Big{]}dt$

where the set $X(j)$ contains the indices of all the neighboring sites of the
site of index $j$ (there are four nearest neighbors in our square $L\times L$
lattice with periodic boundary conditions) and $dW_{c}$ is a complex-valued
Wiener noise term with $\langle|dW_{c}|^{2}\rangle=dt$. Single trajectories of
(4) with a given noise realization correspond closely to experimental runs
[25, 26], while quantum expectation values are obtained by averaging the Weyl
symbols over many trajectories (numbering 240-2400 here). For example, the
expectation value of the local density is
$\langle\widehat{n}_{j}\rangle=\langle|\beta_{j}|^{2}-\tfrac{1}{2}\rangle$.
Remarkably the approach scales only linearly with the number of sites, giving
access to large systems (the $64\times 64$ lattices treated here are still far
away from any significant limitations [23]), while retaining the majority of
quantum fluctuation effects. The method is found to be accurate for weak
onsite interactions and relatively large boson filling [19, 27, 25], as
recently demonstrated for the driven-dissipative Bose Hubbard model [23].
Advantages of the TW include the possibility to reach very long times because
the method is not inherently limited by the entanglement growth of tensor
network methods [14, 13, 15]. TW also gives an easy access to full
distributions of single experimental shots (as shown in Fig. 1) and
correlations Fig. 2). For details of our implementation see [24].

Figure 2: Superfluid phase correlation function (5), versus evolution time for
a $16\times 16$ lattice averaged over 800 trajectories; Curves red to blue
show four successively increasing distances $r=1$ to $4$. The coupling is
$\Omega=2.4J$ and other parameters are as in Fig. 2. Note the metastable
plateau at times $tJ\sim{\mathcal{O}}(100)$.

In the following we explore a regime in which the energy scales of the cavity
and of the matter components are of the same order so that fluctuation,
retardation or dissipation effects can play a significant role. We take
$\kappa=J/2$, $\delta=2J$, $U=0.2J$ and $n_{b}=N/L^{2}=5$. The latter is
roughly twice the average boson filling used in the ETH experiment [4], but we
choose it to increase the accuracy of the TW method. Evolution times are
$tJ=10^{4}$ (i.e. about 200 times longer than in [14]) unless otherwise noted,
starting from an initial state with atoms in a condensate and vacuum in the
cavity photon mode. Deviations from the adiabatic elimination estimate are
indeed observed very clearly for these parameters[24], being
${\mathcal{O}}(1)$ for $\widehat{a}$.

_Metastability and self-organization._ Metastable dynamics is observed for
many $\Omega$ values, particularly strongly in the range
$2\lesssim\Omega/J\lesssim 5$. A typical example is shown in Fig. 2 where the
time dependence of the first order (phase) correlation function

$g(r)=\frac{1}{N}\sum_{j}\textrm{Re}(\langle\widehat{b}^{\dagger}_{j}\widehat{b}_{j+r}\rangle)$
(5)

is plotted for several distances $r$ ( the index $j+r$ is modulo $L$ to
conform to the periodic boundary conditions). Here, a quasi stationary state
appears at the times $tJ\sim{\mathcal{O}}(20)$ reached by previous studies,
while true convergence occurs only at $tJ\approx 10^{4}$.

As the coupling $\Omega$ is increased, three phases are encompassed, with
typical observables in single shots shown in Fig. 1. For small values of
$\Omega$, the density (top) shows strong site-to-site fluctuations, while the
phase (lower) is essentially uniform, corresponding to a SF phase. As $\Omega$
increases, the density distribution becomes a checkerboard, signalling the
appearance of a density wave order and antibunching correlations between
neighboring sites. The actual order is self-organized spontaneously, as
evidenced by the $\Delta/N$ distribution peaked at $\pm 1$, corresponding to
the occupation of even or odd sites. At the same time, the phase distribution
shows sizeable fluctuations, along with the formation of scattered vortex-
antivortex pairs, that we identify by calculating the phase changes
accumulated in a plaquette of nearest neighbors. Their binding implies that
the system possesses both SF and charge order, the hallmark of the SS phase. A
further increase of the Rabi frequency (right) leads to a more pronounced
staggering of the density that is accompanied by a massive proliferation of
vortices and antivortices, which act to destroy the SF order. We identify this
normal phase as a CDW insulator and argue that the SS-CDW phase transition
should obey the BKT scenario, according to which superfluidity in 2D systems
with $U(1)$ symmetry is lost due to the unbinding of the topological
excitations.

Figure 3:  Long-time behavior of the 1st order correlation (5) at $tJ=10^{4}$,
versus the lattice distance $r=|j_{x}-i_{x}|$ along one direction, averaged
over 240 trajectories for a $64\times 64$ lattice. More precisely,
${\widetilde{r}}=(L/\pi)\sin(\pi r/L)$ is the chord distance. Notice the
change from logarithmic to linear scale in the horizontal axis of panel d),
emphasizing the transition from power-law decay QLRO (a-c), to exponential
decay (d) in the proximity of the BKT critical point (see text). We extract
the exponents $\eta$ from fits (red solid lines).

_Quasi long-range order (QLRO)._ To characterize the atom superfluidity in the
long-time regime, we inspect the spatial dependence of the first order
correlation (5) in Fig. 3. In the SF and SS phases it displays a power-law
decay of the form $g(r)=(a+(-1)^{r}b\,)/\tilde{r}^{\eta}$, where $b$ is
nonzero in the presence of solid order and ${\widetilde{r}}=(L/\pi)\sin(\pi
r/L)$ is the “chord” distance, which accounts to leading order for the
periodic boundary conditions [28]. For $\Omega=8J$ (c), where the system is
close to the SS-CDW transition (see below), a fit to the data, with $a,b,\eta$
as free fitting parameters, yields $\eta=0.261(5)$, which is slightly larger
than the critical anomalous exponent $\eta_{c}^{XY}=1/4$, holding for the XY
universality class at equilibrium. We stress that the cavity-induced QLRO
observed here contrasts with the true long-range SF order observed in previous
ground state calculations of the 2D EBH model, reflecting the shortcomings of
the cavity adiabatic elimination. Interestingly, QLRO has also been reported
[16] for 2D microcavity polaritons coherently driven into the optical
parametric oscillator regime, with $\eta$ taking values even larger than one.
In the CDW phase (d) the function $g(r)$ decays exponentially, consistently
with the BKT scenario.

Figure 4: Main plot: static structure factor $S_{\pi,\pi}(L)$ versus $\Omega$
across the superfluid-supersolid phase transition, showing finite size
scaling. Black points represent the $L\to\infty$ thermodynamic limit
$A(\Omega)$, according to (6). The dashed black line is a linear fit of
$A(\Omega)$ in the above threshold phase, giving the Dicke transition point as
$\Omega_{s}=2.265(5)$. The inset shows the Binder cumulant $U_{L}$ for the
“magnetization” $\Delta/N$ across the transition. Here the crossing of
$U_{L}=0$ occurs at $\Omega_{s}=2.265(2)$. Data at $tJ=10^{4}$, 800
trajectories.

_Identification of critical points._ We next employ the finite-size scaling
approach to pinpoint the critical values $\Omega_{s}$ and $\Omega_{c}$ where
the two phase transitions occur and characterize their nature. At
$\Omega=\Omega_{s}$ the system undergoes the Dicke superradiant transition,
where the cavity mode becomes macroscopically occupied and, concomitantly, the
solid order appears [29, 30]. We characterize the latter via the static
structure factor
$S_{q}=\sum_{ij}e^{iq(i-j)}(\langle\widehat{n}_{i}\widehat{n}_{j}\rangle-\langle\widehat{n}_{i}\rangle\langle\widehat{n}_{j}\rangle)/L^{4}$,
with $q=(\pi,\pi)$. Fig. 4 displays the $\Omega$-dependence of $S_{\pi,\pi}$
for increasing system sizes. The cavity occupation shows the same features, as
shown in [24]. In the SF phase the structure factor scales to zero as
$1/L^{2}$, while in the SS phase it converges to a finite value. From a
finite-size scaling analysis according to [24]

$S_{\pi,\pi}(\Omega,L)=A(\Omega)+\frac{C(\Omega)}{L^{2}}$ (6)

and then following the very visibly linear trend of
$A(\Omega)\propto\Omega-\Omega_{s}$ we obtain $\Omega_{s}=2.265(5)$. The
linear in $\Omega$ trend close to the critical point, also observed for 1D
systems [14], confirms that the Dicke superradiant transition can be described
with Landau (mean-field) theory [31].

The TW method also allows us to calculate the Binder cumulant
$U_{L}(\Delta)=1-\langle\Delta^{4}\rangle/3\langle\Delta^{2}\rangle^{2}$ for
the effective magnetization $\Delta$ – something that is ruled out in the mean
field or bad cavity approximations and therefore has not been investigated to
date. It is shown in Fig. 4(inset) and provides the more accurate value
$\Omega_{s}=2.265(2)$ for the critical point at which the jump from $U_{L}=0$
to $U_{L}=2/3$ appears in the $L\to\infty$ limit. Before the transition there
is an interesting regime in which $U_{L}<0$, that appears to stem from the
appearance of fatter tails in the probability distribution of $\Delta$ as a
precursor of the phase transition [24].

Figure 5: Main panel: correlation ratio $r_{L}=g(L/2)/g(L/4)$ at $tJ=10^{4}$,
$2400$ trajectories, as a function of $\Omega$ for different system sizes
(data symbols). The solid lines are fits to the data curves with Padé
approximants $p/q$, where $p$ and $q$ have polynomial degrees $n_{p}=3$ and
$n_{q}=4$. Inset: Identification of the BKT critical point based on Eq. (8)
with $R=0.8$.

We now turn to the BKT transition. For periodic boundary conditions, a useful
indicator is provided by the correlation ratio [32] $r_{L}=g(L/2)/g(L/4)$, In
Fig. 5 we plot the long-time behavior of $r_{L}$ against $\Omega$ for
increasing system sizes. For $\Omega\lesssim 8J$, the correlation ratio is
essentially size independent, reflecting the QLRO. In contrast, for
$\Omega\gtrsim 8J$ the different curves spray out, consistent with the fact
that the correlation ratio must scale to zero in the CDW phase, suggesting
that $\Omega_{c}\approx 8J$. For a more accurate estimate, we use the fact
that the correlation ratio is described by a single parameter scaling
function, $r_{L}=f(L/\xi)$, where $\xi$ is the correlation length in the
normal phase $(\Omega>\Omega_{c})$

$\xi=\exp\left(\frac{c}{\sqrt{\Omega/\Omega_{c}-1}}\right),$ (7)

with $c>0$ being a model-dependent constant. Following [32], we define
$\Omega_{KT}(L)$ as the value of $\Omega$ at which $r_{L}=0.8=R$. From Eq. (7)
it follows that in the scaling regime $\Omega_{KT}(L)$ satisfies
$L\exp(-c/\sqrt{u})=1/b$, where $u=\Omega_{KT}(L)/\Omega_{c}-1$ and $b$ is
related to $R$ via $f(1/b)=R$, so that [32]

$\Omega_{KT}(L)=\Omega_{c}+\frac{c^{2}\Omega_{c}}{\ln^{2}bL}.$ (8)

We use Eq. (8) to fit our data, with $\Omega_{c},c$ and $b$ being fitting
parameters. The fit is displayed in the inset of Fig. 5, confirming that the
linear dependence is well satisfied, with the intercept
$\Omega_{c}/J=8.03(4)$. We have verified that using a different choice of $R$,
with $0.6\lesssim R<0.8$, yields $\Omega_{c}$ consistent within error with the
above.

_Conclusions._ We have characterized the long-time phases of 2D boson-cavity
lattices in a regime where the adiabatic elimination approximation is not
justified. We have unveiled that the SF and SS phases exhibit QLRO, as
compared to the true long-range order suggested by the EBH model at zero
temperature. The QLRO in turn disappears at the SS-CDW phase boundary
according to the BKT scenario. The precise position of the critical points of
the Dicke and of the BKT phase transitions have been identified. Moreover we
demonstrated that the TW method allows one to investigate the long-time
properties of large 2D boson-cavity lattices, without being hindered by the
bad cavity approximation or by specific assumptions on the nature of the
atomic mixed state. We confirm that, as suspected from the small number of
dissipative degrees of freedom, the system is prone to a very slow evolution
that can be missed by time-limited approaches, which are only able to reach
the initial plateau. Such long timescales are nevertheless very relevant for
ongoing experiments [4] that typically deal with times $tJ\sim 10^{4}-10^{5}$.

We are grateful to I. Carusotto, C. Ciuti, F. Hebert, S.B. Jager, King Lun Ng,
and T. Roscilde for helpful explanations and discussions. This research was
co-financed by the Polish National Agency for Academic Exchange, grant No.
BPN/BFR/2022/1/00027 and by the HPC Polonium grant No. 49317RK. This work was
granted access to the HPC resources of TGCC under the allocation AD010513635
and AD010513635R1 supplied by GENCI (Grand Equipement National de Calcul
Intensif). This research was also funded by National Science Centre (Poland)
under the OPUS call within the WEAVE programme 2021/43/I/ST3/01142 (J.Z.).

## References

  * Forn-Díaz  _et al._ [2019] 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).
  * Ritsch  _et al._ [2013] H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger, Cold atoms in cavity-generated dynamical optical potentials, Rev. Mod. Phys. 85, 553 (2013).
  * Farokh Mivehvar and Ritsch [2021] T. D. Farokh Mivehvar, Francesco Piazza and H. Ritsch, Cavity qed with quantum gases: new paradigms in many-body physics, Advances in Physics 70, 1 (2021), https://doi.org/10.1080/00018732.2021.1969727 .
  * Landig  _et al._ [2016] R. Landig, L. Hruby, N. Dogra, M. Landini, R. Mottl, T. Donner, and T. Esslinger, Quantum phases from competing short- and long-range interactions in an optical lattice, Nature 532, 476 (2016).
  * Dreon  _et al._ [2022] 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 608, 494 (2022).
  * Dogra  _et al._ [2016] N. Dogra, F. Brennecke, S. D. Huber, and T. Donner, Phase transitions in a Bose-Hubbard model with cavity-mediated global-range interactions, Phys. Rev. A 94, 023632 (2016).
  * Sundar and Mueller [2016] B. Sundar and E. J. Mueller, Lattice bosons with infinite-range checkerboard interactions, Phys. Rev. A 94, 033631 (2016).
  * Himbert  _et al._ [2019] L. Himbert, C. Cormick, R. Kraus, S. Sharma, and G. Morigi, Mean-field phase diagram of the extended Bose-Hubbard model of many-body cavity quantum electrodynamics, Phys. Rev. A 99, 043633 (2019).
  * Flottat  _et al._ [2017] T. Flottat, L. d. F. de Parny, F. Hébert, V. G. Rousseau, and G. G. Batrouni, Phase diagram of bosons in a two-dimensional optical lattice with infinite-range cavity-mediated interactions, Phys. Rev. B 95, 144501 (2017).
  * Sharma  _et al._ [2022] S. Sharma, S. B. Jäger, R. Kraus, T. Roscilde, and G. Morigi, Quantum critical behavior of entanglement in lattice bosons with cavity-mediated long-range interactions, Phys. Rev. Lett. 129, 143001 (2022).
  * Chanda  _et al._ [2021] T. Chanda, R. Kraus, G. Morigi, and J. Zakrzewski, Self-organized topological insulator due to cavity-mediated correlated tunneling, Quantum 5, 501 (2021).
  * Chanda  _et al._ [2022] T. Chanda, R. Kraus, J. Zakrzewski, and G. Morigi, Bond order via cavity-mediated interactions, Phys. Rev. B 106, 075137 (2022).
  * Halati  _et al._ [2020a] C.-M. Halati, A. Sheikhan, and C. Kollath, Theoretical methods to treat a single dissipative bosonic mode coupled globally to an interacting many-body system, Phys. Rev. Research 2, 043255 (2020a).
  * Halati  _et al._ [2020b] C.-M. Halati, A. Sheikhan, H. Ritsch, and C. Kollath, Numerically exact treatment of many-body self-organization in a cavity, Phys. Rev. Lett. 125, 093604 (2020b).
  * Bezvershenko  _et al._ [2021] A. V. Bezvershenko, C.-M. Halati, A. Sheikhan, C. Kollath, and A. Rosch, Dicke transition in open many-body systems determined by fluctuation effects, Phys. Rev. Lett. 127, 173606 (2021).
  * Dagvadorj  _et al._ [2015] G. Dagvadorj, J. M. Fellows, S. Matyjaśkiewicz, F. M. Marchetti, I. Carusotto, and M. H. Szymańska, Nonequilibrium phase transition in a two-dimensional driven open quantum system, Phys. Rev. X 5, 041028 (2015).
  * Norrie  _et al._ [2005] A. A. Norrie, R. J. Ballagh, and C. W. Gardiner, Quantum turbulence in condensate collisions: An application of the classical field method, Phys. Rev. Lett. 94, 040401 (2005).
  * Steel  _et al._ [1998] M. J. Steel, M. K. Olsen, L. I. Plimak, P. D. Drummond, S. M. Tan, M. J. Collett, D. F. Walls, and R. Graham, Dynamical quantum noise in trapped Bose-Einstein condensates, Phys. Rev. A 58, 4824 (1998).
  * Sinatra  _et al._ [2002] A. Sinatra, C. Lobo, and Y. Castin, The truncated Wigner method for Bose-condensed gases: limits of validity and applications, Journal of Physics B: Atomic, Molecular and Optical Physics 35, 3599 (2002).
  * Wouters and Savona [2009] M. Wouters and V. Savona, Stochastic classical field model for polariton condensates, Phys. Rev. B 79, 165302 (2009).
  * Chiocchetta and Carusotto [2014] A. Chiocchetta and I. Carusotto, Quantum langevin model for nonequilibrium condensation, Phys. Rev. A 90, 023633 (2014).
  * Ruostekoski and Martin [2013] J. Ruostekoski and A. D. Martin, The truncated Wigner method for Bose gases, in  _Quantum Gases_ (Imperial College Press, 2013) Chap. 13, pp. 203–214.
  * Deuar  _et al._ [2021] P. Deuar, A. Ferrier, M. Matuszewski, G. Orso, and M. H. Szymańska, Fully quantum scalable description of driven-dissipative lattice models, PRX Quantum 2, 010319 (2021).
  * [24] See Supplemental Material (below) for additional information concerning the TW approach, the accuracy of the bad cavity approximation, and details on the phase transitions studied.
  * Blakie  _et al._ [2008] P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh, and C. W. Gardiner, Dynamics and statistical mechanics of ultra-cold Bose gases using c-field techniques, Advances in Physics 57, 363 (2008).
  * Lewis-Swan  _et al._ [2016] R. J. Lewis-Swan, M. K. Olsen, and K. V. Kheruntsyan, Approximate particle number distribution from direct stochastic sampling of the Wigner function, Phys. Rev. A 94, 033814 (2016).
  * Norrie  _et al._ [2006] A. A. Norrie, R. J. Ballagh, and C. W. Gardiner, Quantum turbulence and correlations in Bose-Einstein condensate collisions, Phys. Rev. A 73, 043617 (2006).
  * Gerster  _et al._ [2016] M. Gerster, M. Rizzi, F. Tschirsich, P. Silvi, R. Fazio, and S. Montangero, Superfluid density and quasi-long-range order in the one-dimensional disordered Bose–Hubbard model, New Journal of Physics 18, 015015 (2016).
  * Dicke [1954] R. H. Dicke, Coherence in spontaneous radiation processes, Phys. Rev. 93, 99 (1954).
  * Hepp and Lieb [1973] K. Hepp and E. H. Lieb, On the superradiant phase transition for molecules in a quantized radiation field: the dicke maser model, Annals of Physics 76, 360 (1973).
  * Kirton  _et al._ [2019] P. Kirton, M. M. Roses, J. Keeling, and E. G. Dalla Torre, Introduction to the dicke model: From equilibrium to nonequilibrium, and vice versa, Advanced Quantum Technologies 2, 1800043 (2019).
  * Tomita and Okabe [2002] Y. Tomita and Y. Okabe, Finite-size scaling of correlation ratio and generalized scheme for the probability-changing cluster algorithm, Phys. Rev. B 66, 180401 (2002).
  * Drummond and Corney [2001] P. D. Drummond and J. F. Corney, Quantum noise in optical fibers. I. Stochastic equations, J. Opt. Soc. Am. B 18, 139 (2001).
  * Cahill and Glauber [1969] K. E. Cahill and R. J. Glauber, Density operators and quasiprobability distributions, Phys. Rev. 177, 1882 (1969).
  * Deuar [2021] P. Deuar, Multi-time correlations in the positive-P, Q, and doubled phase-space representations, Quantum 5, 455 (2021).
  * Plimak  _et al._ [2001] L. I. Plimak, M. K. Olsen, M. Fleischhauer, and M. J. Collett, Beyond the Fokker-Planck equation: Stochastic simulation of complete Wigner representation for the optical parametric oscillator, Europhysics Letters (EPL) 56, 372 (2001).
  * Polkovnikov [2003] A. Polkovnikov, Quantum corrections to the dynamics of interacting bosons: Beyond the truncated Wigner approximation, Phys. Rev. A 68, 053604 (2003).
  * Drummond [2014] P. D. Drummond, Fundamentals of higher order stochastic equations, J. Phys. A 47, 335001 (2014).
  * Olsen and Bradley [2009] M. Olsen and A. Bradley, Numerical representation of quantum states in the positive-P and Wigner representations, Optics Communications 282, 3924 (2009).
  * Kinsler and Drummond [1991] P. Kinsler and P. D. Drummond, Quantum dynamics of the parametric oscillator, Phys. Rev. A 43, 6194 (1991).
  * Drobný  _et al._ [1997] G. Drobný, A. Bandilla, and I. Jex, Quantum description of nonlinearly interacting oscillators via classical trajectories, Phys. Rev. A 55, 78 (1997).
  * Deuar and Drummond [2007] P. Deuar and P. D. Drummond, Correlations in a BEC collision: First-principles quantum dynamics with 150 000 atoms, Phys. Rev. Lett. 98, 120402 (2007).
  * Carrasquilla and Rigol [2012] J. Carrasquilla and M. Rigol, Superfluid to normal phase transition in strongly correlated bosons in two and three dimensions, Phys. Rev. A 86, 043629 (2012).

Supplementary material for:
Self-organized cavity bosons beyond the adiabatic elimination approximation

Giuliano Orso1, Jakub Zakrzewski2,3, and Piotr Deuar4

1Université Paris Cité, Laboratoire Matéeriaux et Phéenomènes Quantiques
(MPQ), CNRS, F-75013, Paris, France
2 Institute of Theoretical Physics, Jagiellonian University in Krakow, ul.
Łojasiewicza 11, 30-348 Kraków, Poland
3 Mark Kac Complex Systems Research Center, Jagiellonian University in Krakow,
Łojasiewicza 11, 30-348 Kraków, Poland
4 Institute of Physics, Polish Academy of Sciences, Al. Lotników 32/46, 02-668
Warsaw, Poland

###### Contents

  1. S1 Truncated Wigner Methodology
    1. A Overview of the method and its derivation
    2. B Initial conditions and observables
    3. C Simulation details
    4. D Assessment of TW accuracy
  2. S2 Accuracy of the bad cavity approximation
  3. S3 Details of the Dicke SF / SS transition
  4. S4 Condensate fraction near the BKT phase transition

Note: Citation numbers in square brackets refer to the bibliography in the
main paper.

##  S1 Truncated Wigner Methodology

###  A Overview of the method and its derivation

The “truncated” Wigner method is an adaptation of the widely known Wigner
phase-space quasiprobablity distribution $W(x,p)$ for representing the state
of a quantum mode, made to allow its application to very large many-mode
systems. It is quite widely used e.g. in ultracold atom and quantum optics
systems [18, 33, 19, 17, 22] and also has been applied to polariton and open
systems of the Bose-Hubbard type [20, 21, 23].

For $M$ sites, instead of studying the quasiprobability distribution
$W(\vec{x},\vec{p})$ as a distribution in a $2M$-dimensional space, values of
the amplitude vector $\vec{\alpha}=\vec{x}+i\vec{p}$ are sampled according to
the distribution $W(\vec{\alpha})$. Here amplitudes $\alpha_{j}=x_{j}+ip_{j}$,
corresponding to space and momentum operators
$\widehat{x}_{j}=(\widehat{a}_{j}+\widehat{a}^{\dagger}_{j})/2$ and
$\widehat{p}_{j}=(\widehat{a}_{j}-\widehat{a}^{\dagger}_{j})/2i$ in an
abstract “harmonic oscillator space” for each site $j$ are similar (though not
equal) to coherent state amplitudes. These amplitudes then obey stochastic
differential equations with noise amplitudes appropriate to represent the
dissipative quantum mechanical evolution. Appropriate averages over the
ensemble of $S$ samples or “trajectories” correspond then better and better to
the quantum state as $S\to\infty$.

More precisely, the density matrix $\widehat{\rho}$ can be expanded in terms
of the Wigner quasiprobaility distribution as

$\widehat{\rho}=\int
d\vec{\alpha}\,W(\vec{\alpha})\widehat{\Lambda}_{W}(\vec{\alpha}),$ (S1)

where $\widehat{\Lambda}_{W}$ is an operator kernel [34, 35] dependent on the
set of complex phase space amplitudes $\vec{\alpha}=\\{\alpha,\beta_{j}\\}$,
corresponding to the 2nd quantized operators $\widehat{a},\widehat{b}_{j}$ for
the photon mode and boson sites, respectively. With the correct Wigner
representation kernel, operators obey the differential identities:

$\displaystyle\widehat{a}\widehat{\Lambda}_{W}$ $\displaystyle=$
$\displaystyle\left[\alpha-\frac{1}{2}\frac{\partial}{\partial\alpha^{*}}\right]\widehat{\Lambda}_{W}$
(S2a) $\displaystyle\widehat{a}^{\dagger}\widehat{\Lambda}_{W}$
$\displaystyle=$
$\displaystyle\left[\alpha^{*}+\frac{1}{2}\frac{\partial}{\partial\alpha}\right]\widehat{\Lambda}_{W}$
(S2b) $\displaystyle\widehat{b}_{j}\widehat{\Lambda}_{W}$ $\displaystyle=$
$\displaystyle\left[\beta_{j}-\frac{1}{2}\frac{\partial}{\partial\beta_{j}^{*}}\right]\widehat{\Lambda}_{W}$
(S2c) $\displaystyle\widehat{b}^{\dagger}_{j}\widehat{\Lambda}_{W}$
$\displaystyle=$
$\displaystyle\left[\beta_{j}^{*}+\frac{1}{2}\frac{\partial}{\partial\beta_{j}}\right]\widehat{\Lambda}_{W},$
(S2d)

and their adjoints. With their help, using standard phase-space representation
procedure, we can rewrite the master equation (3) as an integral equation

$\displaystyle\int d\vec{\alpha}\,\frac{dW}{dt}\widehat{\Lambda}_{W}$
$\displaystyle=\int
d\vec{\alpha}\,W\left[\sum_{k}A_{k}\frac{\partial}{\partial\alpha_{k}}+\sum_{kk^{\prime}}\frac{D_{kk^{\prime}}}{2}\frac{\partial^{2}}{\partial\alpha_{k}\partial\alpha_{k^{\prime}}}+\dots\right]\widehat{\Lambda}_{W}$

where $\alpha_{k}$ label all variables in $\vec{\alpha}$ and “$\dots$”
indicate higher order partial derivatives (up to 3rd order in our case, as it
turns out). Integrating by parts and assuming vanishing boundary terms at
$|\alpha|,|\beta|_{j}\to\infty$, this gives $\int
d\vec{\alpha}\,\widehat{\Lambda}_{W}\frac{dW}{dt}=\int\widehat{\Lambda}_{W}\left[\circ\circ\circ\right]W$,
which has at least one solution in the form of a Fokker-Planck equation (FPE)
with additional third order terms:

$\displaystyle\frac{dW}{dt}$ $\displaystyle=$
$\displaystyle\left[\circ\circ\circ\right]=\Bigg{[}\sum_{k}\frac{\partial}{\partial\alpha_{k}}(-A_{k})+\sum_{kk^{\prime}}\frac{\partial^{2}}{\partial\alpha_{k}\partial\alpha_{k^{\prime}}}\frac{D_{kk^{\prime}}}{2}$
(S4)
$\displaystyle\qquad-\sum_{kk^{\prime}k"}\frac{\partial^{3}}{\partial\alpha_{k}\partial\alpha_{k^{\prime}}\partial\alpha_{k"}}T_{kk^{\prime}k"}\Bigg{]}W.$

In our particular case, in hopefully self-evident notation, the nonzero
elements are

$\displaystyle A_{\alpha}$ $\displaystyle=$
$\displaystyle-i\delta\alpha-\kappa\alpha+i\frac{\Omega}{\sqrt{N}}\sum_{j}f_{j}\left(|\beta_{j}|^{2}-\tfrac{1}{2}\right)$
$\displaystyle A_{\beta_{j}}$ $\displaystyle=$ $\displaystyle-
iU(|\beta_{j}|^{2}-1)\beta_{j}+i\frac{\Omega}{\sqrt{N}}(\alpha+\alpha^{*})f_{j}\beta_{j}$
$\displaystyle+iJ\sum_{r\in X(j)}\beta_{r}$ $\displaystyle
D_{\alpha\alpha^{*}}$ $\displaystyle=$ $\displaystyle
D_{\alpha^{*}\alpha}=\kappa$ (S5) $\displaystyle
T_{\beta_{j}\beta_{j}\beta^{*}_{j}}$ $\displaystyle=$ $\displaystyle
i\frac{U\beta_{j}}{2}$ $\displaystyle T_{\beta_{j}\beta_{j}^{*}\alpha}$
$\displaystyle=$ $\displaystyle-i\frac{\Omega f_{j}}{4\sqrt{N}},$

as well as $A_{\alpha_{k}^{*}}=(A_{\alpha_{k}})^{*}$ and
$T_{\alpha_{k}^{*}\alpha_{k}^{\prime*}\alpha_{k}"^{*}}=(T_{\alpha_{k}\alpha_{k}^{\prime}\alpha_{k}"})^{*}$.

The $A$ and $D$ terms lead to drift and diffusion processes for samples
$\vec{\alpha}$ of the distribution. Although some notable attempts at
including third order terms like $T$ numerically have been made [36, 37, 38],
there are no known well scalable implementations. The standard truncated
Wigner procedure is to now remove these offending third order terms but remain
in a regime where their effect is negligible. Studies have shown that this is
justified in ultracold atom systems when some conditions are met, especially
that the mean occupation per mode is greater than
${\mathcal{O}}(\tfrac{1}{2})$ (single mode occupations need not) [27]. In our
case the mean occupation is $n_{b}=5$. TW accuracy is looked at in more detail
in S1 D.

The Fokker-Planck equation (S4), after removal of the $T_{kk^{\prime}k"}$,
corresponds to the following stochastic Langevin equations in the Ito
calculus:

$d\alpha_{k}=A_{\alpha_{k}}dt+\sum_{l}B_{\alpha_{k},l}\,dW_{l}$ (S6)

where the diffusion is decomposed as $D=BB^{T}$ using noise matrices $B$ with
matrix elements $B_{\alpha_{k},l}$. The $dW_{l}$ are real Wiener increments
with $\langle dW_{l}\rangle=0,\langle
dW_{l}dW_{l^{\prime}}\rangle=\delta_{ll^{\prime}}$. In our case the noise
matrix is very simple, and one can take

$B_{\alpha,1}=B_{\alpha^{*},1}=\sqrt{\frac{\kappa}{2}},B_{\alpha,2}=-B_{\alpha^{*},2}=i\sqrt{\frac{\kappa}{2}}.$
(S7)

as the only nonzero elements. Substitution of (A) and (S7) into the general
form (S6) yields the truncated Wigner evolution equations for the system, (4).

While the evolution equations (4) bear a lot of resemblance to the c-field or
“Gross-Pitaevskii” formulation, obtained by removing the noises $dW_{l}$, the
interpretation of the amplitudes is different. The Wigner distribution of a
vacuum state is a Gaussian
$W(\alpha)=\frac{2}{\pi}\exp\left[-2|\alpha|^{2}\right]$ with mean
$\langle|\alpha|^{2}\rangle=\tfrac{1}{2}$, and of a coherent state
$\alpha_{0}$ the same Gaussian centered at $\alpha_{0}$. Therefore, it
includes “virtual vacuum noise” that encodes the Heisenberg-limited minimum
uncertainty between non-commuting operators such as $\widehat{a}$ and
$\widehat{a}^{\dagger}$ via the randomness of the ensemble.

###  B Initial conditions and observables

Our initial conditions are vacuum in the cavity photon mode and a condensate
with mean occupation $n_{b}=5$ per site in the atomic lattice. i.e. [39]

$\alpha_{j}(0)=\frac{\xi}{\sqrt{2}}\qquad;\qquad\beta_{j}(0)=\sqrt{n_{b}}+\frac{\chi_{j}}{\sqrt{2}}.$
(S8)

Here $\xi$ and $\chi_{j}$ are independent complex random variables, whose real
and imaginary parts are normal random variables with zero mean and standard
deviation equal to $1/\sqrt{2}$, so that
$\langle|\chi_{j}|^{2}\rangle=\langle|\xi|^{2}\rangle=1$ and other 1st and 2nd
order stochastic averages are zero. We use a coherent initial state because
the Wigner distribution does not describe a Fock state conveniently (though it
is possible to do if needed, see [39]).

In concert with the virtual vacuum fluctuations, physical observables must
also be calculated in a more nuanced way than simply replacing
$\widehat{b}_{j}\to\beta_{j}$ as one would do with mean field treatments.
Expectation values of an observable $\widehat{O}$ are found by taking the
average of its corresponding Weyl symbol (which can be obtained from
application of (S2a) and (S1) on Tr[$\widehat{O}\widehat{\rho}$]). In
particular,

$\displaystyle\langle\widehat{b}^{\dagger}_{j}\widehat{b}_{i}\rangle$
$\displaystyle=$
$\displaystyle\langle\beta_{j}^{*}\beta_{i}\rangle-\frac{\delta_{ij}}{2},$
(S9a) $\displaystyle\langle\widehat{b}_{j}\rangle$ $\displaystyle=$
$\displaystyle\langle\beta_{j}\rangle$ (S9b)
$\displaystyle\langle\widehat{n}_{j}\widehat{n_{i}}\rangle$ $\displaystyle=$
$\displaystyle\left\langle\left(|\beta_{j}|^{2}-\frac{1}{2}\right)\left(|\beta_{i}|^{2}-\frac{1}{2}\right)\right\rangle-\frac{\delta_{ij}}{4},\quad$
(S9c)

and analogous expressions for the photon mode. (S9) are then used to calculate
observables such as occupations, phases, $g(r)$, $S_{\pi,\pi}$, and $\Delta$.

Single trajectories of a TW simulation are known to correspond to a
significant degree to experimental measurements, on the same basis that Gross-
Pitaevskii simulations are a good predictor of single experimental runs. Only
here, quantum fluctuation effects are additionally taken into account [25].
The correspondence is generally up to the level of minimum Heisenberg-limited
uncertainty which is reflected by the virtual noise in the initial conditions.
That is – uncertainty on the order of $\tfrac{1}{2}$ particle on one site,
$1/\sqrt{2n_{j}}$ on the local phase at site $j$, or $\sqrt{L^{2}}/2$ for
density observables on the whole system. As a reflection of this, for example,
single trajectories do not capture $\pm 1$ modulation of particle number in
e.g. squeezed states [26], though ensemble averages remain fine. For
$\Delta/N$ such uncertainty/washing out comes out of the order of
$1/(2Ln_{b})$. As we can see, sizeable site occupations such as $n_{j}\sim
n_{b}$ are expected to lead to a qualitatively correct representation of
single runs under most conditions.

On the other hand non-observable quantities such as entanglement entropy or
purity of the state are not accessible, except through emulating quantum
tomography measurements.

###  C Simulation details

As per the main text, physical parameters are $\kappa=J/2$, $\delta=2J$,
$U=0.2J$ and $n_{b}=N/L^{2}=5$, the mean number of particles per site. In the
SF and CDW phases, local mean density can approach $n_{j}\approx 10$, with
neighbors $n_{j+1}\approx{\mathcal{O}}(1)$. The atoms are on a square 2d
lattice with periodic boundary conditions, of dimensions $L\times L$. Total
evolution times (except for Fig. 2) are $tJ=10^{4}$, which is also a typical
waiting time used in experiments [4]. For expectation values, 240-2400
trajectories were averaged, depending on the case and the signal to noise
ratio.

We have used two kinds of integration scheme, both of which give mutually
consistent results.

(i) the Runge-Kutta 4 method with the time step $dt$ chosen such that the
numerical relative error in the conservation of the total boson number $N$ is
sufficiently small, $|N-\sum_{j}\langle\hat{n}_{j}\rangle|/N<10^{-4}$ over the
evolution time; For the range of parameters explored this yields a time step
varying between $10^{-4}J^{-1}$ and $10^{-3}J^{-1}$.

(ii) a symmetric split step method in which the terms of $d\beta_{j}$ in $U$
and $\Omega$ are performed exponentially in real space over $dt$,

$\beta_{j}\to\beta_{j}\exp\left\\{-idt\left[U(|\beta_{j}|^{2}-1)-2\frac{\Omega}{\sqrt{N}}{\rm
Re}[\alpha]f_{j}\right]\right\\}$ (S10)

while the tunneling term in $J$ is performed in two half-steps of $dt/2$ in
k-space (reached by a Fourier transform of $\beta_{j}$ to
${\widetilde{\beta}}_{\vec{k}}$), using the exponential factors

${\widetilde{\beta}}_{\vec{k}}\to{\widetilde{\beta}}_{\vec{k}}\exp\left[-i\left(\frac{dt}{2}\right)2J\sum_{l=1,2}\left(1-\cos\left(\frac{2\pi
k_{l}}{L}\right)\,\right).\right]$ (S11)

Here $\vec{k}=[k_{1},k_{2}]$ are the lattice indices in k-space, in dimensions
$l=1,2$. Double half-steps help to bring the integration error in local
$n_{j}$ down to $\propto dt^{3}t$. This split-step approach conserves total
atom numbers exactly, and requires a $dt\sim 10^{-3}J^{-1}$.

Figure S1: Nearest neighbor density-density correlation
$g^{(2)}(r=1)=\langle\widehat{b}^{\dagger}_{0}\widehat{b}^{\dagger}_{1}\widehat{b}_{0}\widehat{b}_{1}\rangle/n_{b}^{2}$
in the long-time state at $tJ=10^{4}$ as a function of the cavity coupling
$\Omega$ for a $16\times 16$ lattice.  Figure S2: Distribution of site
occupations $n_{j}=|\beta_{j}|^{2}-\tfrac{1}{2}$ at $tJ=10^{4}$, $L=16$ for
several characteristic $\Omega$ values. 800 trajectories.

###  D Assessment of TW accuracy

The conditions under which the third order terms in the Fokker-Planck Equation
(S4) can be omitted have been the focus of quite a number of investigations,
e.g. [19, 27, 40, 41, 36, 37, 42], especially in the context of the simulation
of quantum optical and ultracold atom systems. A fair amount is known about
what conditions correspond to accurate simulations, and what symptoms mark a
simulation as suspicious. Several practical conditions for a reliable
simulation have emerged in the literature: [25, 19, 27].

  1. 1.

The mean occupation per mode should be greater than $\approx\tfrac{1}{2}$ [19,
27]. However, notably, single mode occupations need not meet the criterion
provided the overall system does [27].

  2. 2.

Local observables such as $n_{j}$ and the density-density correlations
$g^{(2)}(r)=\langle\widehat{n}_{0}(\widehat{n}_{r}-\delta_{0r})\rangle/(n_{0}n_{r})$
should remain within physical limits such as $n_{j}>0$ and $g^{(2)}(r)>0$.

  3. 3.

For Hamiltonian systems that redistribute mode occupations under internal
Hamiltonian dynamics, problems consisting of incorrect redistribution have
been seen when the maximum excitation energies allowed exceed several times
$k_{B}T_{\rm eff}$. Here $T_{\rm eff}$ is the final effective temperature
reached at long times [19]. This condition is correlated with condition 1, in
that allowing very many high energy modes above thermal energies, so that they
are poorly occupied, decreases the mean occupation per mode.

How do our system and simulations compare to the above criteria?

(1) ✓The mean occupation is $n_{b}=5$, meeting the average occupation
criterion well.

(2) With regard to local observables,

  * •

✓mean site occupations $n_{j}$ are $n_{b}=5$ so there are no problems with
negative densities on average. Note that some $n_{j}<0$ in particular
trajectories are expected, since they are essential to describe a vacuum.

  * •

✓With regard to density-density correlation, the lowest values in our system
occur for nearest neighbors ($r=1$). Fig. S1 shows the dependence of this
correlation $g^{(2)}(r=1)$ on $\Omega$. There is no indication of nonphysical
$g^{(2)}<0$ occurring here, though the exact values of $g^{(2)}(r)$ values
very close to zero should be taken with a grain of salt.

  * •

The distribution of site occupations $n_{j}$ is shown in Fig. S2 in the main
regimes studied. In all cases the overwhelming majority of modes have
occupations above 1, and most sites in the low density part of the
checkerboard pattern also keep to $n_{j}>1$. Single shot values are not
expected to conform to $n_{j}>0$ anyhow, since single trajectory values of
such low occupations do not correspond precisely to measurements and e.g. a
vacuum has half shots with $n_{j}<0$ and is represented fine. ✓

Overall we see no indications of problems with observables keeping to physical
limits.

(3) While temperature is not easy to quantify for our system, we can consider
single particle energy, and whether there are many unoccupied high energy
modes. As best seen perhaps by comparing (S11) to free space, which would have
the exponential factor exp[$-i(dt/2)(\hbar k)^{2}/2m$], the basic single-
particle energies are $\epsilon_{\vec{k}}=2J(2-\cos(2\pi k_{1}/L)-\cos(2\pi
k_{2}/L)\,)$. The charge density wave, which is significantly occupied for all
regimes except for the SF, occupies in fact the highest energy kinetic
excitations at $k_{1}=k_{2}=L/2$. Therefore the situation of many unoccupied
higher energy modes does not take place there. Moreover, our system is not
strictly Hamiltonian, therefore one does not expect redistribution problems in
modes whose dissipation timescale is shorter than the rather long
thermalization time. ✓

In conclusion, our TW simulations meet the conditions considered necessary for
good accuracy, and provide for the inclusion of spontaneous and quantum
fluctuating effects at least to leading orders, and in a non-perturbative way.
There are also no direct indications of any appreciable inaccuracies, though
some care should be taken with results that depend primarily on low occupied
modes.

Figure S3: Data from 800 single trajectories of the $L=16$ full model:
Im[$\alpha$] vs. $\Delta$ (points) at $tJ=10^{4}$, compared to the adiabatic
elimination estimate (S12) (dashed line).

Figure S4: Comparison of the actual cavity mode amplitude (turqoise) to that
assumed by the adiabatic elimination (“bad cavity”) approximation as given by
(S12) (dark blue). At $t=10^{4}$ $L=16$ for several characteristic $\Omega$ in
the superfluid (a) and supersolid (b-c) phases. Panels (d-e) show the
difference $\delta a=\alpha-a_{\rm est}$ in individual trajectories.

##  S2 Accuracy of the bad cavity approximation

The bad cavity (adiabatic elimination) approximation assumes that the
dissipative timescale is short enough for the lossy mode (here $\widehat{a}$)
to collapse onto its steady state faster than other system quantities evolve.
Therefore it takes the steady state value of its evolution equation when other
quantities (here $\beta_{j},\Delta$) are assumed constant. Taking the steady
state limit of (4) for the cavity mode by imposing $dW=d\alpha=0$, one obtains
$(\delta-i\kappa)\alpha-\frac{\Omega}{\sqrt{N}}\sum_{j}f_{j}\left(|\beta_{j}|^{2}-\frac{1}{2}\right)=0$.
Solving for $\alpha$ gives the estimate $\alpha\to a_{\rm est}$ where

$a_{\rm
est}=\frac{\Omega}{\sqrt{N}}\left(\frac{\Delta}{\delta-i\kappa}\right).$ (S12)

Substitution of $\alpha=a_{\rm est}$ into the evolution equations (4) for
$\beta_{j}$ gives the TW equations for the EBH model. The EBH model in its
operator form is obtained by making the substitution

$\widehat{a}\to\frac{\Omega}{\sqrt{N}}\left(\frac{\widehat{\Delta}}{\delta-i\kappa}\right)$
(S13)

in the Hamiltonian $\hat{H}$ and by discarding the Linblad term in the
evolution equation (3). This gives

$\widehat{H}_{bc}+\widehat{H}_{c}=-\frac{\Omega^{2}}{N}\,\frac{\delta}{\delta^{2}+\kappa^{2}}\,\widehat{\Delta}^{2},$
(S14)

which is a long-distance interaction since
$\widehat{\Delta}^{2}=\sum_{ij}f_{i}f_{j}\widehat{n}_{i}\widehat{n}_{j}$.

Fig. S3 shows the scatter of the cavity mode $\alpha$ in single realizations
with respect to the one-to-one relationship (S12) indicated in the bad cavity
limit. Fig. S4 shows the relationship in more detail – for each realization it
compares the approximation (S12) made when using the atomic $\Delta$ with the
actual cavity mode amplitude. A random valued discrepancy is evident in both
regimes: superfluid (panels a,d) and supersolid (panels b-c,e). It is
particularly strong in the super-fluid, while in the density correlated phases
the adiabatic approximation $a_{\rm est}$ retains the overall sign and
approximate phase seen in the full model, but collapses the photon amplitude
to a single line. The error in $\alpha$ is of the order of ${\mathcal{O}}(1)$
for our parameters, with a more or less random phase. It corresponds to an
error of order $\sqrt{N_{a}}$ in the cavity photon number and of order
$1/\sqrt{N_{a}}$ in the cavitye phase.

##  S3 Details of the Dicke SF / SS transition

First, Fig. S5 shows the very close relationship between the static structure
factor and the cavity mode occupation. In particular
$S_{\pi,\pi}/(n_{b}L)^{2}\propto N_{a}/\Omega^{2}$ is quite closely followed
in the whole range of $\Omega$ we have studied.

Figure S5: Behavior of the static structure factor $S_{\pi,\pi}$ and cavity
occupation $N_{a}$ and the relationship between them. Shown for the system
with the usual parameters at $tJ=10^{4}$, on a $16\times 16$ lattice, averaged
over 800 trajectories.

Figure. S6 shows the convergence of finite-size values of $S_{\pi,\pi}$ as $L$
grows. The plot is tailored to a $1/L^{2}$ scaling, according to (6) which –
if correct – would correspond to a linear behavior of data in Fig. S6. Looking
at the data, the trend holds well at large $L$, trending slightly down (i.e.
to a smaller power of $1/L$) at small system sizes on the right of the plot.
As expected from critical coarsening, this trend is most pronounced for
$\Omega$ values closest to the transition point. There a given $L$ is
effectively further from the thermodynamic limit than other $\Omega$ values.
Fig. S7, shows the results of fits to the ansatz (6).

Figure S6: Finite size scaling towards $L\to\infty$ for the $S_{\pi,\pi}$
observable at several values of $\Omega$ near the Dicke SS/SF phase
transition. Circles are simulation data with connecting colored lines. The
black dashed lines are linear fits to the three highest $L$ values. Parameters
and data are the same as shown in Fig. 4.

The growth of the density modulation order parameter (static structure factor
$S_{\pi,\pi}$) is very well described as a linear growth at rate $G$ above the
transition point:

$S_{\pi,\pi}=\left\\{\begin{array}[]{c@{\qquad\text{if}\
}l}G(\Omega-\Omega_{s})&\Omega>\Omega_{s}\\\
0&\Omega\leq\Omega_{s}\end{array}\right.$ (S15)

This is shown in Fig. S7a, along with a fit. The fit gives $G=51.3$,
$\Omega_{s}=2.265$. An error estimate can be obtained from the difference when
compared to a fit over only the four rightmost data points $\Omega\geq 2.285$
instead of five. Taking twice this difference as the error estimate, one
arrives at

$\Omega_{s}=2.265(5)\qquad;\qquad G=51.3(9).$ (S16)

The scale coefficient $C\Omega)$ is peaked around the transition point as
expected.

Figure S7: Fitted parameters to Eq. (6) around the SF/SS phase transition.
Top: $A(\Omega)$, the thermodynamic limit of $S_{\pi,\pi}$, as grey squares
with interpolation between as a blue line. The dashed line is a fit to (S15)
using the rightmost five data points $\Omega\geq 2.275$. The bottom panel
plots the finite scale susceptibility $C(\Omega)$.

Details of the behavior of the magnetization $\Delta$ across the transition
point are shown in Fig. S8. The distribution in panel b at $\Omega=2.26$, just
below the transition corresponds to negative Binder cumulants. The cause of
$U_{L}<0$ here are the fatter tails compared to the Gaussian distribution seen
at lower $\Omega=2.25$. This means that a few realizations are already
achieving non-negligible ordering $|\Delta|\gtrsim 0.2N$. The panels c-d which
proceed across the transition show a broadening of the distribution to include
successively more ordered realizations with $|\Delta|>0.2N$ but still without
the exclusion of intermediate unordered cases that is only seen once the
transition has been unambiguously crossed at $\Omega=2.325$ in panel e.

Figure S8: Behavior of the $\Delta$ effective magnetization near the Dicke
transition. (a-e) show the $\Delta$ distributions for $L=32$ around the
transition, while (f) shows detail of the Binder cumulant $U_{L}$ in the
vicinity.

##  S4 Condensate fraction near the BKT phase transition

The condensate fraction is defined as the ratio $N_{0}/N$, where
$N_{k=0}=\sum_{r}\langle b^{\dagger}_{0}b_{r}\rangle$ is the number of bosons
in the zero momentum state, forming the condensate. It is the analogue of the
susceptibility in spin models and has been measured in current experiments
with cold atoms in optical cavities [4]. In Fig. S9 we display the behavior of
the condensate fraction at long times in the vicinity of the BKT phase
transition as a function of the Rabi frequency for different system sizes.
Differently from the correlation ratio, the condensate fraction scales to zero
with system size _even_ in the QLRO phase. This is consistent with the
Hohemberg-Mermin-Wagner theorem, stating that no true condensation can exist
in the thermodynamic limit for 2D systems with continuous U(1) symmetry. This
result is well known for models at thermal equilibrium, like the XY model or
the Bose-Hubbard model, but it has also been verified for non-equilibrium 2D
systems, like driven-dissipative photon-polaritons [16].

Figure S9: Main panel: condensate fraction $N_{0}/N$ as a function of the Rabi
frequency for increasing system sizes. The solid lines are Padé approximant
fits to the numerical data that we use to extract the derivative of the data
curve, from which we extract inflection points $\Omega^{*}(L)$. Inset: scaling
behavior of the inflection points based on Eq. (8), with $\Omega_{KT}(L)$
replaced by $\Omega^{*}(L)$, yielding $\Omega_{c}/J=8.17(14)$.

Importantly, we see from Fig. S9 that the data curves become steeper and
steeper in the vicinity of $\Omega\approx 8J-9J$ as the system size increases.
In the thermodynamic limit the condensate fraction would vanish, with an
infinite derivative at the BKT critical point. As previously shown for the 2D
Bose-Hubbard model at finite temperature [43], one can approximately estimate
the position of the critical point by studying the scaling behavior of the
inflection points $\Omega^{*}(L)$ in the data curves of Fig. S9. The
correlation length $\xi$ at the inflection point must be of the order of the
system size, that is $\xi(\Omega=\Omega^{*}(L))=bL$, where $b$ is a constant
of order unity. By replacing $\xi$ with the BKT form from Eq. (7), we recover
Eq. (8) with $\Omega_{KT}(L)$ replaced by $\Omega^{*}(L)$. We can therefore
extrapolate the critical Rabi frequency by fitting the inflection points data
based on Eq. (8). The results obtained are shown in the inset of Fig. S9. In
particular we find $\Omega_{c}/J=8.17(14)$, which is consistent with, but less
precise, than the estimate based on $r_{L}$ reported in the main text.