# Emergent finite frequency criticality of driven-dissipative correlated lattice bosons

Orazio Scarlatella, Rosario Fazio, Marco Schiró

► To cite this version:

Orazio Scarlatella, Rosario Fazio, Marco Schiró. Emergent finite frequency criticality of driven-dissipative correlated lattice bosons. Physical Review B, 2019, 99 (6), pp.064511. 10.1103/Phys-RevB.99.064511. cea-04542791

HAL Id: cea-04542791

https://cea.hal.science/cea-04542791

Submitted on 11 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Emergent finite frequency criticality of driven-dissipative correlated lattice bosons

Orazio Scarlatella, Rosario Fazio, Marco Schiró

► To cite this version:

Orazio Scarlatella, Rosario Fazio, Marco Schiró. Emergent finite frequency criticality of driven-dissipative correlated lattice bosons. Physical Review B, 2019, 99 (6), pp.064511. 10.1103/Phys-RevB.99.064511. cea-04542791

HAL Id: cea-04542791

https://cea.hal.science/cea-04542791v1

Submitted on 11 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Emergent Finite Frequency Criticality of Driven-Dissipative Correlated Lattice Bosons

Orazio Scarlatella,¹ Rosario Fazio,²,³ and Marco Schiró¹

¹Institut de Physique Théorique, Université Paris Saclay, CNRS, CEA, F-91191 Gif-sur-Yvette, France

²ICTP, Strada Costiera 11, 34151 Trieste, Italy

³NEST, Scuola Normale Superiore & Istituto Nanoscienze-CNR, I-56126, Pisa, Italy

Critical points and phase transitions are characterized by diverging susceptibilities, reflecting the tendency of the system toward spontaneous symmetry breaking. Equilibrium statistical mechanics bounds these instabilities to occur at zero frequency, giving rise to static order parameters. In this work we argue that a prototype model of correlated driven-dissipative lattice bosons, of direct relevance for upcoming generation of circuit QED arrays experiments, exhibits a susceptibility sharply diverging at a finite non-zero frequency, which is an emerging scale set by interactions and non-equilibrium effects. In the broken-symmetry phase the corresponding macroscopic order parameter becomes non-stationary and oscillates in time without damping, thus breaking continuous time-translational symmetry. Our work, connecting breaking of time translational invariance to divergent finite frequency susceptibilities, which are of direct physical relevance, could potentially be extended to study other time-domain instabilities in non-equilibrium quantum systems, including Floquet time crystals and quantum synchronization.

# I. INTRODUCTION

Classical and quantum phase transitions in systems at thermal equilibrium are characterized, according to the Landau paradigm, by the emergence of a static order parameter which spontaneously breaks a symmetry of the system, such as spin rotational invariance for magnetism or spatial translational invariance for crystals¹,². The resulting criticality is described in terms of an instability of the normal symmetric phase, characterized by a singularity of a static susceptibility. For classical systems far away from thermal equilibrium, such as in presence of external forcing and dissipation, the variety of instabilities can be far richer, with both finite momentum and finite frequency modes going unstable and resulting in the formation of patterns, propagating fronts, spatio-temporal chaos, synchronization or other oscillatory behaviors³⁻⁵. In the quantum world, the question of whether finite frequency modes can become critical, giving rise to time-domain instabilities of the quantum dynamics and to an associated breaking of time-translational invariance, is much less explored. Experimental breakthroughs have brought forth a number of platforms which naturally probe non-equilibrium quantum many body dynamics, ranging from ultra cold atoms⁶, trapped ions⁷ and arrays of non-linear circuit QED cavities⁸⁻¹¹, thus making the question of experimental relevance. In this respect, quantum many body systems in presence of both driving and dissipation mechanisms¹²⁻¹⁴ represent natural platforms to understand and explore such time-domain dynamical instabilities. A well known example is provided by exciton-polariton condensates where superfluidity emerges with an order parameter oscillating in time¹⁵⁻¹⁸. Yet those systems are successfully described by semiclassical theories such as driven-dissipative Gross-Pitaevski equations, leaving open the question about the quantum nature of this phenomenon. More recently the attention has shifted toward strongly correlated quantum lattice models with drive and dissipation, where several works have revealed

the existence of limit cycles, i.e. non-stationary solutions of the quantum dynamics for a macroscopic order parameter, at least at the mean field level¹⁹⁻²⁴

In this work we focus on a paradigmatic model of driven-dissipative interacting bosons on a lattice, which is directly relevant for the upcoming generation of circuit QED arrays experiments targeting Mott insulators of polaritons²⁵⁻²⁷. We argue that a dynamical susceptibility of such an open quantum many body system, which in thermal equilibrium is finite and small since non-zero frequency modes are typically damped by interactions, can display a genuine singularity at finite frequency, as a result of strong interactions and non-equilibrium effects. The critical frequency is not fixed a priori, but rather an emerging scale set by the microscopic parameters. Eventually, the system undergoes a dynamical phase transition where the order parameter in the broken symmetry phase becomes non-stationary and oscillates in time without damping, thus breaking the continuous time-translational symmetry. This stationary-state instability is controlled by both dissipative and coherent couplings, in particular by the ratio between hopping and local interaction, thus providing the strongly correlated analogue of weak coupling non-equilibrium bosons condensation. We organize the manuscript in the following way: in section II we introduce the many-body model and we discuss the physics of its single-site. In section III we discuss the stationary-state instability of the normal phase characterized by a diverging finite frequency susceptibility. We present the consequences of this instability on the dissipative evolution in section IV and finally, in section V, we recover this dynamics from the saddle point solution of the effective non-equilibrium field theory which describes the transition from the strongly correlated regime. In appendix A, we discuss results which are related to the specific pump/loss mechanisms and we consider an alternative driving scheme.
Emergent Finite Frequency Criticality of Driven-Dissipative Correlated Lattice Bosons

Orazio Scarlatella,¹ Rosario Fazio,²,³ and Marco Schiró¹

¹Institut de Physique Théorique, Université Paris Saclay, CNRS, CEA, F-91191 Gif-sur-Yvette, France

²ICTP, Strada Costiera 11, 34151 Trieste, Italy

³NEST, Scuola Normale Superiore & Istituto Nanoscienze-CNR, I-56126, Pisa, Italy

Critical points and phase transitions are characterized by diverging susceptibilities, reflecting the tendency of the system toward spontaneous symmetry breaking. Equilibrium statistical mechanics bounds these instabilities to occur at zero frequency, giving rise to static order parameters. In this work we argue that a prototype model of correlated driven-dissipative lattice bosons, of direct relevance for upcoming generation of circuit QED arrays experiments, exhibits a susceptibility sharply diverging at a finite non-zero frequency, which is an emerging scale set by interactions and non-equilibrium effects. In the broken-symmetry phase the corresponding macroscopic order parameter becomes non-stationary and oscillates in time without damping, thus breaking continuous time-translational symmetry. Our work, connecting breaking of time translational invariance to divergent finite frequency susceptibilities, which are of direct physical relevance, could potentially be extended to study other time-domain instabilities in non-equilibrium quantum systems, including Floquet time crystals and quantum synchronization.

# I. INTRODUCTION

Classical and quantum phase transitions in systems at thermal equilibrium are characterized, according to the Landau paradigm, by the emergence of a static order parameter which spontaneously breaks a symmetry of the system, such as spin rotational invariance for magnetism or spatial translational invariance for crystals¹,². The resulting criticality is described in terms of an instability of the normal symmetric phase, characterized by a singularity of a static susceptibility. For classical systems far away from thermal equilibrium, such as in presence of external forcing and dissipation, the variety of instabilities can be far richer, with both finite momentum and finite frequency modes going unstable and resulting in the formation of patterns, propagating fronts, spatio-temporal chaos, synchronization or other oscillatory behaviors³⁻⁵. In the quantum world, the question of whether finite frequency modes can become critical, giving rise to time-domain instabilities of the quantum dynamics and to an associated breaking of time-translational invariance, is much less explored. Experimental breakthroughs have brought forth a number of platforms which naturally probe non-equilibrium quantum many body dynamics, ranging from ultra cold atoms⁶, trapped ions⁷ and arrays of non-linear circuit QED cavities⁸⁻¹¹, thus making the question of experimental relevance. In this respect, quantum many body systems in presence of both driving and dissipation mechanisms¹²⁻¹⁴ represent natural platforms to understand and explore such time-domain dynamical instabilities. A well known example is provided by exciton-polariton condensates where superfluidity emerges with an order parameter oscillating in time¹⁵⁻¹⁸. Yet those systems are successfully described by semiclassical theories such as driven-dissipative Gross-Pitaevski equations, leaving open the question about the quantum nature of this phenomenon. More recently the attention has shifted toward strongly correlated quantum lattice models with drive and dissipation, where several works have revealed

the existence of limit cycles, i.e. non-stationary solutions of the quantum dynamics for a macroscopic order parameter, at least at the mean field level¹⁹⁻²⁴

In this work we focus on a paradigmatic model of driven-dissipative interacting bosons on a lattice, which is directly relevant for the upcoming generation of circuit QED arrays experiments targeting Mott insulators of polaritons²⁵⁻²⁷. We argue that a dynamical susceptibility of such an open quantum many body system, which in thermal equilibrium is finite and small since non-zero frequency modes are typically damped by interactions, can display a genuine singularity at finite frequency, as a result of strong interactions and non-equilibrium effects. The critical frequency is not fixed a priori, but rather an emerging scale set by the microscopic parameters. Eventually, the system undergoes a dynamical phase transition where the order parameter in the broken symmetry phase becomes non-stationary and oscillates in time without damping, thus breaking the continuous time-translational symmetry. This stationary-state instability is controlled by both dissipative and coherent couplings, in particular by the ratio between hopping and local interaction, thus providing the strongly correlated analogue of weak coupling non-equilibrium bosons condensation. We organize the manuscript in the following way: in section II we introduce the many-body model and we discuss the physics of its single-site. In section III we discuss the stationary-state instability of the normal phase characterized by a diverging finite frequency susceptibility. We present the consequences of this instability on the dissipative evolution in section IV and finally, in section V, we recover this dynamics from the saddle point solution of the effective non-equilibrium field theory which describes the transition from the strongly correlated regime. In appendix A, we discuss results which are related to the specific pump/loss mechanisms and we consider an alternative driving scheme.
## II. MODEL AND DRIVING

We consider the Bose-Hubbard (BH) Hamiltonian

$$H = -J \sum_{\langle ij \rangle} (a_i^\dagger a_j + hc) + \sum_i \left( \delta\omega_0 n_i + \frac{U}{2} n_i^2 \right) \quad (1)$$

modelling a lattice of circuit QED resonators$^{25,28,29}$ hosting a single bosonic cavity mode on each site, $a_i, a_i^\dagger$, with frequency $\delta\omega_0$, repulsive local interaction $U$ and hopping rate $J$. We supplement the model with pump and losses, necessary to drive the system into a non-trivial steady state, and we remark that the qualitative features we are going to discuss do not depend on the specific implementation of drive and dissipation. Here we describe the driven-dissipative dynamics in terms of a master equation for the system density matrix

$$i\partial_t \rho = -i[H, \rho] + \mathcal{D}[\rho] \quad (2)$$

where the dissipator takes the form $\mathcal{D}[\rho] = \mathcal{D}_{loss}[\rho] + \mathcal{D}_{pump}[\rho]$. The first term describes single particle losses with rate $\kappa$

$$\mathcal{D}_{loss}[\rho] = \kappa \sum_i \left( a_i \rho a_i^\dagger - \frac{1}{2} \{ a_i^\dagger a_i, \rho \} \right) \quad (3)$$

while the second term

$$\mathcal{D}_{pump}[\rho] = \sum_i f_{in} \tilde{\mathcal{D}}[a_i^\dagger, \tilde{a}_{i\sigma}; \rho] + f_{out} \tilde{\mathcal{D}}[a_i, \tilde{a}_{i\sigma}^\dagger; \rho] \quad (4)$$

accounts for the coupling to an incoherent bath with a finite bandwidth $\sigma$, which injects and removes particles without a well defined phase in each cavity, at rates $f_{in}, f_{out}$. Here we have introduced the modified dissipator$^{30,31}$

$$\tilde{\mathcal{D}}[X, Y] = X\rho Y + Y^\dagger\rho X^\dagger - X^\dagger Y^\dagger\rho - \rho YX \quad (5)$$

and

$$\tilde{a}_{\sigma}^{\dagger} = \sum_{n} \sqrt{n+1} C_{\sigma}^{*} (\Delta_{n+1}) |n+1\rangle\langle n| \quad (6)$$

which is the photon operator dressed by the finite-bandwidth drive, with $\Delta_{n+1} = Un + (\delta\omega_0 + U/2)$ the level spacing of an isolated single site in Eq. (1). All the results in the main text are obtained with a microscopic model of pumping describing an ensemble of driven two-level emitters embedded in each cavity and undergoing a population inversion, as recently proposed$^{24,31,32}$, which results in $f_{in} = f, f_{out} = 0$. We discuss in the appendix A other microscopic models of incoherent pumping$^{30}$ and the effects of changing the driving protocol. We notice that driving the cavities incoherently preserves the $U(1)$ symmetry, the invariance under $a_i \to e^{i\phi}a_i$, of the Hamiltonian in Eq. (1), thus leaving open the possibility for a non-equilibrium phase transition between a Mott-like incoherent phase and a superfluid, which is well understood in equilibrium$^{2,33}$.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Local bosonic occupation in the interacting single-site driven dissipative problem, respectively as a function of the drive bandwidth $\sigma$ (top panel) showing a characteristic staircase structure and drive amplitude $f$ (bottom panel). The value of the local occupation in the plateaux depends strongly on the nature of the driving protocol and can in general be tuned continuously with the drive amplitude $f$. Parameters: loss rate $\kappa = 0.0128U$, resonator frequency $\delta\omega_0 = 0.0$, drive amplitude (top panel) $f = 0.0625U$.</figcaption></figure>

### A. The single-site problem

To gain some insight on the problem, we start considering the limit $J = 0$, where the lattice problem reduces to a collection of interacting, driven-dissipative decoupled sites. While consisting only of a single site, this quantum problem remains quite non-trivial due to non-linearity and dissipation effects and cannot be in general solved analytically, as in the equilibrium case$^2$. The balance between drive, dissipation and interaction results in a finite number of bosons per site, which however remain incoherent, $\langle a_i \rangle = 0$. The boson number in the single-site problem as a function of pump bandwidth $\sigma$, plotted in figure 1 (top panel), shows a starcaise structure characteristic of blockade physics$^{30}$, with a value of $\Delta\sigma \sim U$ required to add extra bosons in the system. This can be understood naturally: to add a boson in an interacting site the drive has to provide extra energy, however since the system is ultimately open and boson number is not exactly conserved, the exact occupancy will be fixed by the ratio of pump and losses. The drive is able to fix the occupancy to almost integer filling, reflecting the fact that the stationary density matrix is almost pure and resembling the ground state physics of a Bose-Hubbard interacting site, but we stress that the average density in this open and dissipative implementation is never exactly integer, due to losses$^{31}$. The width of the steps is set by the interaction while the height can be tuned continuously by changing pump amplitude $f$, as we show in the bottom panel of figure 1.

A key role in this work is played by the single site Green's functions (retarded/advanced/Keldysh) evalu-
## II. MODEL AND DRIVING

We consider the Bose-Hubbard (BH) Hamiltonian

$$H = -J \sum_{\langle ij \rangle} (a_i^\dagger a_j + hc) + \sum_i \left( \delta\omega_0 n_i + \frac{U}{2} n_i^2 \right) \quad (1)$$

modelling a lattice of circuit QED resonators$^{25,28,29}$ hosting a single bosonic cavity mode on each site, $a_i, a_i^\dagger$, with frequency $\delta\omega_0$, repulsive local interaction $U$ and hopping rate $J$. We supplement the model with pump and losses, necessary to drive the system into a non-trivial steady state, and we remark that the qualitative features we are going to discuss do not depend on the specific implementation of drive and dissipation. Here we describe the driven-dissipative dynamics in terms of a master equation for the system density matrix

$$i\partial_t \rho = -i[H, \rho] + \mathcal{D}[\rho] \quad (2)$$

where the dissipator takes the form $\mathcal{D}[\rho] = \mathcal{D}_{loss}[\rho] + \mathcal{D}_{pump}[\rho]$. The first term describes single particle losses with rate $\kappa$

$$\mathcal{D}_{loss}[\rho] = \kappa \sum_i \left( a_i \rho a_i^\dagger - \frac{1}{2} \{ a_i^\dagger a_i, \rho \} \right) \quad (3)$$

while the second term

$$\mathcal{D}_{pump}[\rho] = \sum_i f_{in} \tilde{\mathcal{D}}[a_i^\dagger, \tilde{a}_{i\sigma}; \rho] + f_{out} \tilde{\mathcal{D}}[a_i, \tilde{a}_{i\sigma}^\dagger; \rho] \quad (4)$$

accounts for the coupling to an incoherent bath with a finite bandwidth $\sigma$, which injects and removes particles without a well defined phase in each cavity, at rates $f_{in}, f_{out}$. Here we have introduced the modified dissipator$^{30,31}$

$$\tilde{\mathcal{D}}[X, Y] = X\rho Y + Y^\dagger\rho X^\dagger - X^\dagger Y^\dagger\rho - \rho YX \quad (5)$$

and

$$\tilde{a}_{\sigma}^{\dagger} = \sum_{n} \sqrt{n+1} C_{\sigma}^{*} (\Delta_{n+1}) |n+1\rangle\langle n| \quad (6)$$

which is the photon operator dressed by the finite-bandwidth drive, with $\Delta_{n+1} = Un + (\delta\omega_0 + U/2)$ the level spacing of an isolated single site in Eq. (1). All the results in the main text are obtained with a microscopic model of pumping describing an ensemble of driven two-level emitters embedded in each cavity and undergoing a population inversion, as recently proposed$^{24,31,32}$, which results in $f_{in} = f, f_{out} = 0$. We discuss in the appendix A other microscopic models of incoherent pumping$^{30}$ and the effects of changing the driving protocol. We notice that driving the cavities incoherently preserves the $U(1)$ symmetry, the invariance under $a_i \to e^{i\phi}a_i$, of the Hamiltonian in Eq. (1), thus leaving open the possibility for a non-equilibrium phase transition between a Mott-like incoherent phase and a superfluid, which is well understood in equilibrium$^{2,33}$.

<figure><img src="image_3.png" /><figcaption>FIG. 1. Local bosonic occupation in the interacting single-site driven dissipative problem, respectively as a function of the drive bandwidth $\sigma$ (top panel) showing a characteristic staircase structure and drive amplitude $f$ (bottom panel). The value of the local occupation in the plateaux depends strongly on the nature of the driving protocol and can in general be tuned continuously with the drive amplitude $f$. Parameters: loss rate $\kappa = 0.0128U$, resonator frequency $\delta\omega_0 = 0.0$, drive amplitude (top panel) $f = 0.0625U$.</figcaption></figure>

### A. The single-site problem

To gain some insight on the problem, we start considering the limit $J = 0$, where the lattice problem reduces to a collection of interacting, driven-dissipative decoupled sites. While consisting only of a single site, this quantum problem remains quite non-trivial due to non-linearity and dissipation effects and cannot be in general solved analytically, as in the equilibrium case$^2$. The balance between drive, dissipation and interaction results in a finite number of bosons per site, which however remain incoherent, $\langle a_i \rangle = 0$. The boson number in the single-site problem as a function of pump bandwidth $\sigma$, plotted in figure 1 (top panel), shows a starcaise structure characteristic of blockade physics$^{30}$, with a value of $\Delta\sigma \sim U$ required to add extra bosons in the system. This can be understood naturally: to add a boson in an interacting site the drive has to provide extra energy, however since the system is ultimately open and boson number is not exactly conserved, the exact occupancy will be fixed by the ratio of pump and losses. The drive is able to fix the occupancy to almost integer filling, reflecting the fact that the stationary density matrix is almost pure and resembling the ground state physics of a Bose-Hubbard interacting site, but we stress that the average density in this open and dissipative implementation is never exactly integer, due to losses$^{31}$. The width of the steps is set by the interaction while the height can be tuned continuously by changing pump amplitude $f$, as we show in the bottom panel of figure 1.

A key role in this work is played by the single site Green's functions (retarded/advanced/Keldysh) evalu-
<figure><img src="image_3.png" /><figcaption>FIG. 2. Local Spectral Function and the Emergent Frequency Ω*. Retarded Green's function, see equation (7) of the interacting single-site problem in Eq. (1) with drive and dissipation. The imaginary part, describing the spectral function, shows two peaks separated by a gap of order U. Crucially, the two peaks are not symmetrically placed around zero frequency, i.e. the imaginary part changes its sign at a finite frequency Ω*. This emergent frequency is not fixed a priori but rather fully tunable and depending from the amplitude and bandwidth of the drive. Parameters: Drive amplitude f = 0.0625U, bandwidth σ = 1.5U, loss rate κ = 0.0128U, resonator frequency δω₀ = 0.0.</figcaption></figure>

ated in the stationary state which are defined as

$$G_{\text{loc}}^{R/K}(\omega) = -i \int_{0}^{\infty} dt e^{i\omega t} \langle [a(t), a^{\dagger}(0)]_{\mp} \rangle_{\text{loc}}. \quad (7)$$

and $G_{\text{loc}}^A(\omega) = G_{\text{loc}}^R(\omega)^*$, which we evaluate using their Lehmann representation³⁴ and a numerical diagonalization of the Liouvillian. As we see in the bottom panel of figure 2, when the interaction U is sufficiently strong with respect to dissipation, the lorentzian-shaped spectrum of a driven-dissipative oscillator splits into two atomic-like excitations separated by a large gap U. Remarkably, this spectral function is not constrained to change sign at zero frequency as in thermal equilibrium, but it does it at a non-zero energy scale $\Omega_* \neq 0$ which depends on interactions, drive and dissipation. In figure 2 we plot this quantity as a function of the drive bandwidth σ and drive amplitude f, showing that $\Omega_*$ increases with σ and f, together with the bosonic occupation. This frequency can be interpreted as an emergent chemical potential for the bosons, as we can deduce by considering the bosonic distribution function, which contains information on the occupation of bosonic modes. We can define a bosonic distribution function in analogy with the thermal equilibrium case, i.e.

$$F(\omega) = \frac{G_{\text{loc}}^{K}(\omega)}{G_{\text{loc}}^{R}(\omega) - G_{\text{loc}}^{A}(\omega)} \quad (8)$$

where $G_{\text{loc}}^K(\omega)$ is the Keldysh Green's function. Indeed in thermal equilibrium the fluctuation-dissipation theorem

<figure><img src="image_4.png" /><figcaption>FIG. 3. Effective Distribution Function as a function of frequency around $\Omega_*$ for fixed value of interaction, drive and dissipation. Around the frequency $\Omega_*$ the system develops a singularity which allows to define an effective temperature $T_{\text{eff}}$, whose dependence from the drive bandwidth is plotted in the inset. Parameters: loss rate $\kappa = 0.0128U$, resonator frequency $\delta\omega_0 = 0.0$, drive bandwidth (left panel) $\sigma = 0.625U$, drive amplitude (right panel) $f = 0.0625U$.</figcaption></figure>

constraints the functional form of the distribution function to the canonical bosonic one, $F_{\text{eq}}(\omega) = \coth \beta\omega/2$, which at low frequency (or high temperature) becomes $F_{\text{eq}}(\omega) \approx T/2\omega$. Here, in presence of drive and dissipation, such an identity does not hold and we use Eq.(8) as operational definition of the distribution function. We plot in figure 3 the distribution function for a given value of interaction, drive and dissipation. While its overall shape shows departure from the thermal equilibrium case, we find that around the critical frequency $\Omega_*$ the system develops a singularity of the form $F(\omega) \sim T_{\text{eff}}/(\omega - \Omega_*)$, ultimately arising from the fact that $\text{Im}G_{\text{loc}}^R(\omega)$ has a zero at $\Omega_*$ while the Keldysh component is finite around the same frequency range. This suggests an asymptotic thermalization around the frequency $\Omega_*$, with a small effective temperature $T_{\text{eff}}$ weakly decreasing with the drive bandwidth (see inset of figure 3) and with the frequency $\Omega_*$ playing now the role of an effective chemical potential.

### III. INSTABILITY OF NORMAL PHASE

We now turn to discuss the properties of the lattice model in presence of a finite hopping J. To this extent we introduce a general approach to study the instabilities of the normal incoherent phase of driven-dissipative correlated lattice models, which generalizes the equilibrium strong-coupling approach of Ref³³. We start writing the Keldysh action associated to the many-body quantum master equation (2), which allows to describe the non-equilibrium stationary state and the excitations on top of it, and decouple the hopping term by means of an auxiliary bosonic field playing the role of local order
<figure><img src="image_4.png" /><figcaption>FIG. 2. Local Spectral Function and the Emergent Frequency Ω*. Retarded Green's function, see equation (7) of the interacting single-site problem in Eq. (1) with drive and dissipation. The imaginary part, describing the spectral function, shows two peaks separated by a gap of order U. Crucially, the two peaks are not symmetrically placed around zero frequency, i.e. the imaginary part changes its sign at a finite frequency Ω*. This emergent frequency is not fixed a priori but rather fully tunable and depending from the amplitude and bandwidth of the drive. Parameters: Drive amplitude f = 0.0625U, bandwidth σ = 1.5U, loss rate κ = 0.0128U, resonator frequency δω₀ = 0.0.</figcaption></figure>

ated in the stationary state which are defined as

$$G_{\text{loc}}^{R/K}(\omega) = -i \int_{0}^{\infty} dt e^{i\omega t} \langle [a(t), a^{\dagger}(0)]_{\mp} \rangle_{\text{loc}}. \quad (7)$$

and $G_{\text{loc}}^{A}(\omega) = G_{\text{loc}}^{R}(\omega)^*$, which we evaluate using their Lehmann representation³⁴ and a numerical diagonalization of the Liouvillian. As we see in the bottom panel of figure 2, when the interaction U is sufficiently strong with respect to dissipation, the lorentzian-shaped spectrum of a driven-dissipative oscillator splits into two atomic-like excitations separated by a large gap U. Remarkably, this spectral function is not constrained to change sign at zero frequency as in thermal equilibrium, but it does it at a non-zero energy scale $\Omega_* \neq 0$ which depends on interactions, drive and dissipation. In figure 2 we plot this quantity as a function of the drive bandwidth σ and drive amplitude f, showing that $\Omega_*$ increases with σ and f, together with the bosonic occupation. This frequency can be interpreted as an emergent chemical potential for the bosons, as we can deduce by considering the bosonic distribution function, which contains information on the occupation of bosonic modes. We can define a bosonic distribution function in analogy with the thermal equilibrium case, i.e.

$$F(\omega) = \frac{G_{\text{loc}}^{K}(\omega)}{G_{\text{loc}}^{R}(\omega) - G_{\text{loc}}^{A}(\omega)} \quad (8)$$

where $G_{\text{loc}}^{K}(\omega)$ is the Keldysh Green's function. Indeed in thermal equilibrium the fluctuation-dissipation theorem

<figure><img src="image_5.png" /><figcaption>FIG. 3. Effective Distribution Function as a function of frequency around $\Omega_*$ for fixed value of interaction, drive and dissipation. Around the frequency $\Omega_*$ the system develops a singularity which allows to define an effective temperature $T_{\text{eff}}$, whose dependence from the drive bandwidth is plotted in the inset. Parameters: loss rate $\kappa = 0.0128U$, resonator frequency $\delta\omega_0 = 0.0$, drive bandwidth (left panel) $\sigma = 0.625U$, drive amplitude (right panel) $f = 0.0625U$.</figcaption></figure>

constraints the functional form of the distribution function to the canonical bosonic one, $F_{\text{eq}}(\omega) = \coth \beta\omega/2$, which at low frequency (or high temperature) becomes $F_{\text{eq}}(\omega) \approx T/2\omega$. Here, in presence of drive and dissipation, such an identity does not hold and we use Eq.(8) as operational definition of the distribution function. We plot in figure 3 the distribution function for a given value of interaction, drive and dissipation. While its overall shape shows departure from the thermal equilibrium case, we find that around the critical frequency $\Omega_*$ the system develops a singularity of the form $F(\omega) \sim T_{\text{eff}}/(\omega - \Omega_*)$, ultimately arising from the fact that $\text{Im}G_{\text{loc}}^{R}(\omega)$ has a zero at $\Omega_*$ while the Keldysh component is finite around the same frequency range. This suggests an asymptotic thermalization around the frequency $\Omega_*$, with a small effective temperature $T_{\text{eff}}$ weakly decreasing with the drive bandwidth (see inset of figure 3) and with the frequency $\Omega_*$ playing now the role of an effective chemical potential.

### III. INSTABILITY OF NORMAL PHASE

We now turn to discuss the properties of the lattice model in presence of a finite hopping J. To this extent we introduce a general approach to study the instabilities of the normal incoherent phase of driven-dissipative correlated lattice models, which generalizes the equilibrium strong-coupling approach of Ref³³. We start writing the Keldysh action associated to the many-body quantum master equation (2), which allows to describe the non-equilibrium stationary state and the excitations on top of it, and decouple the hopping term by means of an auxiliary bosonic field playing the role of local order
parameter. The resulting effective action takes the form

$$S_{eff} = \int_C dt \sum_{ij} \psi_i^\dagger J_{ij}^{-1} \psi_j + \sum_i \Gamma[\psi_i^\dagger, \psi_i] \quad (9)$$

where the second term represents the generating functional of the local bosonic Green's functions, $\Gamma[\psi_i^\dagger, \psi_i] = \log\langle T_C e^{i \int_C dt (\psi_i^\dagger a_i + a_i^\dagger \psi_i)} \rangle_{loc}$, with the average taken over the interacting driven-dissipative single-site problem. Our approach therefore combines the strong coupling field theory$^{33}$ with the exact numerical solution of the single-site problem based on a diagonalization of its Liouvillian. As such, it could be applied to lattice models with any scheme of incoherent local drive and dissipation by just solving the appropriate local problem.

Expanding Eq. (9) in the fields $\psi, \psi^\dagger$ and within a gaussian approximation, we obtain the effective action $S_{eff} = \int d\omega d^d q \bar{\Psi}(q, \omega) \hat{G}^{-1}(q, \omega) \Psi(q, \omega)$, where we moved to momentum space, we defined $\bar{\Psi}(q, \omega) = (\psi_c^*(q, \omega) \psi_q^*(q, \omega))$ with $\psi_c/q$ the classical/quantum fields, and we introduced

$$\hat{G}^{-1}(q, \omega) = \begin{pmatrix} 0 & J_q^{-1} - G_{loc}^A(\omega) \\ J_q^{-1} - G_{loc}^R(\omega) & -G_{loc}^K(\omega) \end{pmatrix}$$

In the above expression we have $J_q = -2J \sum_{\alpha=1}^d \cos q_\alpha$, while $G_{loc}^{R/A/K}(\omega)$ are the exact single-site regarded/advanced/Keldysh Green's functions evaluated in the stationary state

From the effective action $S_{eff}$, the susceptibility $\chi^R(q, \omega)$ of the order parameter reads

$$\chi^R(q, \omega) = \frac{1}{J_q^{-1} - G_{loc}^R(\omega)} \quad (10)$$

The gaussian approximation of the effective action, leading to Eq. (10), is well suited within the normal phase, where the order parameter fluctuates around zero. It is equivalent to a strong-coupling resummation of the perturbation theory in the hopping $J$ often referred to as random phase approximation (RPA) around the atomic limit$^{35,36}$ and it is well known to capture qualitatively the instability of the normal phase, as for example discussed in Refs.$^{2,33}$.

In thermal equilibrium, the $U(1)$ susceptibility Eq. (10) is well known to show a zero-frequency singularity at a critical value of the hopping, at which the Mott insulating phase becomes unstable towards superfluidity$^{2,33}$. As we are going to show, the behavior of the same quantity in a non-equilibrium state is remarkably different. In figure 4 we plot the $q=0$ susceptibility, probing the instability of the homogeneous normal phase, for different values of the hopping strength $J$. We find a well defined resonance structure which gets sharper and narrower as the hopping is increased, and eventually turns into a genuine finite frequency pole at $\omega = \Omega_*$ when a critical hopping $J_c$ is reached. Right at $J_c$ the susceptibility diverges as a power law around $\Omega_*$, $\chi^R(q=0, \omega) = \chi_0 / (\omega - \Omega_*)^\alpha$, with $\alpha = 1$. The appearance of a singularity at finite

<figure><img src="image_5.png" /><figcaption>FIG. 4. Normal phase susceptibility in the gaussian approximation, Eq. (10), as a function of frequency and for different values of the hopping $J$. The real part of the susceptibility goes through a zero with a characteristic two-peaks structure which gets sharper as the critical value $J_c$ is approached, turning into a genuine singularity at $\Omega_*$. This emergent scale is controlled by the local occupation of the single site problem (top panel), whose spectral function is plotted in the bottom inset. Parameters: Drive amplitude $f = 0.0625U$, bandwidth $\sigma = 1.5U$, loss rate $\kappa = 0.0064U$, resonator frequency $\delta\omega_0 = 0.0$.</figcaption></figure>

frequency is a remarkable result with no counterpart in systems in thermal equilibrium, where one expects finite frequency modes to be damped by interactions thus acquiring a finite lifetime, ultimately cutting off the singularity of any dynamical susceptibility. Its origin is rooted in the physics of the single-site quantum problem whose spectral function, $\text{Im}G_{loc}^R(\omega)$, enters the susceptibility through Eq. (10). In particular, as we show in the top panel of figure 4, the critical frequency is set by the local occupancy, rather than by the strength of the order parameter as for the oscillations of weakly interacting non-equilibrium superfluids$^{37}$, a fact which highlights the quantum nature of the incoherent phase becoming unstable at $\Omega_*$.

## IV. DISSIPATIVE DYNAMICS

We now discuss the consequences of the finite-frequency, stationary-state instability we have presented so far and investigate the dissipative dynamics (2) of the lattice problem for different values of the hopping $J$. To this extent, we use a time-dependent Gutzwiller (tdG) decoupling of the density matrix, i.e. $\rho(t) = \prod_i \rho_i(t)$ that we further assume homogeneous in space, $\rho_i(t) \equiv \rho_{loc}(t)$. This approximation results in an effective single-site problem $\partial_t \rho_{loc}(t) = -i[H_{eff}(t), \rho_{loc}(t)] + D_{loc}[\rho_{loc}]$ where $H_{eff}(t) = \delta\omega_0 n + Un^2/2 + zJ(a^\dagger\psi(t) + hc)$ with $z = d/2$ the coordination number of the lattice, $D_{loc}$ is
parameter. The resulting effective action takes the form

$$S_{eff} = \int_C dt \sum_{ij} \psi_i^\dagger J_{ij}^{-1} \psi_j + \sum_i \Gamma[\psi_i^\dagger, \psi_i] \quad (9)$$

where the second term represents the generating functional of the local bosonic Green's functions, $\Gamma[\psi_i^\dagger, \psi_i] = \log\langle T_C e^{i \int_C dt (\psi_i^\dagger a_i + a_i^\dagger \psi_i)} \rangle_{loc}$, with the average taken over the interacting driven-dissipative single-site problem. Our approach therefore combines the strong coupling field theory$^{33}$ with the exact numerical solution of the single-site problem based on a diagonalization of its Liouvillian. As such, it could be applied to lattice models with any scheme of incoherent local drive and dissipation by just solving the appropriate local problem.

Expanding Eq. (9) in the fields $\psi, \psi^\dagger$ and within a gaussian approximation, we obtain the effective action $S_{eff} = \int d\omega d^d q \bar{\Psi}(q, \omega) \hat{G}^{-1}(q, \omega) \Psi(q, \omega)$, where we moved to momentum space, we defined $\bar{\Psi}(q, \omega) = (\psi_c^*(q, \omega) \psi_q^*(q, \omega))$ with $\psi_c/q$ the classical/quantum fields, and we introduced

$$\hat{G}^{-1}(q, \omega) = \begin{pmatrix} 0 & J_q^{-1} - G_{loc}^A(\omega) \\ J_q^{-1} - G_{loc}^R(\omega) & -G_{loc}^K(\omega) \end{pmatrix}$$

In the above expression we have $J_q = -2J \sum_{\alpha=1}^d \cos q_\alpha$, while $G_{loc}^{R/A/K}(\omega)$ are the exact single-site regarded/advanced/Keldysh Green's functions evaluated in the stationary state

From the effective action $S_{eff}$, the susceptibility $\chi^R(q, \omega)$ of the order parameter reads

$$\chi^R(q, \omega) = \frac{1}{J_q^{-1} - G_{loc}^R(\omega)} \quad (10)$$

The gaussian approximation of the effective action, leading to Eq. (10), is well suited within the normal phase, where the order parameter fluctuates around zero. It is equivalent to a strong-coupling resummation of the perturbation theory in the hopping $J$ often referred to as random phase approximation (RPA) around the atomic limit$^{35,36}$ and it is well known to capture qualitatively the instability of the normal phase, as for example discussed in Refs.$^{2,33}$.

In thermal equilibrium, the $U(1)$ susceptibility Eq. (10) is well known to show a zero-frequency singularity at a critical value of the hopping, at which the Mott insulating phase becomes unstable towards superfluidity$^{2,33}$. As we are going to show, the behavior of the same quantity in a non-equilibrium state is remarkably different. In figure 4 we plot the $q=0$ susceptibility, probing the instability of the homogeneous normal phase, for different values of the hopping strength $J$. We find a well defined resonance structure which gets sharper and narrower as the hopping is increased, and eventually turns into a genuine finite frequency pole at $\omega = \Omega_*$ when a critical hopping $J_c$ is reached. Right at $J_c$ the susceptibility diverges as a power law around $\Omega_*$, $\chi^R(q=0, \omega) = \chi_0 / (\omega - \Omega_*)^\alpha$, with $\alpha = 1$. The appearance of a singularity at finite

<figure><img src="image_6.png" /><figcaption>FIG. 4. Normal phase susceptibility in the gaussian approximation, Eq. (10), as a function of frequency and for different values of the hopping $J$. The real part of the susceptibility goes through a zero with a characteristic two-peaks structure which gets sharper as the critical value $J_c$ is approached, turning into a genuine singularity at $\Omega_*$. This emergent scale is controlled by the local occupation of the single site problem (top panel), whose spectral function is plotted in the bottom inset. Parameters: Drive amplitude $f = 0.0625U$, bandwidth $\sigma = 1.5U$, loss rate $\kappa = 0.0064U$, resonator frequency $\delta\omega_0 = 0.0$.</figcaption></figure>

frequency is a remarkable result with no counterpart in systems in thermal equilibrium, where one expects finite frequency modes to be damped by interactions thus acquiring a finite lifetime, ultimately cutting off the singularity of any dynamical susceptibility. Its origin is rooted in the physics of the single-site quantum problem whose spectral function, $\text{Im}G_{loc}^R(\omega)$, enters the susceptibility through Eq. (10). In particular, as we show in the top panel of figure 4, the critical frequency is set by the local occupancy, rather than by the strength of the order parameter as for the oscillations of weakly interacting non-equilibrium superfluids$^{37}$, a fact which highlights the quantum nature of the incoherent phase becoming unstable at $\Omega_*$.

## IV. DISSIPATIVE DYNAMICS

We now discuss the consequences of the finite-frequency, stationary-state instability we have presented so far and investigate the dissipative dynamics (2) of the lattice problem for different values of the hopping $J$. To this extent, we use a time-dependent Gutzwiller (tdG) decoupling of the density matrix, i.e. $\rho(t) = \prod_i \rho_i(t)$ that we further assume homogeneous in space, $\rho_i(t) \equiv \rho_{loc}(t)$. This approximation results in an effective single-site problem $\partial_t \rho_{loc}(t) = -i[H_{eff}(t), \rho_{loc}(t)] + D_{loc}[\rho_{loc}]$ where $H_{eff}(t) = \delta\omega_0 n + Un^2/2 + zJ(a^\dagger\psi(t) + hc)$ with $z = d/2$ the coordination number of the lattice, $D_{loc}$ is
<figure><img src="image_6.png" /><figcaption>FIG. 5. Dynamics of the order parameter $\psi(t) = |\psi(t)|e^{i\theta(t)}$ obtained from the tdG method, for different values of the hopping strength $J$. (Top Panel) Normal phase, $J = 0.8J_c$, exponential decay of the order parameter toward an incoherent stationary state $|\psi(t \to \infty)| \sim e^{-t/\tau(J)}$, with oscillations at frequency $\Omega(J)$ (see inset) due to the phase linearly growing in time. (Bottom Panel, log scale) Approaching the critical point $J_c$, the dynamics slows down suggesting a power-law decay right at the transition, as we show analytically. Parameters: Drive amplitude $f = 0.0625U$, bandwidth $\sigma = 1.5U$, loss rate $\kappa = 0.0064U$, resonator frequency $\delta\omega_0 = 0.0$.</figcaption></figure>

the local dissipator including incoherent drive and losses, while $\psi(t) = \text{Tr}\rho_{\text{loc}}(t)a$ is a self-consistent time dependent field. We expect this approximation to capture some qualitative features of the dynamics across the phase transition, at least in high enough dimensions, where its mean field description is supposed to be accurate.

In figure 5 we plot the dynamics of the bosonic order parameter $\psi(t) = \langle a(t) \rangle$ as a function of time for different values of the hopping $J < J_c$. If we introduce a polar representation, $\psi(t) = |\psi(t)|e^{i\theta(t)}$, we see that the absolute value of the order parameter shows an exponential relaxation toward zero, $|\psi(t)| \sim e^{-t/\tau(J)}$, indicating an incoherent stationary state, while the phase grows linearly in time with finite angular velocity $\Omega$, $\theta(t) = \Omega(J)t + \theta_0$. A closer inspection reveals that the characteristic frequency $\Omega(J)$ differs from the value $\Omega_*$ previously identified by an amount $\delta\Omega(J) = |\Omega(J) - \Omega_*|$ which strongly depends on the hopping rate $J$ and vanishes at the critical point $J = J_c$ with a characteristic power law, $\delta\Omega \sim (J_c - J)$, as shown in the top panel of figure 6. Similarly the relaxation time diverges upon approaching the critical hopping $J_c$, $\tau \sim 1/(J_c - J)$ (see figure 6) and the order parameter shows a characteristic critical slowing down, as shown in the bottom panel of figure 5.

## V. NONEQUILIBRIUM FIELD THEORY OF FINITE FREQUENCY CRITICALITY

We now proceed to set up a Keldysh non-equilibrium field theory for the finite-frequency dissipative transition,

<figure><img src="image_7.png" /><figcaption>FIG. 6. Comparison of the results obtained from the tdG method and from Eq. (12) for the two scales $\delta\Omega$ (slow frequency oscillation mode, top panel) and $\tau$ (relaxation time to reach a steady state, bottom) as a function of the distance from the critical point, $J_c - J$. The dissipative dynamical transition is characterized by both energy scales becoming critical. Parameters: Drive amplitude $f = 0.0625U$, bandwidth $\sigma = 1.5U$, loss rate $\kappa = 0.00648U$, resonator frequency $\delta\omega_0 = 0.0$.</figcaption></figure>

which allows us to obtain a complete analytical picture of the mean field dynamics and sets the stage to discuss the role of quantum fluctuations beyond mean field. The starting point is to expand the effective action (9) for $q \to 0$ and $\omega \to \Omega_*$ and then move to a rotating frame where the field is oscillating at frequency $\Omega_*$. Introducing the fields $\tilde{\psi}_{c,q}(x,t) = e^{-i\Omega_*t}\psi_{c,q}(x,t)$, we obtain

$$S_{\text{eff}} = \int dt dx \tilde{\psi}_c^* \left( -r + K_1 i \partial_t + \frac{K_2}{2} \partial_t^2 - K_3 \nabla^2 \right) \tilde{\psi}_q + hc + S_{\text{noise}} + S_{\text{int}} \quad (11)$$

where $r = 1/zJ + \text{Re}G_{\text{loc}}^R(\Omega_*) = (J_c - J)/J_c^2$ is the distance from the dissipative phase transition while $K_3 = 1/zJ^2$ is the bare mass. Differently from the Gross-Pitaevski weak-coupling regime$^{38}$, the effective action (11) features both first and second time-derivative terms, with complex coefficients $K_{1,2} = \partial_{\omega}^{1,2}G_{\text{loc}}^R(\Omega_*)$, a feature of the strong-coupling limit around which we expand. In equilibrium $K_1, K_2$ play a crucial role for the critical behavior of the transition, which changes universality class at the tip of the Mott lobes, where $K_1 = 0$. In the present case we find $K_1 \neq 0$ along the phase boundary, suggesting a single universality class for our driven-dissipative case. We therefore drop $K_2$ in the following, an assumption which is justified at the mean field level, and leave for future studies the investigation on the role of $K_2$ upon including fluctuations. In Eq. (11), $S_{\text{noise}} = \int dx dtdt' \tilde{\psi}_q^*(x,t) G_{\text{loc}}^K(t-t') \tilde{\psi}_q(x,t')$ represents the noise contribution, which depends on the Keldysh Green's function of the single site. Around $\Omega_*$ this is finite which suggests to disregard retardation and obtain a purely time local quadratic action with noise
<figure><img src="image_7.png" /><figcaption>FIG. 5. Dynamics of the order parameter $\psi(t) = |\psi(t)|e^{i\theta(t)}$ obtained from the tdG method, for different values of the hopping strength $J$. (Top Panel) Normal phase, $J = 0.8J_c$, exponential decay of the order parameter toward an incoherent stationary state $|\psi(t \to \infty)| \sim e^{-t/\tau(J)}$, with oscillations at frequency $\Omega(J)$ (see inset) due to the phase linearly growing in time. (Bottom Panel, log scale) Approaching the critical point $J_c$, the dynamics slows down suggesting a power-law decay right at the transition, as we show analytically. Parameters: Drive amplitude $f = 0.0625U$, bandwidth $\sigma = 1.5U$, loss rate $\kappa = 0.0064U$, resonator frequency $\delta\omega_0 = 0.0$.</figcaption></figure>

the local dissipator including incoherent drive and losses, while $\psi(t) = \text{Tr}\rho_{\text{loc}}(t)a$ is a self-consistent time dependent field. We expect this approximation to capture some qualitative features of the dynamics across the phase transition, at least in high enough dimensions, where its mean field description is supposed to be accurate.

In figure 5 we plot the dynamics of the bosonic order parameter $\psi(t) = \langle a(t) \rangle$ as a function of time for different values of the hopping $J < J_c$. If we introduce a polar representation, $\psi(t) = |\psi(t)|e^{i\theta(t)}$, we see that the absolute value of the order parameter shows an exponential relaxation toward zero, $|\psi(t)| \sim e^{-t/\tau(J)}$, indicating an incoherent stationary state, while the phase grows linearly in time with finite angular velocity $\Omega$, $\theta(t) = \Omega(J)t + \theta_0$. A closer inspection reveals that the characteristic frequency $\Omega(J)$ differs from the value $\Omega_*$ previously identified by an amount $\delta\Omega(J) = |\Omega(J) - \Omega_*|$ which strongly depends on the hopping rate $J$ and vanishes at the critical point $J = J_c$ with a characteristic power law, $\delta\Omega \sim (J_c - J)$, as shown in the top panel of figure 6. Similarly the relaxation time diverges upon approaching the critical hopping $J_c$, $\tau \sim 1/(J_c - J)$ (see figure 6) and the order parameter shows a characteristic critical slowing down, as shown in the bottom panel of figure 5.

## V. NONEQUILIBRIUM FIELD THEORY OF FINITE FREQUENCY CRITICALITY

We now proceed to set up a Keldysh non-equilibrium field theory for the finite-frequency dissipative transition,

<figure><img src="image_8.png" /><figcaption>FIG. 6. Comparison of the results obtained from the tdG method and from Eq. (12) for the two scales $\delta\Omega$ (slow frequency oscillation mode, top panel) and $\tau$ (relaxation time to reach a steady state, bottom) as a function of the distance from the critical point, $J_c - J$. The dissipative dynamical transition is characterized by both energy scales becoming critical. Parameters: Drive amplitude $f = 0.0625U$, bandwidth $\sigma = 1.5U$, loss rate $\kappa = 0.00648U$, resonator frequency $\delta\omega_0 = 0.0$.</figcaption></figure>

which allows us to obtain a complete analytical picture of the mean field dynamics and sets the stage to discuss the role of quantum fluctuations beyond mean field. The starting point is to expand the effective action (9) for $q \to 0$ and $\omega \to \Omega_*$ and then move to a rotating frame where the field is oscillating at frequency $\Omega_*$. Introducing the fields $\tilde{\psi}_{c,q}(x,t) = e^{-i\Omega_*t}\psi_{c,q}(x,t)$, we obtain

$$S_{\text{eff}} = \int dt dx \tilde{\psi}_c^* \left( -r + K_1 i \partial_t + \frac{K_2}{2} \partial_t^2 - K_3 \nabla^2 \right) \tilde{\psi}_q + hc + S_{\text{noise}} + S_{\text{int}} \quad (11)$$

where $r = 1/zJ + \text{Re}G_{\text{loc}}^R(\Omega_*) = (J_c - J)/J_c^2$ is the distance from the dissipative phase transition while $K_3 = 1/zJ^2$ is the bare mass. Differently from the Gross-Pitaevski weak-coupling regime$^{38}$, the effective action (11) features both first and second time-derivative terms, with complex coefficients $K_{1,2} = \partial_{\omega}^{1,2}G_{\text{loc}}^R(\Omega_*)$, a feature of the strong-coupling limit around which we expand. In equilibrium $K_1, K_2$ play a crucial role for the critical behavior of the transition, which changes universality class at the tip of the Mott lobes, where $K_1 = 0$. In the present case we find $K_1 \neq 0$ along the phase boundary, suggesting a single universality class for our driven-dissipative case. We therefore drop $K_2$ in the following, an assumption which is justified at the mean field level, and leave for future studies the investigation on the role of $K_2$ upon including fluctuations. In Eq. (11), $S_{\text{noise}} = \int dx dtdt' \tilde{\psi}_q^*(x,t) G_{\text{loc}}^K(t-t') \tilde{\psi}_q(x,t')$ represents the noise contribution, which depends on the Keldysh Green's function of the single site. Around $\Omega_*$ this is finite which suggests to disregard retardation and obtain a purely time local quadratic action with noise
term $S_{noise} = D \int dxdt \tilde{\psi}_q^*(x, t) \tilde{\psi}_q(x, t)$. Here $D$ plays the role of effective diffusion coefficient in the equivalent stochastic (Langevin) dynamics and it is indeed given by $D \simeq G_{loc}^K(\Omega_*) \sim T_{eff}$. Finally $S_{int}$ accounts for the non-linearities and it is completely determined by the multi-particle Green's functions of the driven-dissipative single site problem. If we restrict ourselves to interaction terms with one quantum and three classical fields, which is valid in high enough dimensions according to canonical power counting$^{38}$, we can write this term as $S_{int} = u \int dxdt \tilde{\psi}_q^*(x, t) \tilde{\psi}_c^*(x, t) \tilde{\psi}_c^2(x, t) + hc$. We can now take the saddle point equation $\delta S/\delta \tilde{\psi}_q^*(x, t) = 0$ and obtain the equation of motion

$$(iK_1\partial_t - K_3\nabla^2 - r) \tilde{\psi}_c + u|\tilde{\psi}_c|^2 \tilde{\psi}_c = 0 \quad (12)$$

which takes the form of a complex Ginzburg-Landau equation, well known as a phenomenological description of pattern formation in classical non-equilibrium systems$^{3-5}$. The spatially homogeneous solution of Equation (12) can be obtained in closed form, as we discuss in the next section, and it describes a transition between a phase where $\tilde{\psi}_c \to 0$ for $t \to \infty$ and a phase where the modulus of the order parameter saturates to a finite value. It further allows us to compute the scales $1/\tau$ and $\delta\Omega$ which perfectly match the numerical results that we have found from the Gutzwiller dynamics in the previous section, as shown in Figure 5. We remark that this classical dynamics only describes the evolution in a frame rotating at frequency $\Omega_*$.

## A. Analytical Solution of Saddle Point Dynamics

We discuss here more in detail the solution of Equation (12). Assuming an homogeneous solution and introducing polar coordinates for absolute value and phase of the order parameters $\tilde{\psi}_c(t) = |\tilde{\psi}_c(t)|e^{i\theta(t)}$ one obtains two separate equations for $|\tilde{\psi}_c(t)|$ and $\tilde{\theta}(t)$ whose solution reads

$$|\tilde{\psi}_c(t)| = |\tilde{\psi}_c(0)| \frac{e^{-\tilde{r}_I t}}{\sqrt{1 + \alpha (1 - e^{-2\tilde{r}_I t})}} \quad (13)$$

$$\tilde{\theta}(t) = -\tilde{r}_R t + \tilde{u}_R \int_0^t dt' |\tilde{\psi}_c(t')|^2 \quad (14)$$

where $\tilde{r}, \tilde{u}$ are complex coefficients given by

$$\tilde{r} \equiv \tilde{r}_R + i\tilde{r}_I = r/K_1 \quad (15)$$

$$\tilde{u} \equiv \tilde{u}_R + i\tilde{u}_I = u/K_1 \quad (16)$$

while $\alpha = |\tilde{\psi}_c(0)||\tilde{u}_I|/\tilde{r}_I$. This solution describes a dynamical transition at a critical point $J_c$. Specifically for $J < J_c$ the order parameter shows damped oscillation toward zero

$$\tilde{\psi}_c(t) \sim e^{-t/\tau} e^{-i\delta\Omega t} \quad (17)$$

with a divergent relaxation time $\tau \sim 1/\tilde{r}_I \sim 1/(J_c - J)$ and an oscillation frequency going to zero $\delta\Omega \sim \tilde{r}_R \sim (J_c - J)$. Upon crossing the critical point, for $J > J_c$, the dynamics shows instead an amplification of the order parameter which saturates at long time into a train of finite amplitude oscillations

$$\tilde{\psi}_c(t) \sim |\tilde{\psi}_c(\infty)| e^{-i\delta\Omega_{sf}t} \quad (18)$$

with $|\tilde{\psi}_c(\infty)| = \sqrt{|\tilde{r}_I|/|\tilde{u}_I|} \sim \sqrt{J - J_c}$ and $\delta\Omega_{sf} \sim (\tilde{r}_R + \tilde{r}_I\tilde{u}_R/\tilde{u}_I) \sim (J - J_c)$. In the normal phase, where the non-linearity $u$ in Eq. (12) is essentially irrelevant, the transient dynamics shows harmonic oscillations while in the broken symmetry phase multiple frequencies are present, at least on intermediate time scales, as encoded in the phase dynamics (14). Right at the transition, for $J = J_c$ when $\tilde{r}_{R,I} \to 0$, the amplitude of the order parameter decays towards zero as a power-law$^{39,40}$ while the angular velocity vanishes and the phase grows in time only logarithmically,

$$|\tilde{\psi}_c(t)| \sim 1/\sqrt{t} \quad (19)$$

$$\tilde{\theta}(t) \sim \log(1 + 2\tilde{u}_I |\tilde{\psi}_c(0)|t) \quad (20)$$

## B. Discussion

In the previous section we have shown that disregarding (i) retardation effects in the effective action, i.e. expanding all local correlators around the critical frequency $\Omega_*$, as well as (ii) disregarding terms higher than quadratic in the quantum fields allow to fully reproduce the results obtained by time-dependent Gutzwiller decoupling, describing the finite-frequency dynamical transition at the mean field level. Still the full effective action in Eq. (11) includes the effect of non linearities, noise and quantum fluctuations beyond this semiclassical mean field dynamics and can result in a renormalization of critical behavior and non-mean field exponents$^{38,41,42}$. These can be captured with a renormalization group treatment of the finite-frequency criticality, along the lines discussed for the equilibrium Bose Hubbard model$^{43}$ as well as for weakly interacting non-equilibrium superfluids$^{12,44}$. Particularly interesting in this respect is the role under renormalization of higher order expansion coefficients in the effective action, such as $K_1, K_2$ and the next order noise term controlled by $\partial_\omega G_{loc}^K(\Omega_*)$. The former are related to emergent symmetries, such as particle-hole in the equilibrium Bose-Hubbard model$^2$ or the asymptotic equilibrium symmetry in the driven-dissipative condensation, whose deviation results in a KPZ like critical phase dynamics$^{38}$. A term analogous to the latter was shown instead to give rise to non-trivial critical behavior in diffusively driven one dimensional bosons$^{14}$. An interesting question is whether the inclusion of fluctuations beyond the gaussian level could completely wash away the finite-frequency transition or renormalize the critical frequency $\Omega_*$ down to zero, resulting in a static transition. While
term $S_{noise} = D \int dxdt \tilde{\psi}_q^*(x, t) \tilde{\psi}_q(x, t)$. Here $D$ plays the role of effective diffusion coefficient in the equivalent stochastic (Langevin) dynamics and it is indeed given by $D \simeq G_{loc}^K(\Omega_*) \sim T_{eff}$. Finally $S_{int}$ accounts for the non-linearities and it is completely determined by the multi-particle Green's functions of the driven-dissipative single site problem. If we restrict ourselves to interaction terms with one quantum and three classical fields, which is valid in high enough dimensions according to canonical power counting$^{38}$, we can write this term as $S_{int} = u \int dxdt \tilde{\psi}_q^*(x, t) \tilde{\psi}_c^*(x, t) \tilde{\psi}_c^2(x, t) + hc$. We can now take the saddle point equation $\delta S/\delta \tilde{\psi}_q^*(x, t) = 0$ and obtain the equation of motion

$$(iK_1\partial_t - K_3\nabla^2 - r) \tilde{\psi}_c + u|\tilde{\psi}_c|^2 \tilde{\psi}_c = 0 \quad (12)$$

which takes the form of a complex Ginzburg-Landau equation, well known as a phenomenological description of pattern formation in classical non-equilibrium systems$^{3-5}$. The spatially homogeneous solution of Equation (12) can be obtained in closed form, as we discuss in the next section, and it describes a transition between a phase where $\tilde{\psi}_c \to 0$ for $t \to \infty$ and a phase where the modulus of the order parameter saturates to a finite value. It further allows us to compute the scales $1/\tau$ and $\delta\Omega$ which perfectly match the numerical results that we have found from the Gutzwiller dynamics in the previous section, as shown in Figure 5. We remark that this classical dynamics only describes the evolution in a frame rotating at frequency $\Omega_*$.

## A. Analytical Solution of Saddle Point Dynamics

We discuss here more in detail the solution of Equation (12). Assuming an homogeneous solution and introducing polar coordinates for absolute value and phase of the order parameters $\tilde{\psi}_c(t) = |\tilde{\psi}_c(t)|e^{i\theta(t)}$ one obtains two separate equations for $|\tilde{\psi}_c(t)|$ and $\tilde{\theta}(t)$ whose solution reads

$$|\tilde{\psi}_c(t)| = |\tilde{\psi}_c(0)| \frac{e^{-\tilde{r}_I t}}{\sqrt{1 + \alpha (1 - e^{-2\tilde{r}_I t})}} \quad (13)$$

$$\tilde{\theta}(t) = -\tilde{r}_R t + \tilde{u}_R \int_0^t dt' |\tilde{\psi}_c(t')|^2 \quad (14)$$

where $\tilde{r}, \tilde{u}$ are complex coefficients given by

$$\tilde{r} \equiv \tilde{r}_R + i\tilde{r}_I = r/K_1 \quad (15)$$

$$\tilde{u} \equiv \tilde{u}_R + i\tilde{u}_I = u/K_1 \quad (16)$$

while $\alpha = |\tilde{\psi}_c(0)| |\tilde{u}_I| / \tilde{r}_I$. This solution describes a dynamical transition at a critical point $J_c$. Specifically for $J < J_c$ the order parameter shows damped oscillation toward zero

$$\tilde{\psi}_c(t) \sim e^{-t/\tau} e^{-i\delta\Omega t} \quad (17)$$

with a divergent relaxation time $\tau \sim 1/\tilde{r}_I \sim 1/(J_c - J)$ and an oscillation frequency going to zero $\delta\Omega \sim \tilde{r}_R \sim (J_c - J)$. Upon crossing the critical point, for $J > J_c$, the dynamics shows instead an amplification of the order parameter which saturates at long time into a train of finite amplitude oscillations

$$\tilde{\psi}_c(t) \sim |\tilde{\psi}_c(\infty)| e^{-i\delta\Omega_{sf}t} \quad (18)$$

with $|\tilde{\psi}_c(\infty)| = \sqrt{|\tilde{r}_I|/|\tilde{u}_I|} \sim \sqrt{J - J_c}$ and $\delta\Omega_{sf} \sim (\tilde{r}_R + \tilde{r}_I \tilde{u}_R / \tilde{u}_I) \sim (J - J_c)$. In the normal phase, where the non-linearity $u$ in Eq. (12) is essentially irrelevant, the transient dynamics shows harmonic oscillations while in the broken symmetry phase multiple frequencies are present, at least on intermediate time scales, as encoded in the phase dynamics (14). Right at the transition, for $J = J_c$ when $\tilde{r}_{R,I} \to 0$, the amplitude of the order parameter decays towards zero as a power-law$^{39,40}$ while the angular velocity vanishes and the phase grows in time only logarithmically,

$$|\tilde{\psi}_c(t)| \sim 1/\sqrt{t} \quad (19)$$

$$\tilde{\theta}(t) \sim \log(1 + 2\tilde{u}_I |\tilde{\psi}_c(0)|t) \quad (20)$$

## B. Discussion

In the previous section we have shown that disregarding (i) retardation effects in the effective action, i.e. expanding all local correlators around the critical frequency $\Omega_*$, as well as (ii) disregarding terms higher than quadratic in the quantum fields allow to fully reproduce the results obtained by time-dependent Gutzwiller decoupling, describing the finite-frequency dynamical transition at the mean field level. Still the full effective action in Eq. (11) includes the effect of non linearities, noise and quantum fluctuations beyond this semiclassical mean field dynamics and can result in a renormalization of critical behavior and non-mean field exponents$^{38,41,42}$. These can be captured with a renormalization group treatment of the finite-frequency criticality, along the lines discussed for the equilibrium Bose Hubbard model$^{43}$ as well as for weakly interacting non-equilibrium superfluids$^{12,44}$. Particularly interesting in this respect is the role under renormalization of higher order expansion coefficients in the effective action, such as $K_1, K_2$ and the next order noise term controlled by $\partial_\omega G_{loc}^K(\Omega_*)$. The former are related to emergent symmetries, such as particle-hole in the equilibrium Bose-Hubbard model$^2$ or the asymptotic equilibrium symmetry in the driven-dissipative condensation, whose deviation results in a KPZ like critical phase dynamics$^{38}$. A term analogous to the latter was shown instead to give rise to non-trivial critical behavior in diffusively driven one dimensional bosons$^{14}$. An interesting question is whether the inclusion of fluctuations beyond the gaussian level could completely wash away the finite-frequency transition or renormalize the critical frequency $\Omega_*$ down to zero, resulting in a static transition. While
answering this question certainly deserves further investigation our results suggest that, provided the effective action (11) in the rotating frame admits a non-vanishing $U(1)$ order parameter $\tilde{\psi}_c \neq 0$, which is expected in high enough dimensions, then the broken symmetry phase in the original frame will display undamped oscillations and breaking of time-translational invariance.

## VI. CONCLUSIONS

In this work we have shown that a prototype model of correlated driven-dissipative lattice bosons develops, for a critical value of the hopping rate, a diverging susceptibility at a non-zero frequency $\Omega_*$ . The resulting finite-frequency criticality corresponds to the dissipative dynamics lacking of a stationary state and rather oscillating in time without damping. Writing down the effective Keldysh field theory for this finite frequency transition we have obtained its semiclassical limit which we show to reproduce the results of a time-dependent Gutzwiller decoupling of the density matrix. We emphasize that knowing the critical frequency $\Omega_*$ requires the quantum solution of the single-site dissipative interacting problem and it is therefore not contained in the semi-classical equation of motion which only describes the dynamics in the frame rotating at $\Omega_*$. Our results differ from other studies of limit cycles instabilities in driven-dissipative systems, such as exciton-polariton condensates described by Gross-Pitaevski (GP) types of equation and it could be seen as the strongly correlated version of them. Indeed our transition shares genuine features of a dissipative Mott-Superfluid quantum phase transitions being tuned both by coherent couplings and pump/loss rates. In particular our incoherent phase exists at small hopping even beyond the standard threshold of pump greater than losses, an effect which is genuinely quantum mechanical due to the Hubbard repulsion favouring Fock-like states rather than coherent states. Furthermore the frequency of the limit cycle is set by the local occupation rather than the local coherence as in GP theories.

Our work suggests several interesting future directions. From one side it would be interesting to include dynamical and spatial fluctuations on top of the semiclassical dynamics for the order parameter and study the fate of this dissipative dynamical transition in finite dimensions, following similar investigations done for dynamical transitions in isolated quantum systems$^{45–47}$. Another intriguing open question is whether a similar finite-frequency criticality exists in models of driven-dissipative systems with discrete broken symmetry phases$^{22}$ or even in presence of a purely coherent drive, as for example in the context of optomechanical platforms$^{48,49}$ or coherently driven quantum spin chains$^{20}$.

Finally, while our work focuses on a paradigmatic model of driven-dissipative bosons which is relevant for the upcoming generation of circuit QED arrays experiments$^{25,26}$, it also outlines a generic framework to study

dynamical instabilities in non-equilibrium quantum systems, by focusing on frequency dependent response functions and their divergences. Such a framework has the potential to be applied in a wide range of contexts, including for example driven and isolated Floquet systems, where breaking of discrete time-translational symmetry has been predicted$^{50–52}$ and observed$^{53,54}$, quantum systems undergoing various forms of synchronization$^{48,55,56}$ as well as electronic systems under pump-probe optical-irradiation$^{57}$

*Acknowledgements.* We acknowledge discussions with A. Clerk, M. Goldstein, V. Savona. This work was supported by the CNRS through the PICS-USA-147504, by a grant "Investissements d'Avenir" from LabEx PALM (ANR-10-LABX-0039-PALM) and by a grant IRS-IQUPS of University Paris-Saclay.

## Appendix A: Role of Driving Protocol

The results we have described concerning the finite frequency criticality are extremely robust with respect to the specific driving protocol, provided that bosons are injected in the lattice incoherently. There are however interesting differences in the nature of the normal phase which strongly depend on the nature of the drive, as we are going to discuss below. In this work we have considered two driving protocols which we discuss more in detail here and we address the main qualitative differences in the results for the two cases.

In the main text, we considered a scheme of incoherent pumping, which has been recently proposed$^{24,31,32}$, arising from an ensemble of $N_{at}$ driven two-level emitters embedded in each cavity and having randomly distributed frequencies. In the following we will refer to this scheme as cold drive, for reasons which will become clear in the next sections. The microscopic Hamiltonian for this driving scheme reads

$$H_{\text{pump}}^{\text{cold}} = \sum_{i} \sum_{n=1}^{N_{at}} \omega_{at}^{(n)} \sigma_{i}^{+(n)} \sigma_{i}^{-(n)} + g \sum_{i,n} (a_{i}^{\dagger} \sigma_{i}^{-(n)} + hc) \quad (\text{A1})$$

where the transition frequencies $\omega_{at}^{(n)}$ of the two-level systems are assumed to be uniformly distributed over a finite range and each emitter is incoherently pumped in the excited state$^{31}$.

In this appendix we consider also a second scheme$^{30}$, which we will refer to as hot drive in the following, where a random classical drive is modulated with a coherent tone as described by the time-dependent Hamiltonian

$$H_{\text{pump}}^{\text{hot}}(t) = \sum_{i} \left( e^{i\omega_{L}t} a_{i}^{\dagger}\eta(t) + hc \right) \quad (\text{A2})$$

where $\eta(t)$ is assumed to have gaussian statistics with zero average, $\langle\eta(t)\rangle = 0$, and correlations $\langle\eta(t)\eta(t')\rangle = fC_{\sigma}(t-t')$. We assume the noise spectrum to be box-shaped with a finite bandwidth $\sigma$, i.e. $C_{\sigma}(\omega) = \theta(\sigma-|\omega|)$,
answering this question certainly deserves further investigation our results suggest that, provided the effective action (11) in the rotating frame admits a non-vanishing $U(1)$ order parameter $\tilde{\psi}_c \neq 0$, which is expected in high enough dimensions, then the broken symmetry phase in the original frame will display undamped oscillations and breaking of time-translational invariance.

## VI. CONCLUSIONS

In this work we have shown that a prototype model of correlated driven-dissipative lattice bosons develops, for a critical value of the hopping rate, a diverging susceptibility at a non-zero frequency $\Omega_*$ . The resulting finite-frequency criticality corresponds to the dissipative dynamics lacking of a stationary state and rather oscillating in time without damping. Writing down the effective Keldysh field theory for this finite frequency transition we have obtained its semiclassical limit which we show to reproduce the results of a time-dependent Gutzwiller decoupling of the density matrix. We emphasize that knowing the critical frequency $\Omega_*$ requires the quantum solution of the single-site dissipative interacting problem and it is therefore not contained in the semi-classical equation of motion which only describes the dynamics in the frame rotating at $\Omega_*$. Our results differ from other studies of limit cycles instabilities in driven-dissipative systems, such as exciton-polariton condensates described by Gross-Pitaevski (GP) types of equation and it could be seen as the strongly correlated version of them. Indeed our transition shares genuine features of a dissipative Mott-Superfluid quantum phase transitions being tuned both by coherent couplings and pump/loss rates. In particular our incoherent phase exists at small hopping even beyond the standard threshold of pump greater than losses, an effect which is genuinely quantum mechanical due to the Hubbard repulsion favouring Fock-like states rather than coherent states. Furthermore the frequency of the limit cycle is set by the local occupation rather than the local coherence as in GP theories.

Our work suggests several interesting future directions. From one side it would be interesting to include dynamical and spatial fluctuations on top of the semiclassical dynamics for the order parameter and study the fate of this dissipative dynamical transition in finite dimensions, following similar investigations done for dynamical transitions in isolated quantum systems$^{45–47}$. Another intriguing open question is whether a similar finite-frequency criticality exists in models of driven-dissipative systems with discrete broken symmetry phases$^{22}$ or even in presence of a purely coherent drive, as for example in the context of optomechanical platforms$^{48,49}$ or coherently driven quantum spin chains$^{20}$.

Finally, while our work focuses on a paradigmatic model of driven-dissipative bosons which is relevant for the upcoming generation of circuit QED arrays experiments$^{25,26}$, it also outlines a generic framework to study

dynamical instabilities in non-equilibrium quantum systems, by focusing on frequency dependent response functions and their divergences. Such a framework has the potential to be applied in a wide range of contexts, including for example driven and isolated Floquet systems, where breaking of discrete time-translational symmetry has been predicted$^{50–52}$ and observed$^{53,54}$, quantum systems undergoing various forms of synchronization$^{48,55,56}$ as well as electronic systems under pump-probe optical-irradiation$^{57}$

*Acknowledgements.* We acknowledge discussions with A. Clerk, M. Goldstein, V. Savona. This work was supported by the CNRS through the PICS-USA-147504, by a grant "Investissements d'Avenir" from LabEx PALM (ANR-10-LABX-0039-PALM) and by a grant IRS-IQUPS of University Paris-Saclay.

## Appendix A: Role of Driving Protocol

The results we have described concerning the finite frequency criticality are extremely robust with respect to the specific driving protocol, provided that bosons are injected in the lattice incoherently. There are however interesting differences in the nature of the normal phase which strongly depend on the nature of the drive, as we are going to discuss below. In this work we have considered two driving protocols which we discuss more in detail here and we address the main qualitative differences in the results for the two cases.

In the main text, we considered a scheme of incoherent pumping, which has been recently proposed$^{24,31,32}$, arising from an ensemble of $N_{at}$ driven two-level emitters embedded in each cavity and having randomly distributed frequencies. In the following we will refer to this scheme as cold drive, for reasons which will become clear in the next sections. The microscopic Hamiltonian for this driving scheme reads

$$H_{\text{pump}}^{\text{cold}} = \sum_{i} \sum_{n=1}^{N_{at}} \omega_{at}^{(n)} \sigma_{i}^{+(n)} \sigma_{i}^{-(n)} + g \sum_{i,n} (a_{i}^{\dagger} \sigma_{i}^{-(n)} + hc) \quad (\text{A1})$$

where the transition frequencies $\omega_{at}^{(n)}$ of the two-level systems are assumed to be uniformly distributed over a finite range and each emitter is incoherently pumped in the excited state$^{31}$.

In this appendix we consider also a second scheme$^{30}$, which we will refer to as hot drive in the following, where a random classical drive is modulated with a coherent tone as described by the time-dependent Hamiltonian

$$H_{\text{pump}}^{\text{hot}}(t) = \sum_{i} \left( e^{i\omega_{L}t} a_{i}^{\dagger}\eta(t) + hc \right) \quad (\text{A2})$$

where $\eta(t)$ is assumed to have gaussian statistics with zero average, $\langle\eta(t)\rangle = 0$, and correlations $\langle\eta(t)\eta(t')\rangle = fC_{\sigma}(t-t')$. We assume the noise spectrum to be box-shaped with a finite bandwidth $\sigma$, i.e. $C_{\sigma}(\omega) = \theta(\sigma-|\omega|)$,
<figure><img src="image_8.png" /><figcaption>FIG. 7. Effective temperature $T_{eff}$ as a function of the drive bandwidth $\sigma$ for the two driving protocols considered. We notice that while in the cold drive case $T_{eff}$ depends only weakly on $\sigma$ and stay small, for the hot drive case it substantially increases as $\sigma$ (and therefore the number of photons) increases. Parameters: loss rate $\kappa = 0.0128U$, resonator frequency $\delta\omega_0 = 0.0$, drive bandwidth (left panel) $\sigma = 0.625U$, drive amplitude (right panel) $f = 0.0625U$.</figcaption></figure>

and amplitude $f$, although the results we obtain do not depend qualitatively from the exact shape of $C_{\sigma}(\omega)$. Treating the incoherent driving at the master-equation level we obtain, in both cases of Eqs. (A2-A1), the contribution to the dissipator in Eq. (4) which we report here for clarity$^{30,31}$

$$ \mathcal{D}_{\text{pump}}[\rho] = \sum_{i} f_{\text{in}} \tilde{D}[a_{i}^{\dagger}, \tilde{a}_{i\sigma}; \rho] + f_{\text{out}} \tilde{D}[a_{i}, \tilde{a}_{i\sigma}^{\dagger}; \rho] \quad (\text{A3}) $$

where we have already introduced the modified dissipator $\tilde{D}$ and the photon operator dressed by the finite bandwidth drive $\tilde{a}_{\sigma}^{\dagger}$ in section II. In the random noise case, Eq. (A2), we obtain $f_{\text{in}} = f_{\text{out}} \equiv f$, namely the drive acts both as a source and as a sink of particles, much like a finite-temperature bath. Instead, in the case of inverted random emitters we have $f_{\text{in}} = f$ and $f_{\text{out}} = 0$, namely there are no additional losses of particles associated to the drive.

## 1. Local Bosonic Occupation and Effective Temperature

As we discussed in section II, figure 1, the boson number in the single-site problem as a function of pump bandwidth $\sigma$, shows a starcaise structure characteristic of blockade physics$^{30}$, with a value of $\Delta\sigma \sim U$ required to add extra bosons in the system. In this respect, an important difference between the two driving protocols already appear, namely the cold drive is able to fix the occupancy to almost integer filling$^{31}$, while the hot drive to half-integer filling, reflecting the fact that the stationary density matrix is almost pure in the cold drive case, while it has a box-shaped distribution of populations in the hot

<figure><img src="image_9.png" /><figcaption>FIG. 8. Stationary State Phase diagram as a function of hopping strength $J$ and drive bandwidth $\sigma$, for the two driving protocols considered namely cold drive (top panel) and hot drive (bottom panel). Parameters: loss rate $\kappa = 0.0128U$, resonator frequency $\delta\omega_0 = 0.0$, drive amplitude (left panel) $f = 0.0625U$.</figcaption></figure>

case. The dependence of $T_{eff}$, introduced in section II A, from the drive bandwidth $\sigma$ is plotted in figure 7 and reveals a rather substantial difference in the two driving protocols for what concerns the effective heating properties of the system. Indeed upon increasing the drive bandwidth the effective temperature increases in the hot drive case while decreases (slightly) in the cold case one. For the hot drive, this can be understood because the density matrix would reach an infinite temperature one by sending the bandwidth of the drive to infinity. This offers an alternative perspective on the recent proposed scheme to engineer effective ground state phases of interacting photons through the use of non-markovian reservoirs$^{31,32}$.

## 2. Phase Diagram

Finally we conclude presenting, in figure 8, the stationary phase diagram of the finite-frequency phase transition for the two driving protocols we have discussed so far. In the main text we considered a specific value of the drive bandwidth $\sigma$ ($\sigma = 1.5U$) while now we present the phase boundary in the $\sigma, J$ plane. As we discussed in the main text $\sigma$ controls both the local density and the effective chemical potential $\Omega_*$, in step like fashion that resembles the equilibrium ground-state physics of the problem. It is therefore a natural choice for a tuning parameter in the phase diagram. For both driving protocols we generically find a similar behavior, namely a small hopping phase $J < J_c(\sigma)$ which has a stable stationary state fully incoherent and a large hopping regime $J > J_c(\sigma)$ where the stationary state becomes unstable toward an oscillating regime and the system develops a $U(1)$ order parameter at finite frequency. It is neverthe-
<figure><img src="image_9.png" /><figcaption>FIG. 7. Effective temperature $T_{eff}$ as a function of the drive bandwidth $\sigma$ for the two driving protocols considered. We notice that while in the cold drive case $T_{eff}$ depends only weakly on $\sigma$ and stay small, for the hot drive case it substantially increases as $\sigma$ (and therefore the number of photons) increases. Parameters: loss rate $\kappa = 0.0128U$, resonator frequency $\delta\omega_0 = 0.0$, drive bandwidth (left panel) $\sigma = 0.625U$, drive amplitude (right panel) $f = 0.0625U$.</figcaption></figure>

and amplitude $f$, although the results we obtain do not depend qualitatively from the exact shape of $C_{\sigma}(\omega)$. Treating the incoherent driving at the master-equation level we obtain, in both cases of Eqs. (A2-A1), the contribution to the dissipator in Eq. (4) which we report here for clarity$^{30,31}$

$$ \mathcal{D}_{\text{pump}}[\rho] = \sum_{i} f_{\text{in}} \tilde{D}[a_{i}^{\dagger}, \tilde{a}_{i\sigma}; \rho] + f_{\text{out}} \tilde{D}[a_{i}, \tilde{a}_{i\sigma}^{\dagger}; \rho] \quad (\text{A3}) $$

where we have already introduced the modified dissipator $\tilde{D}$ and the photon operator dressed by the finite bandwidth drive $\tilde{a}_{\sigma}^{\dagger}$ in section II. In the random noise case, Eq. (A2), we obtain $f_{\text{in}} = f_{\text{out}} \equiv f$, namely the drive acts both as a source and as a sink of particles, much like a finite-temperature bath. Instead, in the case of inverted random emitters we have $f_{\text{in}} = f$ and $f_{\text{out}} = 0$, namely there are no additional losses of particles associated to the drive.

## 1. Local Bosonic Occupation and Effective Temperature

As we discussed in section II, figure 1, the boson number in the single-site problem as a function of pump bandwidth $\sigma$, shows a starcaise structure characteristic of blockade physics$^{30}$, with a value of $\Delta\sigma \sim U$ required to add extra bosons in the system. In this respect, an important difference between the two driving protocols already appear, namely the cold drive is able to fix the occupancy to almost integer filling$^{31}$, while the hot drive to half-integer filling, reflecting the fact that the stationary density matrix is almost pure in the cold drive case, while it has a box-shaped distribution of populations in the hot

<figure><img src="image_10.png" /><figcaption>FIG. 8. Stationary State Phase diagram as a function of hopping strength $J$ and drive bandwidth $\sigma$, for the two driving protocols considered namely cold drive (top panel) and hot drive (bottom panel). Parameters: loss rate $\kappa = 0.0128U$, resonator frequency $\delta\omega_0 = 0.0$, drive amplitude (left panel) $f = 0.0625U$.</figcaption></figure>

case. The dependence of $T_{eff}$, introduced in section II A, from the drive bandwidth $\sigma$ is plotted in figure 7 and reveals a rather substantial difference in the two driving protocols for what concerns the effective heating properties of the system. Indeed upon increasing the drive bandwidth the effective temperature increases in the hot drive case while decreases (slightly) in the cold case one. For the hot drive, this can be understood because the density matrix would reach an infinite temperature one by sending the bandwidth of the drive to infinity. This offers an alternative perspective on the recent proposed scheme to engineer effective ground state phases of interacting photons through the use of non-markovian reservoirs$^{31,32}$.

## 2. Phase Diagram

Finally we conclude presenting, in figure 8, the stationary phase diagram of the finite-frequency phase transition for the two driving protocols we have discussed so far. In the main text we considered a specific value of the drive bandwidth $\sigma$ ($\sigma = 1.5U$) while now we present the phase boundary in the $\sigma, J$ plane. As we discussed in the main text $\sigma$ controls both the local density and the effective chemical potential $\Omega_*$, in step like fashion that resembles the equilibrium ground-state physics of the problem. It is therefore a natural choice for a tuning parameter in the phase diagram. For both driving protocols we generically find a similar behavior, namely a small hopping phase $J < J_c(\sigma)$ which has a stable stationary state fully incoherent and a large hopping regime $J > J_c(\sigma)$ where the stationary state becomes unstable toward an oscillating regime and the system develops a $U(1)$ order parameter at finite frequency. It is neverthe-
less quite interesting to discuss the different shapes of the
phase boundary, which instead rather strongly depend on
the protocol. We notice that for the cold drive case(top
panel) the boundary resembles the ground state one, with
a lobe-like structure for different values of the local den-
sity and a critical hopping $J_c$ which decreases as the local
filling increases. Viceversa, in the hot drive case (bottom
panel) we find a rather opposite effect, namely the critical
hopping increases with $\sigma$ and the region of normal phase

stability expands. We can understand this effect from the
discussion on the occupation of the bosonic mode and the
effective temperature: indeed in the hot drive case, in-
creasing the bandwidth $\sigma$ has the effect of both changing
the local occupation (see figure 1) and of increasing the
effective temperature (see figure A), with the result of
shrinking the broken symmetry region due to effectively
increased thermal fluctuations.

¹ L. D. Landau and E. M. Lifshitz, *Statistical Physics, Part I* (1995).

² S. Sachdev, *Quantum Phase Transitions* (Cambridge University Press, 1998).

³ M. Cross and P. Hohenberg, Rev. Mod. Phys. **65**, 851 (1993).

⁴ M. Cross and H. Greenside, *Pattern Formation and Dynamics in Nonequilibrium Systems*, 1st ed. (Cambridge University Press, United Kingdom, 2009).

⁵ W. van Saarloos, Physics Reports **386**, 29 (2003).

⁶ I. Bloch, J. Dalibard, and S. Nascimbène, Nature Physics **8**, 267 EP (2012).

⁷ R. Blatt and C. F. Roos, Nature Physics **8**, 277 EP (2012).

⁸ A. Wallraff, D. Schuster, A. Blais, L. Frunzio, M. J. Huang, R.S., S. Kumar, S. M. Girvin, and R. J. Schoelkopf, Nature **431**, 162 (2004).

⁹ A. A. Houck, H. E. Tureci, and J. Koch, Nature Physics **8** (2012).

¹⁰ S. Schmidt and J. Koch, Annalen der Physik **525**, 395 (2013).

¹¹ K. L. Hur, L. Henriet, A. Petrescu, K. Plekhanov, G. Roux, and M. Schiró, Comptes Rendus Physique **17**, 808 (2016).

¹² L. M. Sieberer, S. D. Huber, E. Altman, and S. Diehl, Phys. Rev. Lett. **110**, 195301 (2013).

¹³ M. F. Maghrebi and A. V. Gorshkov, Phys. Rev. B **93**, 014307 (2016).

¹⁴ J. Marino and S. Diehl, Phys. Rev. Lett. **116**, 070407 (2016).

¹⁵ J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymanska, R. Andre, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, Nature **443**, 409 (2006).

¹⁶ M. H. Szymańska, J. Keeling, and P. B. Littlewood, Phys. Rev. Lett. **96**, 230602 (2006).

¹⁷ H. Deng, H. Haug, and Y. Yamamoto, Rev. Mod. Phys. **82**, 1489 (2010).

¹⁸ I. Carusotto and C. Ciuti, Rev. Mod. Phys. **85**, 299 (2013).

¹⁹ J. Jin, D. Rossini, R. Fazio, M. Leib, and M. J. Hartmann, Phys. Rev. Lett. **110**, 163605 (2013).

²⁰ C.-K. Chan, T. E. Lee, and S. Gopalakrishnan, Phys. Rev. A **91**, 051601 (2015).

²¹ R. M. Wilson, K. W. Mahmud, A. Hu, A. V. Gorshkov, M. Hafezi, and M. Foss-Feig, Phys. Rev. A **94**, 033801 (2016).

²² M. Schiró, C. Joshi, M. Bordyuh, R. Fazio, J. Keeling, and H. E. Türeci, Phys. Rev. Lett. **116**, 143603 (2016).

²³ E. T. Owen, J. Jin, D. Rossini, R. Fazio, and M. J. Hartmann, ArXiv e-prints (2017), arXiv:1711.11309 [quant-ph].

²⁴ A. Biella, F. Storme, J. Lebreuilly, D. Rossini, R. Fazio, I. Carusotto, and C. Ciuti, Phys. Rev. A **96**, 023839 (2017).

²⁵ M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch, and A. A. Houck, Phys. Rev. X **7**, 011016 (2017).

²⁶ R. Ma, C. Owens, A. Houck, D. I. Schuster, and J. Simon, Phys. Rev. A **95**, 043811 (2017).

²⁷ R. Ma, B. Saxberg, C. Owens, N. Leung, Y. Lu, J. Simon, and D. I. Schuster, ArXiv e-prints (2018), arXiv:1807.11342 [cond-mat.quant-gas].

²⁸ J. Raftery, D. Sadri, S. Schmidt, H. E. Türeci, and A. A. Houck, Phys. Rev. X **4**, 031043 (2014).

²⁹ S. Hacohen-Gourgy, V. V. Ramasesh, C. De Grandi, I. Siddiqi, and S. M. Girvin, Phys. Rev. Lett. **115**, 240501 (2015).

³⁰ A. J. Hoffman, S. J. Srinivasan, S. Schmidt, L. Spietz, J. Aumentado, H. E. Türeci, and A. A. Houck, Phys. Rev. Lett. **107**, 053602 (2011).

³¹ J. Lebreuilly, A. Biella, F. Storme, D. Rossini, R. Fazio, C. Ciuti, and I. Carusotto, Phys. Rev. A **96**, 033828 (2017).

³² J. Lebreuilly, M. Wouters, and I. Carusotto, Comptes Rendus Physique **17**, 836 (2016), polariton physics / Physique des polaritons.

³³ M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B **40**, 546 (1989).

³⁴ O. Scarlatella, A. A. Clerk, and M. Schirò, arXiv e-prints, arXiv:1811.03518 (2018), arXiv:1811.03518 [quant-ph].

³⁵ S. Schmidt and G. Blatter, Phys. Rev. Lett. **103**, 086403 (2009).

³⁶ J. Koch and K. Le Hur, Phys. Rev. A **80**, 023811 (2009).

³⁷ K. Dunnett and M. H. Szymańska, Phys. Rev. B **93**, 195306 (2016).

³⁸ L. M. Sieberer, M. Buchhold, and S. Diehl, Reports on Progress in Physics **79**, 096001 (2016).

³⁹ S. Diehl, A. Tomadin, A. Micheli, R. Fazio, and P. Zoller, Phys. Rev. Lett. **105**, 015702 (2010).

⁴⁰ A. Tomadin, S. Diehl, and P. Zoller, Phys. Rev. A **83**, 013611 (2011).

⁴¹ C. De Dominicis, E. Brézin, and J. Zinn-Justin, Phys. Rev. B **12**, 4945 (1975).

⁴² P. Hohenberg and B. Halperin, Rev. Mod. Phys. **49**, 435 (1977).

⁴³ A. Rançon and N. Dupuis, Phys. Rev. B **84**, 174513 (2011).

⁴⁴ J. Marino and S. Diehl, Phys. Rev. B **94**, 085150 (2016).
less quite interesting to discuss the different shapes of the
phase boundary, which instead rather strongly depend on
the protocol. We notice that for the cold drive case(top
panel) the boundary resembles the ground state one, with
a lobe-like structure for different values of the local den-
sity and a critical hopping $J_c$ which decreases as the local
filling increases. Viceversa, in the hot drive case (bottom
panel) we find a rather opposite effect, namely the critical
hopping increases with $\sigma$ and the region of normal phase

stability expands. We can understand this effect from the
discussion on the occupation of the bosonic mode and the
effective temperature: indeed in the hot drive case, in-
creasing the bandwidth $\sigma$ has the effect of both changing
the local occupation (see figure 1) and of increasing the
effective temperature (see figure A), with the result of
shrinking the broken symmetry region due to effectively
increased thermal fluctuations.

¹ L. D. Landau and E. M. Lifshitz, *Statistical Physics, Part I* (1995).

² S. Sachdev, *Quantum Phase Transitions* (Cambridge University Press, 1998).

³ M. Cross and P. Hohenberg, Rev. Mod. Phys. **65**, 851 (1993).

⁴ M. Cross and H. Greenside, *Pattern Formation and Dynamics in Nonequilibrium Systems*, 1st ed. (Cambridge University Press, United Kingdom, 2009).

⁵ W. van Saarloos, Physics Reports **386**, 29 (2003).

⁶ I. Bloch, J. Dalibard, and S. Nascimbène, Nature Physics **8**, 267 EP (2012).

⁷ R. Blatt and C. F. Roos, Nature Physics **8**, 277 EP (2012).

⁸ A. Wallraff, D. Schuster, A. Blais, L. Frunzio, M. J. Huang, R.S., S. Kumar, S. M. Girvin, and R. J. Schoelkopf, Nature **431**, 162 (2004).

⁹ A. A. Houck, H. E. Tureci, and J. Koch, Nature Physics **8** (2012).

¹⁰ S. Schmidt and J. Koch, Annalen der Physik **525**, 395 (2013).

¹¹ K. L. Hur, L. Henriet, A. Petrescu, K. Plekhanov, G. Roux, and M. Schiró, Comptes Rendus Physique **17**, 808 (2016).

¹² L. M. Sieberer, S. D. Huber, E. Altman, and S. Diehl, Phys. Rev. Lett. **110**, 195301 (2013).

¹³ M. F. Maghrebi and A. V. Gorshkov, Phys. Rev. B **93**, 014307 (2016).

¹⁴ J. Marino and S. Diehl, Phys. Rev. Lett. **116**, 070407 (2016).

¹⁵ J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymanska, R. Andre, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, Nature **443**, 409 (2006).

¹⁶ M. H. Szymańska, J. Keeling, and P. B. Littlewood, Phys. Rev. Lett. **96**, 230602 (2006).

¹⁷ H. Deng, H. Haug, and Y. Yamamoto, Rev. Mod. Phys. **82**, 1489 (2010).

¹⁸ I. Carusotto and C. Ciuti, Rev. Mod. Phys. **85**, 299 (2013).

¹⁹ J. Jin, D. Rossini, R. Fazio, M. Leib, and M. J. Hartmann, Phys. Rev. Lett. **110**, 163605 (2013).

²⁰ C.-K. Chan, T. E. Lee, and S. Gopalakrishnan, Phys. Rev. A **91**, 051601 (2015).

²¹ R. M. Wilson, K. W. Mahmud, A. Hu, A. V. Gorshkov, M. Hafezi, and M. Foss-Feig, Phys. Rev. A **94**, 033801 (2016).

²² M. Schiró, C. Joshi, M. Bordyuh, R. Fazio, J. Keeling, and H. E. Türeci, Phys. Rev. Lett. **116**, 143603 (2016).

²³ E. T. Owen, J. Jin, D. Rossini, R. Fazio, and M. J. Hartmann, ArXiv e-prints (2017), arXiv:1711.11309 [quant-ph].

²⁴ A. Biella, F. Storme, J. Lebreuilly, D. Rossini, R. Fazio, I. Carusotto, and C. Ciuti, Phys. Rev. A **96**, 023839 (2017).

²⁵ M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch, and A. A. Houck, Phys. Rev. X **7**, 011016 (2017).

²⁶ R. Ma, C. Owens, A. Houck, D. I. Schuster, and J. Simon, Phys. Rev. A **95**, 043811 (2017).

²⁷ R. Ma, B. Saxberg, C. Owens, N. Leung, Y. Lu, J. Simon, and D. I. Schuster, ArXiv e-prints (2018), arXiv:1807.11342 [cond-mat.quant-gas].

²⁸ J. Raftery, D. Sadri, S. Schmidt, H. E. Türeci, and A. A. Houck, Phys. Rev. X **4**, 031043 (2014).

²⁹ S. Hacohen-Gourgy, V. V. Ramasesh, C. De Grandi, I. Siddiqi, and S. M. Girvin, Phys. Rev. Lett. **115**, 240501 (2015).

³⁰ A. J. Hoffman, S. J. Srinivasan, S. Schmidt, L. Spietz, J. Aumentado, H. E. Türeci, and A. A. Houck, Phys. Rev. Lett. **107**, 053602 (2011).

³¹ J. Lebreuilly, A. Biella, F. Storme, D. Rossini, R. Fazio, C. Ciuti, and I. Carusotto, Phys. Rev. A **96**, 033828 (2017).

³² J. Lebreuilly, M. Wouters, and I. Carusotto, Comptes Rendus Physique **17**, 836 (2016), polariton physics / Physique des polaritons.

³³ M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B **40**, 546 (1989).

³⁴ O. Scarlatella, A. A. Clerk, and M. Schirò, arXiv e-prints, arXiv:1811.03518 (2018), arXiv:1811.03518 [quant-ph].

³⁵ S. Schmidt and G. Blatter, Phys. Rev. Lett. **103**, 086403 (2009).

³⁶ J. Koch and K. Le Hur, Phys. Rev. A **80**, 023811 (2009).

³⁷ K. Dunnett and M. H. Szymańska, Phys. Rev. B **93**, 195306 (2016).

³⁸ L. M. Sieberer, M. Buchhold, and S. Diehl, Reports on Progress in Physics **79**, 096001 (2016).

³⁹ S. Diehl, A. Tomadin, A. Micheli, R. Fazio, and P. Zoller, Phys. Rev. Lett. **105**, 015702 (2010).

⁴⁰ A. Tomadin, S. Diehl, and P. Zoller, Phys. Rev. A **83**, 013611 (2011).

⁴¹ C. De Dominicis, E. Brézin, and J. Zinn-Justin, Phys. Rev. B **12**, 4945 (1975).

⁴² P. Hohenberg and B. Halperin, Rev. Mod. Phys. **49**, 435 (1977).

⁴³ A. Rançon and N. Dupuis, Phys. Rev. B **84**, 174513 (2011).

⁴⁴ J. Marino and S. Diehl, Phys. Rev. B **94**, 085150 (2016).
45 B. Sciolla and G. Biroli, Phys. Rev. B **88**, 201110 (2013).

46 A. Maraga, A. Chiocchetta, A. Mitra, and A. Gambassi, Phys. Rev. E **92**, 042151 (2015).

47 A. Chiocchetta, M. Tavora, A. Gambassi, and A. Mitra, Phys. Rev. B **94**, 134311 (2016).

48 M. Ludwig and F. Marquardt, Phys. Rev. Lett. **111**, 073603 (2013).

49 B. A. Levitan, A. Metelmann, and A. A. Clerk, New Journal of Physics **18**, 093014 (2016).

50 V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. **116**, 250401 (2016).

51 D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. **117**, 090402 (2016).

52 C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi, Phys. Rev. B **94**, 085112 (2016).

53 J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I. D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe, Nature **543**, 217 (2017).

54 S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, Nature **543**, 221 (2017).

55 D. Witthaut, S. Wimberger, R. Burioni, and M. Timme, Nature Communications **8**, 14829 EP (2017).

56 N. Lörch, S. E. Nigg, A. Nunnenkamp, R. P. Tiwari, and C. Bruder, Phys. Rev. Lett. **118**, 243602 (2017).

57 T. Nag, R.-J. Slager, T. Higuchi, and T. Oka, arXiv e-prints, arXiv:1802.02161 (2018), arXiv:1802.02161 [cond-mat.str-el].
45 B. Sciolla and G. Biroli, Phys. Rev. B **88**, 201110 (2013).

46 A. Maraga, A. Chiocchetta, A. Mitra, and A. Gambassi, Phys. Rev. E **92**, 042151 (2015).

47 A. Chiocchetta, M. Tavora, A. Gambassi, and A. Mitra, Phys. Rev. B **94**, 134311 (2016).

48 M. Ludwig and F. Marquardt, Phys. Rev. Lett. **111**, 073603 (2013).

49 B. A. Levitan, A. Metelmann, and A. A. Clerk, New Journal of Physics **18**, 093014 (2016).

50 V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. **116**, 250401 (2016).

51 D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. **117**, 090402 (2016).

52 C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi, Phys. Rev. B **94**, 085112 (2016).

53 J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I. D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe, Nature **543**, 217 (2017).

54 S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, Nature **543**, 221 (2017).

55 D. Witthaut, S. Wimberger, R. Burioni, and M. Timme, Nature Communications **8**, 14829 EP (2017).

56 N. Lörch, S. E. Nigg, A. Nunnenkamp, R. P. Tiwari, and C. Bruder, Phys. Rev. Lett. **118**, 243602 (2017).

57 T. Nag, R.-J. Slager, T. Higuchi, and T. Oka, arXiv e-prints, arXiv:1802.02161 (2018), arXiv:1802.02161 [cond-mat.str-el].