Dynamical Mean-Field Theory for Markovian Open Quantum Many-Body Systems

Open quantum many body systems describe a number of experimental platforms relevant for quantum simulations, ranging from arrays of superconducting circuits to ultracold atoms in optical lattices. Their theoretical understanding is hampered by their large Hilbert space and by their intrinsic nonequilibrium nature, limiting the applicability of many traditional approaches. In this work we extend the nonequilibrium bosonic Dynamical Mean Field Theory (DMFT) to Markovian open quantum systems. Within DMFT, a Lindblad master equation describing a lattice of dissipative bosonic particles is mapped onto an impurity problem describing a single site embedded in its Markovian environment and coupled to a self-consistent field and to a non-Markovian bath, where the latter accounts for finite lattice connectivity corrections beyond Gutzwiller mean-field theory. We develop a non-perturbative approach to solve this bosonic impurity problem, which treats the non-Markovian bath in a non-crossing approximation. As a first application, we address the steady-state of a driven-dissipative Bose-Hubbard model with two-body losses and incoherent pump. We show that DMFT captures hopping-induced dissipative processes, completely missed in Gutzwiller mean-field theory, which crucially determine the properties of the normal phase, including the redistribution of steady-state populations, the suppression of local gain and the emergence of a stationary quantum-Zeno regime. We argue that these processes compete with coherent hopping to determine the phase transition towards a non-equilibrium superfluid, leading to a strong renormalization of the phase boundary at finite-connectivity. We show that this transition occurs as a finite-frequency instability, leading to an oscillating-in-time order parameter, that we connect with a quantum many-body synchronization transition of an array of quantum van der Pol oscillators.

Open quantum many body systems describe a number of experimental platforms relevant for quantum simulations, ranging from arrays of superconducting circuits hosting correlated states of light, to ultracold atoms in optical lattices in presence of controlled dissipative processes. Their theoretical understanding is hampered by the exponential scaling of their Hilbert space and by their intrinsic nonequilibrium nature, limiting the applicability of many traditional approaches. In this work we extend the nonequilibrium bosonic Dynamical Mean Field Theory (DMFT) to Markovian open quantum systems. Within DMFT, a Lindblad master equation describing a lattice of dissipative bosonic particles is mapped onto an impurity problem describing a single site embedded in its Markovian environment and coupled to a self-consistent field and to a non-Markovian bath, where the latter accounts for fluctuations beyond Gutzwiller mean-field theory due to the finite lattice connectivity.We develop non-perturbative approach to solve this bosonic impurity problem, which dresses the impurity, featuring Markovian dissipative channels, with the non-Markovian bath, in a self-consistent scheme based on a resummation of non-crossing diagrams. As a first application of our approach, we address the steady-state of a driven-dissipative Bose-Hubbard model with two-body losses and incoherent pump. We show that DMFT captures hopping-induced dissipative processes, completely missed in Gutzwiller mean-field theory, which crucially determine the properties of the normal phase, including the redistribution of steady-state populations, the suppression of local gain and the emergence of a stationary quantum-Zeno regime. We argue that these processes compete with coherent hopping to determine the phase transition towards a non-equilibrium superfluid, leading to a strong renormalization of the phase boundary at finite-connectivity. We show that this transition occurs as a finite-frequency instability, leading to an oscillating in time order parameter, that we connect with a quantum many-body synchronization transition of an array of quantum van der Pol oscillators.

I. INTRODUCTION
Developments in quantum science and quantum engineering have brought forth a variety of platforms to store, control and process information at genuine quantum levels. Examples include trapped atoms and ions [1], quantum cavity-QED systems [2], superconducting qubits [3] or quantum optomechanical systems [4]. These architectures are not only of great relevance for quantum technologies but also for the quantum simulation of emergent collective many-body phenomena. We now have several experimental quantum simulators worldwide running on a variety of architectures, from ultracold atoms in optical lattices [5], Rydberg gases [6], trapped ions [7] and coupled cavity arrays [8]. Such simulators have led to significant advances in our understanding of quantum manybody phases and offer us an opportunity to address deep unanswered questions concerning the behavior of quantum matter in novel unexplored regimes, particularly far away from thermal equilibrium.
Many of these systems are typically characterized by excitations with a finite lifetime, due to unavoidable losses, dephasing and decoherence processes originating from their coupling to external environments and therefore feature an intrinsic nonequilibrium nature. Arrays of circuit QED cavities, for example, are emerging as a unique platform for the quantum simulation of drivendissipative quantum many body systems [8][9][10][11], where the balance of drive and loss processes and the presence of strong non-linearities, allows one to reach interesting non-equilibrium stationary states. Experiments have recently started to show a variety of dissipative quantum phases and phase transitions [12][13][14][15], including the recent experimental realization of a dissipatively stabilised Mott insulator [16]. On a different front, recent advances with ultracold atoms in optical lattices allow the engineering of controlled dissipative processes, such as correlated losses [17,18] or heating due to spontaneous emission [19,20] and to study the resulting quantum many body dynamics over long time scales. These experimental settings offer fresh new inputs for quantum simulation of strongly correlated driven-dissipative quantum manybody systems, at the interface between quantum optics and condensed matter physics [21][22][23][24][25][26].
From a theoretical standpoint these experimental platforms can be well described as open Markovian quantum systems, of either fermionic/bosonic particles or quantum spins, evolving through a Lindblad master equation for the many body density matrix ρ [27] The crucial aspect of this problem lies in the interplay between unitary dynamics and the dissipative evolution controlled by jump operators L α , out of which non trivial stationary states can be engineered [28][29][30][31]. Open Markovian quantum many body systems are predicted to display a broad range of new transient dynamical phenomena [32][33][34][35] and dissipative quantum phase transitions [36][37][38][39][40][41].
Solving the Lindblad equation for many-body systems is extremely challenging. Exact numerical solutions based on the diagonalization of the Lindbladian superoperator, or direct time evolution, are limited to very small systems, and only slightly larger systems can be addressed with quantum trajectories [42]. For one dimensional systems an efficient representation of the problem in the language of matrix product operator is possible [43,44] and has been successfully used in the recent past [45], however its extension to higher dimensional systems poses problems, although some recent results have been obtained [46][47][48]. As a result, a number of theoretical approaches have been developed in recent years to tackle driven-dissipative lattice systems [47,[49][50][51][52][53][54][55][56][57][58][59]. Driven-dissipative correlated bosons, such as those described by Bose-Hubbard and related models, are particularly challenging to tackle numerically due to the unbounded Hilbert space.
A powerful non-perturbative approach to quantum lattice models, based on the limit of large lattice connectivity z [60,61], is provided by the dynamical mean field theory (DMFT) of correlated electrons [62] and bosons [63][64][65]. When applied to equilibrium lattice bosons DMFT captures, non-perturbatively, the 1/z corrections to Gutzwiller mean-field theory through the solution of a quantum impurity model. Originally introduced for equilibrium problems, DMFT has been successfully applied to a variety of nonequilibrium fermionic problems with or without dissipation [66], including Markovian fermions [67], and to study unitary dynamics of correlated bosons [68].
In this work we introduce a new approach to Markovian bosonic quantum many-body systems. Specifically we extend the bosonic non-equilibrium DMFT [68] to systems evolving under the Lindblad master equation (1) and combine it with a strong coupling impurity solver tailored for correlated dissipative processes described by many-body jump operators. In the large connectivity limit DMFT maps the Lindblad dynamics (1) onto a dissipative quantum impurity model describing a single interacting bosonic mode, subject to many-body Markovian dissipation, coupled to a coherent field and a non-Markovian environment both to be determined selfconsistently. The non-Markovian DMFT bath describes the effect of hopping processes from neighboring sites, which are completely missed by Gutzwiller mean-field theory. We will show that these processes are particularly interesting in open quantum systems since they not only introduce coherent effects but also provide new dissipative channels.
Solving the resulting quantum impurity model in presence of multiple many-body jump operators remains highly non-trivial and state of the art impurity solvers for non-equilibrium DMFT cannot be efficiently used in this setting. Here we build upon our recent work on Markovian impurity models [69] to develop and benchmark a DMFT solver for driven-dissipative bosonic systems. This approach is based on a super-operator hybridization expansion and resummation of an infinite class of diagrams known as non-crossing approximation (NCA).
As first application of the DMFT/NCA method we study a Bose-Hubbard model subject to two-particle losses and single particle incoherent pumping. This model is directly relevant for recent experiments with ultracold bosonic atoms in optical lattices under controlled dissipation [18,20] as well as for arrays of superconducting circuits [16]. We predict the emergence of a dissipative phase transition from a normal to a superfluid phase, where above a critical hopping or pump strength the system spontaneously develops a coherent field oscillating in time, and discuss the effect of finite-connectivity corrections to Gutzwiller mean-field. We highlight how this transition can be naturally interpreted as the onset of many body synchronization in an array of quantum Van der Pol oscillators, a phenomenon which recently attracted lots of attention [70][71][72][73][74][75][76][77][78][79][80][81][82]. We show that DMFT allows to predict results which are qualitatively missed by Gutzwiller mean-field theory, including the redistribution of steady-state population and the suppression of gain due to hopping processes, a stationary quantum Zeno regime and a new competition between coherent and dissipative processes towards symmetry breaking. These results reflect in large quantitative corrections to the Gutzwiller mean-field phase diagram.
The paper is organized as follows. In the next section II we summarize the main concepts and results of this work, which will be presented in more detail in the rest of the paper. In Sec. III we introduce DMFT in more detail and discuss its physical content. In Sec. IV we discuss two methods to solve the quantum impurity problem: a strong coupling perturbative approach and a more powerful self-consistent NCA method. In Sec. V we present our results for an interacting Bose-Hubbard driven-dissipative lattice problem. Sec. VI is devoted to conclusions. In the Appendixes we provide technical details on the derivation of DMFT (Appendix A), a non trivial consistency check on the NCA approximation (Appendix B) and various analytical results quoted in the main text (Appendixes C-G).   in the main text) is mapped in the large connectivity limit z 1 onto a single site problem, with interactions, local Markovian drive and losses, coupled to a self-consistent environment (right, top panel). This includes a coherent drive and a non-Markovian bath (characterized by two independent hybridization functions ∆ R (t, t ) and ∆ K (t, t )) which accounts for fluctuations due to neighboring sites. In the infinite coordination limit, z = ∞, one recovers the Gutzwiller mean-field theory (right, bottom panel), where the effect of the lattice is simply encoded in a self-consistent field.

II. SUMMARY OF MAIN RESULTS
In this section we present a summary of the main results of this work, which will be discussed more in detail in the rest of the paper. In particular in Sec. II A we discuss the formulation of DMFT for Markovian Bosons and of the NCA impurity solver, in Sec. II B we discuss our theoretical approach in relation to prior work on open quantum many-body systems, while in Sec. II C we present the application to a driven-dissipative Bose-Hubbard model with two-body losses and single-particle incoherent pump.

A. Dynamical Mean-Field Theory for Open Markovian Bosons
The class of models that we consider describe drivendissipative bosonic particles on a lattice with coordination number z. The many-body density matrix of the system evolves according to a Lindblad master equation with a set of local jump operators for each lattice site L iµ , accounting for dissipative processes, and with coherent evolution governed by the Hamiltonian Here [b i , b † i ] = 1 are bosonic modes localized around the lattice site i, coupled together by a nearest-neighbours hopping term J. H loc [b † i , b i ] is a generic local Hamiltonian, which can contain arbitrary local interactions. In order for the problem to remain well defined in the limit of infinite connectivity z → ∞, to which we will compare our DMFT results, we scale the hopping with a 1/z factor which is the correct scaling for bosons [63,64,83](see also Sec. III A for further details).
The DMFT approach considers the master equation (2) in the limit of large, but finite, lattice connectivity z 1. In fact, when the number of neighbors z of a given lattice site is large, statistical and quantum fluctuations induced by the neighboring sites get small and can be treated in an approximate way, while the local, on-site physics is accounted for exactly. As we discuss in detail in Sec. III, in the z 1 limit the Lindblad lattice problem formally maps onto a quantum impurity model describing an interacting Markovian singlesite, characterized by the same local Hamiltonian H loc and local jump operators L iµ entering Eq. (1), coupled to a time-dependent field Φ eff (t) acting as a coherent drive and a non-Markovian quantum bath (Fig. 1, top panel). These take into account the effect of the neighboring sites and have to be determined self-consistently through the calculation of impurity properties. As a result of the non-equilibrium nature of the problem, the non-Markovian environment is described in terms of two independent bath hybridization functions, the retarded ∆ R (t, t ) and Keldysh ∆ K (t, t ) which in a stationarystate encode information about spectrum and occupation of the single-particle excitations. In a generic nonequilibrium condition these are independent and not related by the fluctuation-dissipation theorem.
To appreciate the physical content of DMFT it is instructive to compare it with the widely used Gutzwiller mean-field theory. As we will show in Sec. III A the latter coincides with the z → ∞ solution of the manybody master equation and corresponds to DMFT when the non-Markovian bath is set to zero. Gutzwiller meanfield theory amounts to decouple the hopping term in the Hamiltonian (3), or equivalently assumes a product-state density matrix over different lattice sites, thus reducing the full many-body problem to a single site coupled to a self-consistent coherent field (Fig. 1, bottom panel). An obvious shortcoming of the Gutzwiller approach is that it cannot capture any coherent or dissipative processes arising from the coupling to neighboring sites, unless the system is in a broken-symmetry phase with a non-vanishing local order parameter, leading to a finite self-consistent field. The result is a particularly poor description of strongly interacting, yet incoherent, normal phases such as bosonic Mott insulators or many-body quantum Zeno phases that we discuss in this work, whose local properties within Gutzwller are completely independent on the hopping and identical to those of an isolated site. In this perspective DMFT accounts non-perturbatively, through the solution of a quantum impurity model with a non-Markovian bath ∆ ∼ 1/z, for finite-connectivity corrections to Gutzwiller mean field theory, thus capturing fluctuations induced by the neighboring sites even in absence of an order parameter.
Although simplified with respect to the full master equation, the solution of DMFT equations and in particular of the quantum impurity problem sketched in Figure 1 still poses tremendous challenges. In particular the presence a Markovian environment containing arbitrary, possibly non-linear, jump operators, in addition to local interactions and the non-Markovian DMFT bath makes this problem hard to be solved efficiently with state of the art approaches for non-equilibrium DMFT. A major result of the present work is the development and benchmark of an impurity solver for Markovian bosonic DMFT based on the super-operator hybridization expansion formulated in Ref. 69 and applied there to a simple noninteracting fermionic resonant level model. As we will discuss more in detail in Sec. IV B this approach fully captures the underlying Markovian dynamics of the impurity problem in Fig 1 and accounts for the non-Markovian bath ∆ through the resummation of an infinite class of diagrams in the impurity-bath coupling known as non-crossing approximation (NCA).
As we will discuss further on in the paper, the self-consistent nature of the non-crossing approximation (NCA) we use, as opposed to bare perturbative expansions to which we will compare our results, is crucial to fully capture the non-trivial correlations associated to the DMFT bath. We give a more complete picture of the DMFT formalism, including a discussion of the basic equations and of impurity solvers in Sec. III and Sec. IV.

B. Relation to Prior Works
Here we wish to relate our approach with respect to previous works on nonequilibrium dissipative many-body systems. In the context of fermionic nonequilibrium DMFT, dissipation at single particle level (i.e. tunneling to external metallic leads) has been included before in several works, focusing for example on steadystate transport [84][85][86][87][88][89], Floquet driven systems [90][91][92] or photodoping [93]. We note that this type of dissipation is straightforward to handle within DMFT and pose no additional methodological challenges since it can be included within any impurity solver used for nonequilibrium DMFT in absence of dissipation. On the other hand many-body dissipative processes, such as those we focus here in the Lindblad context or those modeling the coupling between fermions and bosonic baths, are more challenging to handle since they induce effective interactions. Up to date these have been included in nonequilibrium DMFT studies of dissipative problems only through perturbative expansions [94][95][96][97][98]. In this respect our work goes beyond those studies in that all Markovian dissipative couplings (single and many-body) are treated on the same footing and encoded in the local Lindbladian of the impurity model, which opens up the possibilities for non-perturbative treatments of those couplings. This strategy is similar in spirit to what done for Markovian fermionic systems in Ref. 67, where a discretization of the DMFT bath was used to solve the impurity problem with exact diagonalization. Here instead we use the NCA impurity solver which works directly in the thermodynamic limit of an infinite bath and does not introduce any discretization, which would be particularly severe for bosonic problems such as the one we focus here.
In the context of Markovian quantum many-body systems there have been recent methodological developments to include correlations beyond mean-field theory. Although a precise comparison with DMFT is beyond the scope of this work, it is worth to discuss some of those methods here. The Cluster Mean-Field Theory [51], developed for driven-dissipative quantum spin models, is similar to DMFT in that it adds short-range correlations on top of a Gutzwiller mean-field, although this is achieved through the exact solution of a finite-size cluster, rather than through an infinite (non-interacting) bath. The Corner-Space Renormalization Method [49] performs a diagonalization of the Lindbladian in a cor-ner of the full Hilbert space, whose size is iteratively increased. As opposed to DMFT which works in the thermodynamic limit, this is a finite-size method, which can however produce converged results for sizes much larger than brute force methods [99]. Both those approaches focus naturally on static correlations encoded in the stationary state density matrix while DMFT is constructed around the frequency-resolved Green's functions.
C. Application to a Driven-Dissipative Bose Hubbard Lattice In this work, we apply DMFT to a lattice model of driven-dissipative interacting bosons by specifying the local Hamiltonian and local jump operators entering Equations (2)(3). We consider for the former i.e. a characteristic frequency ω 0 and on-site Kerr non linearity of strength U while for the latter we consider two kinds (µ = 1, 2) of jump operators for each lattice site i, We emphasize the correlated nature of the dissipative process encoded by L i2 which acts only on states with multiple bosonic occupancy. This term will play a key role for our results. Interestingly this kind of loss process can be realized both with ultracold atoms [17,18] as well with superconducting circuits [100]. The resulting lattice model, Eq. (3-6), therefore describes a drivendissipative realisation of the Bose-Hubbard model [83] whose properties in presence of incoherent drive and dissipation has been the subject of tremendous attention recently [32,40,54,[101][102][103][104][105][106][107][108][109]. The specific form of dissipation we consider in Eq. (5) is rather unexplored in a many-body context, although few results are available in the literature that we will recall briefly here. The many body master equation (1,4-6) has a global U (1) symmetry, associated with the invariance of the Liouvillian with respect to the transformation b i → b i e iθ , and is time translational invariant (TTI). In the limit of a large number of bosons per site one expects a semiclassical description to work. The model reduces then to a discretized version of Gross-Pitaevskii equation, largely studied in the context of exciton-polariton condensation [36], which predicts a coherent phase of bosons for any r > 0, independently of J/U . This phase, which spontaneously breaks both the U (1) and TTI symmetry corresponds to a nonequilibrium superfluid. In the opposite regime of uncoupled sites, J/U = 0, the steady-state density matrix is known analytically from Refs. [110,111] and describes an incoherent state: it is a statistical mixture of Fock states with b i = 0, as might be expected given the lack of any coherent driving. How those two different phases are connected upon increasing J/U , in the quantum regime of few bosons per site and finite lattice connectivity, is one of the main focus of this paper. In Figure 2 we plot the DMFT phase diagram for our Bose-Hubbard model as a function of r (the dimensionless pump-to-loss ratio) and J/U , for different values of the lattice connectivity z, together with the Gutzwiller mean field phase boundary corresponding to the z = ∞ limit [112]. For a given fixed value of z we find a critical line r c (J) separating a small-hopping normal phase where the bosons remain incoherent, b i = 0, from a large hopping phase where the system develops a local order parameter breaking the U (1) symmetry of the master equation. A first striking result that clearly appears from Figure 2 is that upon decreasing the connectivity, i.e. increasing the strength of fluctuations on top of the Gutzwiller mean field result, the phase diagram changes substantially. In particular the broken symmetry phase shrinks and moves toward larger values of pump and hopping. Interestingly, the DMFT corrections to the phase diagram turn out to be substantially larger than for equilibrium lattice bosons [64,65,68]. The effect of finiteconnectivity fluctuations is however not only quantitative. As we are going to discuss below, and more extensively in Sec. V, these corrections arise from a qualitatively new physics that is captured by the DMFT/NCA description of the normal phase and completely missed by Gutzwiller mean-field theory.
As we will discuss in more details in Sec. V A-V B, the normal phase of our model come with unique nonequilibrium properties, inherited from the local many-body physics of the single site problem [113]. Above a pump threshold r ndos the system develops a negative density of states (NDOS) at positive frequencies, a signature of incipient gain, i.e. energy emission in response to a weak coherent drive. Upon further increasing the pump to loss ratio above r inv > r ndos the steady-state reduced density matrix shows population inversion, namely higher energy states become more occupied than lower energy ones. Within Gutzwiller mean-field theory, which describes the normal phase as a product state of single sites, those scales are independent from the hopping J. DMFT instead shows that fluctuations due to finite connectivity reshape completely the spectral and distribution properties of the normal phase, leading in particular to a suppression of NDoS and population inversion. This arises from hopping-induced losses, a hallmark of the interplay between coherent and dissipative dynamics in our model, which are the key physics captured by DMFT through the non-Markovian bath.
In Sec. V C we show that an interplay of NDoS and sufficiently strong hopping J controls the true normal phase instability for values of the pump above r c (J), when the system develops full phase coherence and enters the superfluid phase. In particular we find that the unstable mode is modulated in time and that the system displays a finite frequency phase transition corresponding to an order parameter which develops undamped oscillations, thus breaking TTI [109,114,115]. The large reduction of the normal phase at finite connectivity can be interpreted as an effect due to hopping-induced losses arising from the non-Markovian DMFT bath, which is able to wipe out the NDoS and absorb part of the energy emitted by the system, as we discuss more in detail in Sec. V C 2.This mechanism for the destruction of an ordered phase is of genuine nonequilibrium origin and cannot be understood in terms of an effective heating. Indeed, as we show in Sec. V C 3 while an effective thermalisation is captured by DMFT through an effective temperature, this remains comparable to the Gutzwiller mean-field result up to small values of the connectivity and therefore cannot explain by itself the substantial reshape of the phase diagram.
The finite-frequency transition towards an oscillating nonequilibrium superfluid shares similarities with phenomena such as lasing, limit cycles and synchronization. As we show more in detail in Sec. V D, the drivendissipative Bose-Hubbard (3)(4)(5) reduces in the semiclassical limit to an array of coupled classical van der Pol (vdP) oscillators [116][117][118], which admits a stable limit cycle phase, a coherent phase with an order parameter oscillating in time at finite frequency, for any finite pump r > 0 and any coupling J. In the quantum regime of few bosons per site, the picture qualitatively changes and a FIG. 3. Emergence of Quantum Zeno Regime at fixed drive rate f /U = rη/U = 0.013 and increasing two-particle loss rate η. Parameters: tmax = 20, dt = 0.002, dimH = 6. Probability of having one boson per site in the steady-state ρ1 as a function of η/U , rescaled by its value at η/U = 2.67 for convenience. It follows a non-monotonic behaviour with a minimum at η = U , which is more pronounced the larger the hopping J. This behaviour is a manifestation of quantum Zeno physics, controlled by the hopping-induced loss rate Γ eff as we show clearly in the inset, which is instead completely missed by Gutzwiller.
Finally, in order to highlight the role of hoppinginduced dissipative processes and the qualitative differences between DMFT/NCA and Gutzwiller, in Sec. V E we consider the limit of large two body losses η J for our Bose-Hubbard model. We note that this regime is experimentally accessible with ultracold gases in presence of inelastic collisions [17,18]. For large two-body losses, and in absence of any external pump, perturbation theory shows that one can restrict the Bose-Hubbard dynamics to a subspace made by hard-core bosonic states, the dark states of the local dissipator [120]. Within this quantum Zeno subspace [121,122] the dominant dissipative processes left are those among neighboring sites, controlled by a hopping-induced loss rate Γ eff = 2(J/z) 2 U 2 +η 2 η. This scale was shown to affect the transient dynamics [17,120,123], featuring a power-law decay towards the trivial zerodensity steady-state. Here we use DMFT/NCA to show that in presence of an additional small residual pump f η one can stabilize a Quantum Zeno stationary state in which physical properties are controlled by the scale Γ eff (see Sec. V E for a detailed discussion). An example is provided in Figure 3, where we plot the DMFT/NCA steady-state occupation probability of the n = 1 bosonic state versus η/U , for different values of J/U and compare it with Gutzwiller results. The latter shows, independently of J, a very weak dependence from η which would be completely absent if not for the small residual pump. DMFT results instead show a clear non-monotonous behavior as J/U increases, with a minimum at η ∼ U . This is a signature of the emergence of the Zeno scale Γ eff controlling the physics, as we clearly reveal in the inset, where a scaling collapse is shown. We emphasize that the emergence of a Quantum Zeno stationary state represents a stringent test for the ability of DMFT of capturing hopping-induced losses, which are instead completely missed by Gutzwiller mean-field theory.
We give a more complete picture of our results for the Bose-Hubbard model in Sec. V.

III. DYNAMICAL MEAN-FIELD THEORY FOR MARKOVIAN BOSONS
In this section we present the formalism of DMFT for Markovian bosons, including the basic self-consistency equations, its formal relation with Gutzwiller mean field theory and its physical interpretation, leaving its derivation to Appendix A. The starting point is to cast the Lindblad master equation (2) in the language of nonequilibrium Keldysh field theory, as discussed extensively in the literature [52]. The result is an action written in terms of bosonic coherent fieldsb iα , b iα on each lattice site i and on the upper/lower Keldysh contours, α = ±, which takes the form where the Lindbladian L is given by with H α , L iµα are the expectation values of the Hamiltonian (3) and of the jump operators (5), expressed in terms of creation and annihilation operators belonging to the α contour, taken on bosonic coherent states. The full solution of the Keldysh action in Eq. (7) is of course not possible in general, due to the coupling between many interacting modes and the presence of interaction, drive and dissipation.
The key idea of DMFT is to write down an effective Keldysh action for the boson field of a single site of the lattice, obtained after integrating out all its neighbors [68]. As we show in Appendix A, in the limit of large lattice connectivity, z 1, this effective action has the closed form where α, β are Keldysh contour indices and we have dropped the site index from the local boson field for simplicity (we assume translational invariance) and grouped together creation/annihilation fields into a Nambu field The above local Keldysh action describes a drivendissipative quantum impurity model [69]. The first term in Eq.
is the local, on-site, contribution of the original lattice action (7)(8) and therefore includes interactions, as well as Markovian incoherent drive and dissipation leading to off-diagonal terms in Keldysh space. The second and third terms describe the feedback of the rest of lattice onto the site through its neighbors, in terms of an effective coherent drive Φ † eff α (t) and an effective non-Markovian bath with hybridization function ∆ αβ (t, t ). Both these quantities have to be determined self-consistently, in particular the effective coherent drive reads and has two contributions, the first coming from the average of the bosonic field as in Gutzwiller mean field theory and the second one coming from the non-Markovian bath, a non-trivial finite z correction accounting for the feedback of neighboring sites on the local effective field [64,65,68]. This latter term, whose origin will be discussed more in detail in Appendix A, plays a key role within DMFT, in particular for what concerns the corrections to phase diagram as we discuss in Sec. V C-V C 2. The self-consistency relation for the Green's function depends on the specific choice of the lattice. In the following we will use the simplified relation for the Bethe lattice [68] which directly relates the hybridization function of the non-Markovian bath to the impurity connected Green's function The DMFT solution of the original Markovian lattice problem thus requires one to solve the Keldysh action (9), computing in particular the impurity Green's function (14) and the average of the bosonic field (12), for given values of the non-Markovian bath ∆ and effective field Φ, and to iterate (11-13) until self-consistency.

A. Limit of infinite coordination number: Gutzwiller Mean Field Theory
It is instructive at this point to take explicitly the limit of infinite coordination number z → ∞. In this limit, the DMFT effective action (15) becomes completely local in time since the non-Markovian bath scales as 1/z (See Eq. 13), and as such can be unfolded back into a master equation for a single-site density matrix ρ, which satisfies where L is the local part of the Lindbladian and the feedback from the neighboring sites is carried by . This corresponds to a factorized Gutzwiller-like ansatz for the lattice many body density matrix where i is the site index and ρ i (t) ≡ ρ(t) because of translational invariance. In other words we have explicitly shown that, as for equilibrium or closed systems [64,68] also for driven-dissipative lattice systems the infinite connectivity limit of bosons coincides with Gutzwiller meanfield theory. We note that when Φ = 0, this meanfield describes completely uncoupled sites, while DMFT (z < ∞) captures the feedback from neighboring sites through the self-consistent bath ∆. In the following section we are going to add some physical intuition on how the DMFT action (9) describes the effect of neighboring sites through a fictitious non-Markovian bath.

B. DMFT Effective Action in the Classical/Quantum basis
We now give a physical interpretation to the various terms entering the DMFT effective action in Eq. (9), in particular to the non-Markovian term. It is useful to introduce the so called classical and quantum components of the bosonic field in terms of which we can re-write the Keldysh action as In this basis only two independent combinations of the non-Markovian kernels ∆ αβ enter, namely the retarded component , which couples the classical and quantum components of the field and encodes the spectral properties of the bath resulting in a frequency dependent damping for the bosonic mode and the Keldysh component encoding the occupation of the bath and resulting in a frequency dependent noise for the bosonic mode. It is worth stressing that the above Keldysh action contains quadratic, non-Markovian terms in the classical/quantum fields as well as non-linearities and higher powers of the classical/quantum fields included in the local part of the action S loc [b † cl/q , b cl/q ]. While the structure of Eq. (18) is a generic feature of DMFT, the local part of the action depends on the particular form of local interaction and jump operators.
Finally, we can express also the impurity Green's functions Eq (14) in this basis to obtain the retarded Green's function and the Keldysh one Those correlation functions contain crucial physical information about the local physics of the driven-dissipative lattice problem. The retarded Green's function in particular encodes information about the local excitation spectrum of the system and it is known to be a crucial probe for the transition from delocalization to Mottness in strongly correlated systems at equilibrium [62]. For open Markovian quantum systems the retarded Green's function contains, much like for closed equilibrium systems, information on the structure of the excitations on top of the stationary state [113] and it directly probes dissipative phase transitions where those excitations become unstable. Its poles correspond to eigenvalues of the Lindbladian in the single particle sector, which come with a characteristic frequency and lifetime, and their (possibly complex) weight. The retarded Green's function has also a directly physical meaning: it describes the linear response of the expectation b(t) to a weak, classical field h(t ), which couples linearly to b † . In the case where b describes a photonic cavity mode, G R (t) can be directly measured by weakly coupling the cavity to an input-output waveguide and measuring the reflection of a weak probe tone (see e.g. [124,125]).
The Keldysh Green's function on the other hand contains information about how the finite frequency excita-tions above the stationary state are populated. In thermal equilibrium those two functions are not independent, but constrained to satisfy the fluctuation-dissipation theorem [126].

C. Computing Lattice Quantities
Solving the DMFT effective action and computing the impurity Green's functions (19) gives direct information on the local properties of the driven-dissipative lattice problem. Furthermore one can access non-local quantities, such as momentum distribution or non local correlation functions, through the lattice Green's functions at momentum k These satisfies a Dyson equation with a lattice self-energy Σ αβ (t, t ), that within DMFT is momentum independent [61,66], and coincides with the self-energy of the impurity problem where in the above equations ⊗ indicates time convolutions, g αβ (t, t ) are the Green's functions of the quantum impurity problem with no interactions, but including the non-Markovian bath ∆ and g αβ k (t, t ) are the noninteracting lattice Green's functions.

IV. QUANTUM IMPURITY SOLVERS
The main challenge behind our DMFT approach is to solve the Markovian quantum impurity model described by the Keldysh action (9), computing in particular the Green's functions. We stress that this remains a difficult task due to the presence of interactions on the impurity site, non-linear jump operators (such as our two-body losses) and the non-Markovian DMFT bath. While several impurity solvers have been developed in recent years for non-equilibrium DMFT [66], none of them can be efficiently applied in our case (See Sec. II B for a detailed discussion). To make progress we take explicit advantage of the Markovian structure of the impurity, which allows to treat non-linear jump operators as dissipative couplings of a local Lindbladian. This unleashes the possibility of developing strong-coupling impurity solvers for bosonic Markovian problems, which treat exactly the local Lindblad problem and include the effect of the non-Markovian DMFT bath through perturbative or non-perturbative schemes. We note that for non-equilibrium closed systems these strong-coupling methods represent the current state of the art of DMFT impurity solvers [66]. Here we develop two such schemes for bosonic Markovian systems, the Hubbard-I approximation and the more powerful Non-Crossing-Approximation, that we both present below. We comment in the Sec. VI on possible methodological extensions.

A. Hubbard-I Approximation
The simplest approximation to solve the impurity problem (9) is based on perturbation theory in the non-Markovian bath kernel ∆, and its lowest order is known as the Hubbard-I approximation [62,127]. As we will see this approach already gives a hopping dependence of correlation functions which goes beyond Gutzwiller meanfield theory, but misses important correlations due to the non-Markovian bath. Our DMFT approach will be based on the more powerful non-crossing approximation solver which we will introduce in the next section, but we will use Hubbard-I results for comparison and to motivate the need of a more powerful solver.
For simplicity, we formulate Hubbard-I in the normal phase, where Φ = 0 and anomalous correlation functions vanish, thus we can restrict to the first Nambu component and refer to it with non-bold symbols, e.g. G αβ = G αβ 11 where α, β = ± are Keldysh indexes. We also focus on the stationary state regime, where Green's functions depend on time differences and we can move to the frequency domain i.e. G αβ (ω), which is the case we will consider in our application in Sec. V.
The impurity Green's function obeys a Dyson equation, see Eq. (22), in terms of a self-energy Σ αβ (ω) which contains the effect of interaction, incoherent drive and dissipation and which is in general a functional of the non-Markovian bath kernel ∆ αβ (ω). Hubbard-I consists in approximating the impurity self-energy by its value for ∆ = 0, i.e. in absence of the bath, when it can be written as Here G αβ 0 (ω) is the Green's function of the impurity site with interaction, incoherent drive and dissipation, but without the bath (the latter condition is indicated by the index 0); it can be computed numerically [113]. In contrast g αβ 0 (ω) corresponds to the Green's function of the impurity site in absence of the bath and without interactions (lowercase letter), which is known analytically. Plugging this self-energy back in the Dyson equation (22), and using the self-consistency condition on the Bethe Lattice, we obtain a closed matrix equation for the Keldysh components of the local lattice Green's functions The expressions of the retarded and Keldysh components are given explicitly in appendix E. In the appendix, we also show that the Hubbard-I approximation, despite introducing a beyond mean-field hopping dependence of Green's functions, still yields the same phase diagram as mean-field, motivating the need for a more powerful solver.

B. Super-Operator Hybridization Expansion and Non-Crossing Approximation
To go beyond the Hubbard-I approximation, we build upon the method recently formulated in Ref. [69] and applied so far only to a simple toy model fermionic system, to develop a DMFT/NCA impurity solver for bosonic Markovian systems. The idea is to perform a diagrammatic expansion in powers of the non-Markovian bath ∆ and to resum an infinite set of diagrams by solving a self-consistent Dyson-type equation. We remark that this expansion is carried out around an interacting problem, the single-site Markovian impurity, hence it is not based on Wick's theorem as in weak-coupling perturbation theories. As such, working directly with Green's functions is not convenient and the more natural formulation is in terms of evolution super-operators, that we will denote in the following with a hat [69]. We start by defining the evolution super-operatorV of the reduced density matrix of the impurity formally obtained by tracing out the bath degrees of freedom. We note that Eq. (25) assumes that at time t = 0 the non-Markovian bath is not entangled with the impurity site, i.e. in the original lattice problem the initial condition corresponds to the limit of decoupled sites.
Since the bath degrees of freedom are treated as noninteracting, only the single-particle Green's function of the bath enters the reduced dynamics, the hybridization function ∆ introduced in Eq. (9). Expanding the super-operatorV(t, 0) in powers of ∆ we obtain a series which can be represented diagrammatically as shown in Fig.  4, where bold solid lines describe the propagatorV(t, 0), dashed lines represent the hybridization function ∆ while simple solid lines represent the bare Markovian evolution This diagrammatic representation allows to define the self-energyŜ of the series as the sum of one-particleirreducible (1PI) diagrams, which cannot be cut into two disconnected parts by removing a solid line, and thus to formally resum the series into the Dyson equation We remark thatV 0 ,V andŜ here are super-operators and that the self-energyŜ is a functional of the propa-gatorV whose closed form is not known in general. The resulting series (26) generalizes to the case of Markovian impurities the hybridization expansion obtained for unitary quantum impurity models [128][129][130][131][132]. For the latter, exact resummation techniques based on diagrammatic Monte Carlo [133] have been employed but generically suffer from the so called sign problem, especially out of equilibrium, limiting the propagation time. Here instead we adopt a self-consistent approximation for the self-energyŜ. This can be written in general as a systematic expansion in diagrams with an increasing number of crossing hybridization lines, an approach which has been extensively used for unitary quantum impurity models [134][135][136][137]. The lowest order contribution is given by non-crossing diagrams, e.g. in Fig. 4, giving an explicit expression for the NCA self-energŷ In the above expression α, β = ± are Keldysh indices and a, b = {1, 2} are Nambu indices. Thus ∆ αβ ab is a given component of the bath hybridization function introduced in Eq. (9) , i.e. ∆ αβ ab = ∆ αβ ab . We also introduce the super-operators analogues of the Nambu fields of Eq. (10), that we define aŝ and denote their a Nambu component asb αa in Eq. (27). The Keldysh index α = ± for a super-operator specifies whether it should act from the left or the right of its argument, i.e.b and similarly forb † α . We notice that the self-energy depends on the full propagatorV, rather than on the bare oneV 0 , thus containing diagrams to all orders in ∆. Corrections to the NCA can be obtained systematically including self-energy diagrams with higher number of crossings, although the resulting computational cost increases. In this work, where we focus on the normal phase and its instability, we will limit ourselves to the NCA scheme, while we expect that to access the superfluid phase or for lower values of the connectivity higher order corrections would become important. We note in fact that self-energy corrections including higher number of crossings diagrams come with higher powers of the DMFT bath, which for bosons is of order J 2 /z (see Eq. (13)) and therefore subleading at least for large to moderate values of the connectivity.
Once the self-energyŜ is known in closed form, the propagatorV can be obtained numerically by solving Eqs. (26) and (27). To use this NCA impurity solver in our DMFT approach, we need to compute the oneparticle Green's functions of the impurity, Eq. (14). This can be obtained by taking the functional derivative with respect to ∆ of the partition function Z = tr [ρ imp (∞)] in Eq. (25) and using the Dyson equation forV (see Appendix G). The final result reads where as before we have written explicitly both the Keldysh indices α, β and the Nambu ones a, b and where Φ αa (t) = tr b αa ρ imp (t) . We notice that this result, which resembles a quantum regression theorem [27] for the non-Markovian mapV(t, t ) is only valid within NCA, while including higher order diagrams into the self-energy would lead to further terms which can be interpreted as vertex corrections. Finally, we conclude by emphasizing that the solver introduced in this section is different from other NCA approaches to quantum impurity models with or without dissipation [68,96,98,136,138], which treat at the non-crossing level all couplings to the baths. Here, by formulating the hybridization expansion at the superoperator level, we are able to fully capture the underlying local Markovian dynamics, resorting to an NCA approximation only for the non-Markovian DMFT bath. This introduces several differences with respect to the NCA literature, including the way the Green's functions are evaluated (see Eq. (30)) and in the way the stationarystate theory is constructed, as we discuss further in the next section.

Stationary state DMFT/NCA
While the formalism introduced so far allows to compute the whole transient dynamics, in this section we show how to directly address the stationary state properties of the system within our DMFT/NCA approach. At stationarity we expect the local Green's functions (14) and, through the self-consistent condition (13), the bath hybridization function ∆ αβ to depend only on time differences. We can then solve the NCA Dyson equation (26) for the stationary state propagatorV(t − t ) which also depends only on time differences. This allows to significantly reduce the computational cost for time-propagating this equation from O(t 3 max ) to O(t 2 max ), where t max is the maximum integration time.
A complete steady-state DMFT/NCA procedure requires to compute, in addition to the stationary state propagator, also the steady-state density matrix of the impurity ρ s ≡V(∞, 0)ρ imp (0), which is needed to evaluate the impurity Green's functions (See Eq. (30)). While in principle this would require to perform the full transient dynamics from an arbitrary initial condition, here we show how to obtain ρ s directly from the stationary state propagatorV(t − t ). We note that for Markovian open quantum systems the stationary-state density matrix can be directly obtained as zero-eigenvalue of the Lindblad supeoperator generating the dynamics. This argument however does not directly apply to the present case, since the DMFT bath makes the mapV(t, t ) non-Markovian. A generalized stationarity condition for the non-Markovian map (26) can be obtained [69] by requiring the derivative of Eq. (26) to vanish at long times, i.e. lim t→∞ ∂ tV (t, 0)ρ 0 = 0. This equation however still requires the knowledge of the full transient propagator. A major simplification arises in DMFT if the system reaches a stationary state becoming time-translational invariant. Then the condition for the impurity density matrix simplifies to (see Appendix (F)) where the self-energyŜ(τ ) depends only on the steadystate propagatorV(τ ) and not on the transient dynam-ics, allowing to compute ρ s in a steady-state DMFT procedure. Equation (31) is analogous to the well known condition for the stationary state of Markovian master equations, with an additional contribution of the non-Markovian bath given by the time-integral of the NCA self-energy. In practice, to solve DMFT/NCA for the stationary state, we solve the Dyson equation (26) for V(t) starting from an initial ansaz for ∆(t − t ), Φ(t ) and ρ s . As an initial condition we usually compute these quantities from the steady-state solution of the singlesite problem. Then we compute the updated stationary density matrix ρ s using Eq. (31) and the updated ∆(t − t ), Φ(t ) from Eqs. (11), (13) and iterate until convergence is reached. We conclude by noting that in principle the stationary state approximation could break down, leading to oscillatory behaviors. It is therefore important to study the stability of the steady state, which is encoded in the retarded Green's function as we discuss more in detail in Sec. V C.

V. DMFT RESULTS FOR A DRIVEN-DISSIPATIVE BOSE-HUBBARD LATTICE
In this section we discuss our results for the driven-dissipative Bose-Hubbard model introduced in Sec. II C, comparing different impurity solvers (NCA and Hubbard-I approximation) and highlighting the effect of introducing fluctuations beyond Gutzwiller mean-field due to the finite lattice connectivity. We start by discussing the properties of the normal phase at low hopping as encoded in its local spectral function (Sec. V A). We then move on to occupation properties of the nonequilibrium normal phase (Sec. V B) from the point of view of the local density and populations of the stationarystate reduced density matrix. In Sec. V C we discuss the finite-frequency instability of the normal phase, leading to the DMFT/NCA phase diagram, and provide a physical interpretation based on hopping-induced dissipation for the large reduction of the ordered phase found in DMFT with respect to the Gutzwiller mean-field result. In Sec. V D we connect the phase transition in our driven-dissipative Bose-Hubbard model to the physics of an array of quantum Van der Pol oscillators, in particular to the onset of many body synchronization and limit cycles and discuss their fate at finite lattice connectivity. Finally, in Sec. V E we discuss the regime of large two-body losses, where Quantum Zeno physics emerges and the qualitative differences between Gutzwiller and DMFT/NCA results appear even more clearly.
Unless stated otherwise, we work in the regime where the interaction strength dominates the dissipation scale, i.e. we fix η/U = 0.02, and study the model as a functions of the pump/loss ratio r and the hopping to interaction ratio J/U . We set U = 5 and ω 0 = 1, although we note that this latter scale only sets the zero of energy and can be eliminated by going to a rotating frame, so it does not play any role in the physics.
We introduce a cutoff on the local Hilbert space dim H , whose value will be specified for each result. We solve DMFT for the normal phase, where Φ = 0 and the anomalous (Nambu) Green's function components vanish so that the self-consistent bath only retains Keldysh indexes ∆ αβ . The NCA propagator in the stationary regime V(t) is obtained, as described in Sec. IV B 1, by propagating in time the derivative of the Dyson equation (26) assuming time-translational invariance with an implicit second-order Runge-Kutta scheme [66], a propagation time t max = 10 and a time step dt = 0.004. We note that in the regime under consideration in this work the dynamics of the Dyson equation is dominated by the non-Markovian bath rather than by the two particle losses and therefore a t max = 10 = 1/η is sufficient to reach convergence. Convergence of the implicit Runge-Kutta at each time step is assumed to be reached when 1/ (t)| < 10 −5 , being i the iteration index. The convergence of the DMFT scheme is assessed by checking that being i the index of the DMFT iteration. We have checked that our results essentially do not change by decreasing those thresholds or increasing the Hilbert space cutoff.

A. Spectral Function in the Normal Phase
To characterize the properties of the system, we first focus on the local retarded Green's function defined in Eq. (19).
Since in the normal phase all anomalous (Nambu) Green's function components vanish as well as the average of the order parameter, we have only one independent Nambu component G R (t) = −iθ(t) [b(t), b † (0)] . Its imaginary part defines the local spectral function In Fig. 5 we plot the local spectral function in the low pump regime, r = 0.6, for different values of J/U and compare the DMFT/NCA results with those obtained with Hubbard-I impurity solver and Gutzwiller meanfield. The Gutzwiller mean-field spectral function shows a series of narrow peaks, whose broadening is controlled only by the local dissipation. We remark that in this approach, corresponding to the infinite coordination number limit z = ∞, all properties of the normal phase are independent on the hopping and coincide with the singlesite J = 0 limit. Indeed as we have discussed in Sec. III A, for z = ∞ the only feedback from neighboring sites comes through the order parameter Φ, which vanishes in this the normal phase. DMFT instead is able to capture the effect of coherent hopping processes, resulting in a further broadening of the resonances. This finite hopping correction to the spectral function reflects the fact that the stationary density matrix in the normal phase is not a tensor product of single-site density matrices, as predicted by Gutzwiller, but rather includes correlations among neighboring sites encoded within DMFT in the non-Markovian bath.
A comparison between Hubbard-I and NCA, shows that the former largely underestimates the effect of the bath. Indeed within NCA the sharp peaks of the isolated single site problem are largely broadened already for a moderate value of the hopping rate J/U = 0.2, a trend that further increases for larger values of J/U . At the same time the location of the poles is found to be weakly dependent on the hopping rate and, at least for J/U = 0.2, essentially captured already by Hubbard-I and Gutzwiller mean-field. This difference can be understood by noticing that there are two main sources of resonance broadening within our DMFT approach, one coming from the bare non-Markovian bath ∆(ω), the other coming from the Markovian interacting single-site problem, encoded in the self-energy Σ(ω) (see Sec. III C). Within Hubbard-I the latter is independent of J, and only set by pump and losses. NCA on the other hand accounts for many body scattering channels mediated by the bath and results in an imaginary part of Σ(ω), also scaling with the hopping strength and responsible for the larger broadening. We emphasize that while the main effect of the DMFT bath in this regime is to broaden the resonances, this broadening is not uniform in frequency, i.e. it could not be reproduced by treating the DMFT bath with a Markov approximation. In fact, the selfconsistent condition Eq. (13) implies that the spectrum of the DMFT bath is given by the local spectral function itself, namely it comes with a rich multi-peak structure in frequency which prevents the use of a simple Markovian approximation. Overall the spectral function in this low-drive, low-hopping normal region is very reminiscent of an equilibrium Bose Hubbard model in the Mott insulating phase [127], with Hubbard bands, describing doublon/holons and multiparticle excitations, which are partially filled by incoherent pump and dissipation. As we show in the next section, increasing the pump strength reveals a spectral feature which is instead unique to interacting driven-dissipative systems.

Negative Density of States
We now discuss how the spectral features of the normal phase evolve upon increasing the strength of the drive/loss ratio r. While in the low pump regime all the peaks of the spectral function are positive, see Fig. 5, a novel effect appears at large drives. Above a threshold r ndos the lowest Hubbard band flips sign and a region of Negative Density of States (NDoS) appears in a positive frequency range [139].
We show this in Fig. 6 where we plot the spectral function obtained within NCA/DMFT for different values of drive/loss ratio r, at fixed J/U = 0.2. The region of NDoS extends up to ω = Ω 0 , a frequency at which the with γ > 0, while for ω > Ω 0 the conventional positive sign is recovered. As we show in figure 6 the spectral range of NDoS increases with the drive r and so does the frequency Ω 0 (r). We stress that a negative spectral function at positive frequency is a genuine nonequilibrium phenomenon that cannot happen for closed systems in thermal equilibrium [113]. It has a direct physical consequences on the response of the system to a weak local coherent drive oscillating at frequency ω, Indeed for an open system the power absorbed from the perturbation, defined as [140-142]Ẇ = Trρ(t)V can be written within linear response theory as (see appendix D) This expression highlights how the spectral function at frequency ω controls the power absorbed by the system under an external drive. A change in sign of this quantity, i.e. a negative absorbed power, signals the onset of energy emission and gain, a condition which is generally associated to optical amplification and lasing [143][144][145][146].
As we are going to discuss in Sec V C the NDoS effect and the frequency Ω 0 will play a crucial role in the nonequilibrium phase transition from the normal to the superfluid phase. We emphasize that the NDoS effect arises already in the single-site problem, i.e for J = 0 in our model, above a threshold pump r ndos which depends on the strength of Kerr nonlinearity, as discussed in Ref. [113]. As a result, it naturally appears at large drive in the normal-phase spectral function of our lattice model calculated within Gutzwiller mean-field theory as well as DMFT/Hubbard-I, both built out of the exact solution of the single site problem.
We now discuss the dependence of the NDoS effect from the hopping J. Clearly such a question goes beyond Gutzwiller mean-field theory, which as we stressed cannot capture any effect due to coherent hopping within the normal phase. In figure 7 we plot the spectral function obtained with DMFT/Hubbard-I and NCA, for increasing values of J/U and compare with the results obtained from Gutzwiller. We find that the NCA spectral function is strongly affected by the hopping, which broadens the sharp high-energy peaks and decreases the strength of the negative peak around Ω 0 , up to a value of J/U 0.8 at which this peak turns back to positive, washing away the NDoS effect. In other words NCA is able to capture a renormalization of the scale Ω 0 from the hopping J. This is surprising at first since J is a purely coherent energy scale while we have seen in Figure 6 that the strength of the peaks and the NDoS is controlled by the dissipative scales, i.e. the pump to loss ratio. We interpret this effect as a first example of hopping-induced losses, a mechanism that is unique to open quantum systems and that plays a key role in the physics of our model. Importantly, this effect is completely missed by the simple impurity solver Hubbard-I, whose spectral function, also shown in Figure 7, changes very little with respect to the Gutzwiller mean-field one [147].
To summarize, we have seen that changing drive and hopping largely affects the spectral properties of the normal phase. In particular we have identified for positive frequencies 0 < ω < Ω 0 (r, J/U ) a region of NDoS emerging above a threshold drive strength r > r ndos (J/U ). Both these quantities depend within DMFT/NCA from the hopping-to-interaction ratio J/U , an effect which is completely missed by Gutzwiller mean field as well as by Hubbard-I. As we are going to discuss in Sec. V C, these dependencies of the critical frequency Ω 0 and of the threshold drive r ndos from the hopping strength J will have direct consequences on the phase diagram of the model.

B. Steady State Local Density Matrix and Population Inversion
We now discuss the occupation properties of the stationary state distribution in the normal phase. For a lattice problem computing the full many-body density matrix can be done only for very small systems. Nevertheless, within our DMFT/NCA approach, describing the thermodynamic limit of infinitely-many sites, we can compute the reduced steady-state density matrix of a given site of the lattice, say site i = 0, obtained by performing a partial trace on all other sites, namely ρ s = tr j =0 ρ latt,s . This corresponds to the steady-state density matrix of the DMFT self-consistent quantum impurity model (9) and thus of the non-Markovian mapV (25) and it obeys Eq. (31). This reduced on-site stationary density matrix allows to study the change of the local populations of bosons due to hopping processes, which is completely missed by Gutzwiller mean-field theory. Also, for open systems these hopping processes enable new, effective dissipative channels. For example a particle can be injected from the Markovian environment on one site, hop to another site and escape the system, rather than just being created and annihilated on the same site. Those processes are captured by our DMFT approach and mimicked by the non-Markovian environment ∆ and unlock interesting new physics which we are going to discuss here.

Local Occupation vs J
From the knowledge of the on-site reduced stationary state density matrix ρ s we can obtain the average local density n = tr b † bρ s . We notice that the local density can be also obtained from the Green's functions, in particular from the Keldysh component at equal times which gives consistently the same result in our NCA approach.
Within DMFT the local density acquires a dependence from the hopping J, which is obviously missing in Gutzwiller mean field. In Fig. 8 we plot the density as a function of J/U at z = 6 and for different values of the drive, normalised to the mean-field value (z = ∞). We see that quite generically the density decreases smoothly upon increasing the hopping within the normal phase, i.e. for J < J c . This can be understood as an interplay of two particle losses and coherent hopping between neighboring sites, an effect that will be further explained by discussing the stationary-state populations in the next section. Interestingly the rate of decrease of the density with hopping changes quite strongly with the strength of the pump r and in particular we notice in Figure 8 that a large drive seems to make the density more pinned to the single-site value.
The result in Fig. 8 turns out to be a specific feature of dissipative lattices with two-particle losses. In fact one can generically prove that for a driven-dissipative Bose-Hubbard model with only single particle losses and single particle drive the stationary state density matrix is independent of any Hamiltonian parameter [148], leading to a density of particles independent of J (although not necessarily integer, as it would be in the equilibrium Mott ground state of the Bose-Hubbard model ) and only set by drive/loss balance. In App. B we show that this effect is correctly captured by our DMFT/NCA approach, a highly non-trivial benchmark for its validity.

Steady-State Populations and Population Inversion
In this section we discuss the effect of coherent hopping processes on the steady state reduced density matrix, which as we show exhibits richer physics than the local occupancy. In the normal phase this quantity is diagonal in the basis of Fock states with n bosons per site and the steady state populations ρ n are shown in Fig. 9 for different hopping values, z = 6 and r = 3.
First we observe that the single-site model, corresponding to J = 0, shows a non-monotonic behavior of the populations as a function of the number of bosons per site, for drive/loss ratio r > 1 and any value of the Kerr nonlinearity U (which in fact does not affect the stationary state as it is has been long known [111]). This population inversion at J = 0 appears clearly in Fig. 9, where the probability of finding n bosons per site is maximum at n = 1 despite the fact that a finite bosonic occupation costs energy E n = ω 0 n + U n 2 /2 ∼ U and should be therefore thermodynamically suppressed.
Increasing the hopping changes the populations at low occupancy while leaving essentially unaffected the tail at large n. In particular, the coherent hopping from and into the neighboring sites increases the probability of having an empty site at expenses of finite occupation. This is a genuine feature of our dissipative many body lattice problem with local two body losses: starting from a state with average filling n ∼ 1, hopping processes towards neighboring sites creates double occupations which escape at a rate η, reducing the total occupation. This trend goes on upon further increasing J, ultimately suppressing the population inversion above a threshold hopping. This mechanism also explains more in detail the observed overall decrease of average occupation with J, Fig. 8, which we already discussed in the previous section. An interesting question concerns the relation between the NDoS effect discussed in Sec. V A and the population inversion in the reduced stationary density matrix. In closed quantum systems described by unitary evolution the two concepts are directly related, namely a NDoS could only emerge in presence of a inversion of populations where higher energy states are more occupied than lower energy ones. For open quantum systems the situation is more subtle and the two concepts are not in one-to-one correspondence [113]. In Figure 10 we plot the behavior of the threshold for population inversion r inv and for NDoS r ndos as a function of J/U . We notice that those thresholds are independent from the hopping within Gutzwiller mean field (see dashed lines, which coincide with the J = 0 values of DMFT) while they are substantially renormalized in DMFT. In particular the two scales r ndos < r inv further deviates from each other as the hopping is increased. We note that r inv increases monotonically with the hopping strength J in DMFT. Based on closed systems arguments, this hopping-induced suppression of population inversion would suggest that the NDoS is also always suppressed by hopping, as for example Fig. 7 shows. Surprisingly, this is not always the case. Figure 10 shows that r ndos has a non-monotonic behavior with the hopping rate J, namely its behavior changes from small and large hopping values. While for large values of J the NDoS threshold r ndos indeed increases following the behavior of the r inv threshold, as expected from closed system arguments, for small values of J it is actually reduced below the J = 0 threshold, corresponding to the single site. Namely, for small values of J, the non-Markovian bath ∆ actually generates a NDoS, even in a regime where the single site model would not present any signature of this effect. This is a unique feature of dissipative quantum systems, where an NDoS can be generated even in absence of population inversion. For Markovian systems it has been shown that this can be traced back to the structure of excitations on top of the stationary state, which come with characteristic complex weights, leading to anti-lorentzian lineshapes [113].

C. Finite Frequency Instability of the Normal Phase
In this section we discuss how the peculiar spectral and occupation properties of the normal phase contribute to an instability towards a spontaneous breaking of U (1) symmetry. We show that the conventional static superfluid transition of the equilibrium Bose Hubbard model as a function of the hopping to interaction ratio J/U is pushed to finite frequency as a result of drive and dissipation, leading to an order parameter oscillating in time. We emphasize the role of the NDoS for the onset of the phase transition and compare the DMFT/NCA and Gutzwiller phase boundaries. We argue that the effect of finite-connectivity fluctuations is not only quantitative, but rather underlines a qualitatively new physical mechanism for the onset of an ordered phase in open quantum lattices with two-body losses, which cannot be simply interpreted as the destruction of an ordered phase by thermal fluctuations in an effective equilibrium problem.

DMFT Phase Boundary
Within our DMFT approach, we can derive an equation for the phase boundary separating the normal and the broken symmetry phases. We assume to be in the early symmetry-broken phase, where the order parameter Φ(t) = b(t) has just formed and it is small. This implies a small external field Φ eff (t) (11) in the DMFT effective action (9). We also assume to be in a stationary regime at long times, such that two point correlators depend only on time differences, and move to Fourier space. The average value of the bosonic field Φ(ω) ≡ b(ω) is, to linear order in Φ eff where we used the fact that b(ω) Φ eff =0 = 0. A key point now is that at finite z the effective field Φ eff (ω) in DMFT has two contributions, one from the local order parameter itself and the other from neighboring sites encoded in the non-Markovian bath, see Eq. (11) which now reads (using Plugging (38) into (37) and using the DMFT selfconsistency on the Bethe lattice (13) one finally gets The critical coupling J c and critical frequency Ω c needed for a self-consistent broken-symmetry solution, Φ(Ω c ) = 0 corresponding to an order parameter whose phase oscillates in time b(t) ∼ e −iΩct for J > J c , are given by Equation (40), which to the best of our knowledge is an original result of this paper, is generic for bosonic DMFT theories on the Bethe lattice and it holds also for equilibrium problems. Its solution, leading to the phase boundary in Figures 2 and 10, strongly depends on the driven-dissipative nature of the problem, as we are going to discuss now. First, Eq. (40) has real and imaginary parts, which both need to vanish simultaneously, resulting in the two conditions We remark that there is another solution possible, where ImG R (Ω c , J c ) = 0, but this is never realized in our simulations. In thermal equilibrium the first condition Eq. (41) can be only satisfied at zero frequency, where fluctuation-dissipation theorem constraints the imaginary part of a bosonic retarded Green's function to vanish, thus allowing for static symmetry breaking patterns (as in equilibrium superfluids). Far from equilibrium this does not need to be the case [109,149] and indeed we have seen that the normal phase shows, above a threshold drive r ndos , a spectral function vanishing at a positive frequency, corresponding to the formation of a NDoS and the onset of gain in the system. The critical frequency Ω c solving Eq. (41) corresponds therefore to the frequency at which the local spectral function of the normal phase vanishes for a critical value of hopping J c determined by jointly solving Eq. (42). The energy scale Ω 0 for the NDoS is therefore a precursor of the mode that will become unstable at the transition. We conclude that the NDoS effect discussed in previous section is a key, necessary condition for a phase transition into the nonequilibrium superfluid phase. This is clearly shown in Figure 10, where we plot the threshold pump r ndos for NDoS and the critical drive r c obtained from solving Eq. (40) with DMFT/NCA as a function of J/U . We see that generically r ndos < r c , namely the system first develops gain and then becomes truly unstable towards U (1) symmetry breaking. On the other hand, from Figure 10 we see that one can obtain an instability even in absence of population inversion.

Role of finite-connectivity fluctuations and comparison with Gutzwiller
We now go back to an important aspect mentioned at the beginning of our paper, namely the large renormalisation to the phase boundary obtained within DMFT/NCA upon decreasing the connectivity z. As we show in Figure 2 and Figure 11, fluctuations due to the finite-number of neighbors shift the phase boundary towards larger values of the hopping J and pump/loss ratio r. We now provide a physical interpretation of this effect based on the properties of the normal phase discussed so far.
We start by considering the condition for the normal phase instability obtained within Gutzwiller mean-field theory, which corresponds to the z → ∞ limit of the Eq. (40). In this limit the problem reduces to a quantum single site, while the feedback from neighboring sites is treated at a purely classical level (See Sec. III A), in terms of a self-consistent coherent field which reads Φ eff = J c b(t) ∼ Je iΩct near the instability. As such, if we repeat the argument of the previous section we can obtain a condition for the Gutzwiller phase boundary, We further plots the thresholds for NDoS obtained within the two approaches. We see that fluctuations due to finite connectivity reduce the broken symmetry phase, pushing it towards higher values of hopping and drive. We interpret this effect as a signature of hopping-induced losses, which reduce the local gain and prevents the system to become unstable at finite frequency. Parameters :dimH = 10 which reads where the first term is the effective field contribution and G R 0 (ω) is the retarded Green's function of the isolated single-site problem and is therefore independent from the hopping. Within Gutzwiller mean-field theory the hopping J has only the role of triggering the symmetry breaking, through the self-consistent field Φ eff , while the onset of gain is controlled by the pump to loss ratio r. Indeed we known that the single site spectral function develops a NDoS above a constant threshold pump r ndos (see gray dashed line in Figure 11). The feedback from the neighboring sites acts as a seed for a single site which is on the verge of energy emission (negative absorbed power at Ω c , see Eq. (35)) and leads, above a threshold hopping J c shown in Figure 11, to amplification of the local coherent field at frequency Ω c and a spontaneous breaking of U (1) and time-translation symmetry. Interestingly, the Gutzwiller phase boundary is very close to the line r ndos (see Figure 11) suggesting that at large hopping as soon as the system develops gain the symmetry breaking occurs.
DMFT/NCA on the other hand accounts for a more subtle effect of neighboring sites, which are encoded in the non-Markovian quantum bath, in addition to the classical coherent field. As we know this provides an increased effective dissipation, due to hopping processes from lossy neighboring sites, which is responsible for wiping out the NDoS region in the local spectral function, a necessary condition for the onset of the instability. This is clearly shown in Figure 10 where for large values of J/U the threshold for NDoS eventually becomes larger than the Gutzwiller mean-field one (see gray dashed line in Figure 10), leaving a normal phase which would be superfluid at z = ∞. These hopping-induced losses provide therefore a clear physical mechanism behind the finiteconnectivity renormalization of the phase boundary. The picture that emerges from DMFT is the one of a single site on the verge of energy emission, coupled to an oscillating seed field which would favor optical amplification and embedded in a non-Markovian bath which is instead able to absorb part of the emitted power from the system, thus reducing the effective gain and requiring a stronger value of pump to trigger the instability. This mechanism is presumably also effective at intermediate value of the hopping where, as we see in Figure 11, the DMFT/NCA threshold for NDoS remains well below the critical drive responsible for the true many-body instability at intermediate coupling, while approaching it at large hopping.
We conclude that this hopping-induced dissipation is a qualitatively new mechanism for the destruction of an ordered phase, which is unique to open systems settings and one of the hallmark of our DMFT/NCA approach. Quite interestingly this mechanism is not only completely missed by the Gutzwiller mean field approach, but also by a perturbative solver such as the Hubbard-I approximation. In fact we have discussed (see in figure 7) how the NDoS is changed very little within this scheme. In appendix E we also show that the Hubbard-I phase diagram, obtained using the DMFT equation for the critical point (40), still reduces to the Gutzwiller mean-field one. This further highlights the non-perturbative nature of the fluctuations responsible for the observed renormalization of the NDoS and the importance of using a self-consistent scheme such as our NCA approach.

Thermal vs. Non-Thermal Origin of Finite Connectivity Fluctuations
A different perspective on the role of finite-connectivity fluctuations on the phase boundary can be obtained by looking at the occupation of single particle modes at finite frequency, describing excitations on top of the stationary state and their effective thermal character. In fact, in an effective equilibrium picture, one could expect heating to provide an efficient mechanism for the reduction of a broken symmetry phase. To investigate this physics we look at the Keldysh Green's function, defined in Eq (36), which heuristically describes the fluctuations of the observable b. If the system was in true thermal equilibrium, the quantum fluctuation-dissipation theorem (FDT) would constrain the Keldysh and the retarded components to obey the relation [126] FIG. 12. Change of the effective temperature T eff at the critical point (Ωc, Jc) as a function of lattice connectivity z, with respect to its mean-field value at z = ∞, for different values of the pump to loss ratio r. We see that the relative change in T eff remains of the order of few percents up to connectivity z = 50 and slightly increases up to 20% for the lowest connectivity z = 20. This should be compared to the DMFT renormalization of the critical hopping for similar values of drive and connectivity, which is instead much stronger as shown in Figure 2. In particular for z = 30 and r = 3 the critical hopping is pushed to infinity corresponding to a complete destruction of the normal phase. Parameters : dimH = 14.
where T is the system temperature. At low frequency or high temperatures, ω T , one has F eq (ω) ∼ T /ω. In a non-equilibrium system on the contrary there is no well-defined temperature and the FDT does not hold in general. Nonetheless, it is useful to use the left-hand side of the FDT relation in Eq. (45) to define an effective distribution function and to study its frequency dependence. Within the DMFT and Gutzwiller normal phases, for pump above the threshold for NDoS r > r ndos , the spectral function A(ω) vanishes at frequency Ω 0 , with linear corrections (see Eq. (34 ). On the other hand we find that the Keldysh component has a finite non-zero value at Ω 0 , which gives a distribution function of pseudo-equilibrium form at least for the modes around Ω 0 From this expression we can therefore identify an effective temperature T eff , which emerges quite ubiquitously in nonequilibrium quantum systems [52,[150][151][152][153][154]. At z = ∞, when the normal phase is described as a collection of independent sites, an effective temperature emerges due to the interplay of local drive and Kerr interaction [113]. Within DMFT/NCA one could in principle expect important corrections due to the non-Markvovian bath, as we have seen for the NDoS. To assess this point we sit at the phase boundary, J c (r), choose the corresponding critical frequency Ω 0 (J c , r) = Ω c and plot the behavior of T eff , scaled with respect to the z = ∞ value, as a function of the lattice connectivity z and for different values of the pump-to-loss ratio r. We see that at large values of r the effective temperature slightly increases upon decreasing z, while at smaller pumps, i.e. r = 3, it shows a weak non-monotonic behavior with z. Overall, the relative variation of T eff with respect to the Gutzwiller mean-field value remains rather moderate. On the other hand, in the same range of variation of z, the DMFT critical hopping shows instead strong renormalizations. This is particularly true for the small drive regime, where already for r = 3 and z = 30 the ordered phase is completely washed out. This observation suggests that effective thermal fluctuations and heating are not enough to explain the large renormalization of the phase diagram observed in DMFT, which is instead mainly driven by the reduction of local gain by virtual hopping processes, through the mechanism of hoppinginduced losses.

D. Nonequilibrium Superfluidity, Lasing and Many Body Synchronization of Van der Pol Arrays
We now discuss how the nonequilibrium superfluid transition in our Bose-Hubbard model is connected to other dynamical phenomena associated with breaking of time-translation symmetry, such as nonequilibrium Bose-Einstein condensation, lasing and synchronization of Van der Pol oscillators. To appreciate this point it is useful to start from the semiclassical limit of our model (See Appendix C) where non-linear quantum fields contributions) are disregarded. The dynamics of the bosonic field at site j takes the form of a Langevin equation with an effective frequency and damping terms which depend non-linearly on the field itself, i.e.ω 0 (b jcl ) = ω 0 + U |b jcl | 2 /2 and γ(b jcl ) = f /2 − η|b jcl | 2 /2, and where ξ j (t) is a zero average white noise ξ i (t)ξ j (t ) = f δ(t − t )δ ij . In the continuum limit this reduces to a complex Gross-Pitaevski equation [155,156] with pump and non-linear losses which describes a variety of non-linear phenomena from exciton polariton condensates to multi-mode lasers [157]. In absence of any noise the spatially uniform stationary state admits a stable limit cycle, i.e. β(t) = |β|e −iω vdp t for r > 0 and any J with frequency ω vdp ω vdp = ω 0 − J + U |β| 2 /2 (48) and amplitude set by the incoherent drive, |β| = √ r.Equation (48) is therefore the semiclassical version of our condition for finite-frequency instability in the normal phase. A major difference exists between the lasing threshold and our case, namely that the threshold for the onset of gain and the one for development of full coherence is well separated, while it coincides in the usual lasing regime [149]. This can be seen easily by looking at the Green's function which develops a NDoS as soon as the system becomes unstable.
The equation above describe also an array of coupled classical Van der Pol oscillators. Our driven-dissipative Bose-Hubbard model can be therefore also seen as a quantum many-body version of the VdP array. From this perspective, the onset of finite-frequency oscillations at J c (r) described in the last section can be seen as a signature of a quantum synchronization [70][71][72][73][74][75][76][77][78][79][80][81][82], where above a certain coupling J all quantum VdP oscillators enter into a collective limit-cycle phase. As we are going to discuss next, these oscillations share qualitative features with the semiclassical solution at least at large drive values, while deviate significantly for smaller drives where quantum fluctuations are important.
Despite the similarities the semiclassical limit described above is rather different from the large connectivity limit: in the former case one has a deterministic non-linear equation and desynchronization can only appear due to the noise. In the latter instead, corresponding to DMFT, one reduces to a quantum VdP oscillator coupled to a self-consistent field and a quantum bath. This is responsible for the large separation between the onset of gain and the true instability. and drive-to-loss ratio r = 3 within DMFT/NCA. We see that the frequency Ω0 at which the density of states vanishes, decreases with U and eventually disappears for small enough values of U . Parameters: z = 6, dimH = 14.

Critical Frequency versus Drive/Loss Ratio and Kerr
Non-Linearity In this section we discuss the behavior of the critical frequency Ω c , signaling the onset of a quantum synchronized phase, as a function of pump/loss ratio r and Kerr non-linearity U . In Figure 13 we plot this frequency as a function of the incoherent drive amplitude r, both for DMFT (for z = 50) and for Gutzwiller mean field theory (z = ∞). We see that Ω c scales linearly with the drive at large values of r (see fit in Fig. 13), a result which is in agreement with the semiclassical result obtained for the VdP array Eq. (48) where ω vdp /U ∼ r/2.
As the pump is reduced and the number of bosons per site decreases one expects quantum fluctuations to become more important. Indeed we see significant deviations from the semiclassical result at small r, already captured by Gutzwiller mean field but more pronounced for DMFT. Another interesting aspect is the role played by the Kerr non-linearity. In Fig. 14 we plot the local density of states for fixed value of hopping and drive/loss ratio and for different values of the interaction U . A first interesting observation is that the frequency Ω 0 at which NDoS emerges, related to the mode becoming critical at the synchronization transition (see Eq. (43), decreases upon reducing the Kerr non-linearity. A detailed analysis shows Ω c ∼ U . We notice the analogy with the semiclassical result (48), nevertheless we stress that in this picture the critical frequency is proportional to the modulus square of the order parameter | b | 2 = 0, while within DMFT we find Ω c ∼ U at the critical point, where by definition b = 0 . This further highlights the quantum nature of the synchronization transition considered in this work. Indeed we have seen how Ω c is smoothly connected to the frequency Ω 0 where NDoS emerges, a scale that exists already well inside the normal incoherent phase and that is a genuine feature of the quantum impurity model. Furthermore this implies that the Kerr non linearity is crucial in order to push the transition at finite frequency, implying undamped oscillations of the order parameter. We indeed observe that for a related model of all-to-all coupled quantum VdP oscillators recent Gutzwiller mean-field analysis reported a static (first-order) transition in absence of any Kerr non-linearity [70,119]. We conclude by noting that an oscillating phase in the order parameter could be in principle gauged away by going to an appropriate rotating frame. This however could not be done a priori, since as we have seen the critical frequency Ω c itself depends from the many-body physics of the problem. We are considering here only the instability of the normal phase, rather than the dynamics in the full broken symmetry phase which could show a more complex dynamical behavior, whose description goes beyond the scope of the stationary-state oriented approach (see Sec. IV B 1) used here.

Limit Cycles at Finite Connectivity
Limit-cycles emerge ubiquitously within mean field approaches [38,[158][159][160], and indeed even in the present problem the Gutzwiller solution at z = ∞ predicts one. The role of fluctuations on their stability has been discussed before, in particular in the context of a coherently driven anisotropic Heisenberg model with spontaneous decay [161]. There an approach based on self-consistent Mori projection [162] and cluster mean field [51] predicts that limit cycles disappear as the coordination number z is decreased below a threshold value z * which depends on the system parameters.
On the basis of these results, it is particularly interesting to study the fate of our synchronization transition beyond Gutzwiller mean-field theory. In Fig. 15 we plot the behavior of the critical frequency Ω c obtained from DMFT/NCA, as a function of z for different value of r. We observe that finite z corrections tend to reduce the value of Ω c with respect to the mean field value, which nevertheless remains finite down to the lowest value of connectivity at which, for a given value of drive, a synchronization transition exists, consistently with the phase boundary moving to higher values of r for decreasing z (see Fig. 2). In other words we find that, within our treatment, finite connectivity does not destroy the limit cycle phase which is only pushed at higher values of the incoherent drive. Indeed in Fig. 15 we show that for drive r = 7, the highest value that we can numerically access given the constraints on the local Hilbert space truncation, a limit cycle would exists down to z = 20. We therefore expect that at higher drives the synchronized phase would survive down to low connectivity values. The regime of strong drive is however difficult to access within NCA and we leave this question open for future works. We remark nevertheless that our model is different from the one of Ref. [161]. In the present case the existence of a limit cycle phase is tightly related to a U (1) symmetry present in the original Lindblad problem and spontaneously broken at the transition [109].  17. Quantum Zeno regime in the spectral function, increasing two-particle loss rate η at fixed drive rate f /U = rη/U = 0.013. Parameters: tmax = 20, dt = 0.002, dimH = 6. (a) Spectral function at increasing η/J. The first peak corresponds to transitions between states with on-site occupation n = 0 and n = 1, while higher energy peaks are suppressed at large η/J. We fit the first peak as described in the main text in order to extract a measure of its lifetime γ. (b) Partial lifetime γzeno/U of the first peak of the spectral function obtained by removing the zero hopping and zero drive/dissipation contributions. We show that this quantity has a non-monotonic behaviour, analogous to the populations in Fig. (3), manifesting a quantum Zeno behaviour in the lifetime of excitations of the steady-state. (c) γzeno/U depends only on the dissipative scale Γ eff /U , analogously to the occupation probability ρ1 as shown in the inset of Fig. (3).

E. Quantum-Zeno Regime
In the previous sections we have seen that the key feature of DMFT is the ability to capture hopping-induced dissipative processes in the normal phase, which are missed by Gutzwiller mean-field theory and ultimately responsible for the large corrections to the phase boundary due to finite connectivity. In this section we discuss a different parameter regime of our model, corresponding to large two-particle losses η J, which alllows us to highlight even more the qualitative difference between DMFT and Gutzwiller predictions. In the large twoparticle losses limit, and in absence of any external drive, a perturbative analysis has shown that the dissipative dynamics of the Bose-Hubbard model effectively takes place in the subspace with zero and one boson per site, which are dark states of the local dissipator [120]. These emergent hard-core bosons are subject to next-neighbor losses controlled by the scale This effective dissipation is hopping-mediated -namely it is non-zero only at finite hopping J -and remarkably it shows a non-monotonic behaviour as a function of the physical dissipation η: for η U , it increases linearly with η. Instead, when η U , the effective dissipation Γ eff is suppressed by increasing the physical dissipation η. The latter is a quantum-Zeno regime [121,122], in which the coherent hopping dynamics bringing the system out-side of the hard-core bosons subspace is suppressed due to the large coupling to the environment [120]. In absence of an external pump the effective loss rate (49) ultimately leads the system to a zero density state and the Zeno scale 49 manifests itself in the transient dynamics [120,123]. Here instead we demonstrate how, in presence of a small pump, the Quantum Zeno regime emerges in the stationary state properties of the system and how this effect is completely missed by Gutzwiller, even at a qualitative level, while captured by our DMFT/NCA approach.
We consider the regime of η J and a parametrically small pump rate f η. In this case we expect the local occupation of all the states with more than one boson per site to be largely suppressed upon increasing the two-particle losses, since the small residual pump is not sufficient to counterbalance the losses. This is indeed shown in Fig. 16 where we plot the occupation probabilities of the on-site reduced density matrix obtained from DMFT/NCA and see that only the occupation of the states n = 0, 1 remain of order one, while ρ n≥1 is exponentially suppressed. We emphasize that in this case the long-time limit of the problem remains non-trivial within this subspace, due to the interplay between the small residual pump and the coherent dynamics generated by the Hamiltonian. In Figure 3 we have shown that our DMFT/NCA approach is able to capture the emergence of a stationary Zeno regime, through a non-monotonic behavior of the steady-state probability of having exactly one bosons per site, namely ρ 1 , versus η/U , showing a universal collapse when plotted against Γ eff . This is remarkable, since the dissipative scale in Eq. (49) describes particle losses which are non-local in space [120] which a priori could have been expected to be beyond reach for a local approach such as DMFT. Instead, the non-Markovian bath is able to capture this hopping-induced dissipative processes. We emphasize that DMFT goes beyond the effective hard-core model of Ref. 120 in that the full crossover from normal to Zeno phase is captured as η/J is increased. In this sense we can see DMFT as providing a non-perturbative solution of the effective hard-core model in the limit of large lattice connectivity.
The presence of Zeno scale Γ eff in the system not only affects the steady state populations, but is also expected to influence the lifetime of excitations of the steady-state, which can be extracted from the Green functions computed in DMFT. In Fig. 17 (a) we show the behaviour of the spectral function, defined in Eq. (33), for increasing η/J. We work in a regime in which U is large so that different resonances in the spectrum, corresponding to single-particle transitions between states with occupation n and n + 1, are well separated [113]. By increasing η/J all the peaks but the first one decrease in amplitude, as the former are proportional to the occupation of states with n ≥ 2, which decreases with η (see Fig. 16). The amplitude of the lowest frequency peak, instead, increases with η due to the increased probability of being in the subspace with occupation n = 0, 1. We concentrate on this peak in order to extract its width. In order to do so, we fit the peak with a generalized Lorentzian A(ω) = a (ω−ω0)b+1 (ω−ω0) 2 +γ 2 , with fitting parameters a, b, γ, which exactly describes the shape of the peak at zero hopping and which approximates it at small hopping, whose lifetime is given by the fitting parameter γ.
We remark that the width γ of the first peak has three contributions, γ = γ diss + γ coh + γ zeno . The first term comes from dissipation and pump at zero hopping γ diss = γ(J = 0, η, f, U ) and it can be well estimated with perturbation theory [113], which shows it does not depend on η, as the transition n = 0 → n = 1 is unaffected by the two-particle losses at J = 0. The second contribution to the lifetime is given by the hybridization of the energy levels of different sites at finite hopping, which is a purely coherent effect, independent of drive and dissipation. For example, this effect is responsible for the formation of the holon and doublon bands in the local spectral functions of the isolated Bose-Hubbard model, see e.g. [127]. We call this term γ coh = γ(J, η = 0, f = 0, U ).
Finally, we expect a contribution coming from the hopping-induced decay, which we call γ zeno , which is present only for finite hopping and dissipation/drive. In Figure Fig. 17 (b), we show that the latter contribution follows a similar non-monotonic behaviour as Γ eff as a function of η/U , witnessing the presence of quantum-Zeno behaviour in steady-state excitations. Analogously to populations, in Fig. 17 (c) we show that plotting γ zeno /U as a function of Γ eff /U , the data points at different values of J/U all collapse on the same curve, demon-strating that γ zeno only depends on J and η through Γ eff .

VI. CONCLUSIONS
In this paper we formulated the Dynamical Mean Field Theory for bosonic open quantum many-body systems described by a Lindblad master equation on the lattice. This method is based on a systematic expansion in the limit of large lattice connectivity z. Within DMFT, fluctuations due to finite lattice connectivity are treated nonperturbatively through the solution of a self-consistent quantum impurity model, which in our case amounts to a Markovian single-site problem coupled to a coherent field and to a non-Markovian quantum bath mimicking the rest of the lattice. The non-Markovian bath contains the key new ingredient of DMFT, which makes it different from the Gutzwiller mean-field theory, to which it reduces in the z = ∞ limit. In particular, with respect to the former, DMFT is able to capture both coherent and dissipative processes arising from the neighboring sites, the former playing a particularly crucial role in open quantum systems with correlated jump operators.
Using DMFT, together with a non-perturbative bosonic impurity solver based on a super-operator hybridization expansion truncated at the NCA level, we solved a Bose-Hubbard model in presence of two-particle losses and single particle pump, relevant for dissipative ultracold atoms as well as for arrays of superconducting circuits. We have shown that this model features a dissipative phase transition from an incoherent normal phase to a nonequilibrium superfluid, which occurs above a critical hopping or pump strength.
Within DMFT/NCA the phase boundary is strongly renormalized with respect to the Gutzwiller one, and pushed towards higher values of the couplings, leading to an increase of the normal phase. We have identified a new mechanism for this reduction of the broken symmetry phase which is unique to open quantum systems and arises from the suppression of local gain due to hoppinginduced losses. This is a key feature brought forth by NCA, which treats the non-Markovian DMFT bath nonperturbatively, and it is instead missed by a simpler perturbative DMFT solver based on the Hubbard-I approximation. We have further discussed how the increased effective dissipation due to the finite number of lossy neighbors affects all the unusual properties of the DMFT/NCA normal phase, from the renormalization of the gain (NDOS) threshold and steady-state populations due to hopping, to the emergence of a stationary-state quantum Zeno regime for large two-body losses. These effects are all qualitatively missed by Gutzwiller meanfield theory.
Finally, we have shown that the transition into the nonequilibrium superfluid phase occurs at finitefrequency, corresponding to a local order parameter that oscillates in time at a frequency Ω * . Within DMFT/NCA this scale depends on both coherent and dissipative cou-plings, another important aspect which is missed by Gutzwiller mean-field theory. Drawing from the physics of quantum VdP oscillators we interpreted this phenomenon as the onset of quantum many body synchronization and limit cycles at finite lattice connectivity.
Our DMFT holds the promise to be applied to a variety of driven-dissipative quantum many body problems. Different bosonic models or driving schemes can be considered and readily studied with the NCA approach developed here, such as the recently introduced quadratically driven Kerr resonator [40], the coherently driven case explored in the context of quantum bistability [47,54,105,107,108] or models relevant for the dissipatively stabilised Mott insulators of photons [16]. Interesting directions for the future involve the development of other quantum impurity solvers, made possible by the Markovian structure of the quantum impurity problem. In particular, one could go beyond the Non-Crossing Approximation by including higher-order diagrams or resumming the full hybridization expansion [69] with real-time diagrammatic Monte Carlo methods [130,132]. Alternatively, one could map back the non-Markovian Keldysh action into a Lindbladian problem with a finite number of bath levels [86] which could then be solved by exact diagonalisation [163], quantum trajectories, through a matrix product operator representations of the density matrix [164,165] or using an extension of the numerical renormalization group [166]. Further extensions of this approach could include the development of a cluster version of DMFT building upon recent developments for Markovian problems [51].
We present here a sketch of the derivation of the DMFT action and self-consistency conditions, using the quantum cavity method and following and extending Ref [68] to the open case. The main idea is to single out a given site of the lattice, i = 0 in the following, and to write down its effective Keldysh action obtained after integrating out all the other sites.
This effective single-site action describes, in principle, all the effects on site i = 0 due to the coupling to the other sites by the hopping J, which, for our assumption of local jump operators and interactions (Sec. II), is the only term responsible for the coupling between different sites. To proceed, we notice that the full action Eq. (7) can be divided in three therms, S = S 0 + δS + S (0) cav respectively S 0 = S[b 0 , b 0 ] containing only the terms in Eqns. (7)(8) involving fields at site i = 0, a term δS containing the hopping terms to the 2z fields (b j , b j ) at the neighboring sites and a term S (0) cav describing a lattice with a cavity, namely including all the degrees of freedom except those at site i = 0, see Figure 18. For an interacting many-body problem on a finite dimensional lattice integrating over the neighboring sites can only be done formally, as a cumulant expansion in δS, and leads to an effective action containing arbitrary powers of the local fieldsb 0 , b 0 , with coefficients given by the multipoints correlation functions of the cavity problem [62]. In the large connectivity limit, z 1, one can formally organize this expansion in power of 1/z and obtain where, in order to allow for condensed phases, we have introduced the Nambu-spinor notation The coefficients entering the effective action (A3) are related to quantum averages over the cavity action S The final step, in order to obtain the self-consistent conditions, is to relate the average of bosonic fields on the cavity action to those evaluated to the effective action (A3). For what concerns the non-Markovian bath in Eq. (A5) this depends on the lattice geometry: it becomes particularly transparent for a Bethe lattice, a lattice with no loops such that once a cavity is created two neighbors j, k are completely disconnected, where one gets where we have used the property of the Bethe lattice in the first equality, the fact that in thermodynamic limit the local property of the cavity action and the original action must be the same in the second equality and translational invariance in the last step. Plugging this into Eq. (A5) gives the self-consistency condition (13). We notice that simular arguments can be used for a different lattice, to relate averages over cavity and effective action, the only difference would be a more complicated self-consistency relation between bath and local Green's function [62]. Finally, the average of the bosonic field taken on the cavity action can be related to the one on the effective action (A3) as [68] b † jα (0) which can be plugged in Eq. (A4) to give the second selfconsistent condition (11). In this appendix we report a benchmark our DMFT/NCA approach for driven-dissipative many-body master equations defined by Eqs. (2) and (3). We consider a different model from the main text, namely the driven-dissipative Bose-Hubbard model with singleparticle losses and single-particle drive. This is specified by the same Bose-Hubbard Hamiltonian Eqs. (3) and (4) as in the main text, but with the jump operatorŝ where there is single-particle losses instead of the twoparticle losses of the main text.
One can prove [148] that the stationary state density matrix of this model is independent from any Hamiltonian parameter, thus, for example, the on-site occupation is constant with the hopping rate J. Obtaining this constant occupation with J is a highly non-trivial benchmark for our DMFT/NCA approach. Figure 19 shows that this behavior is correctly reproduced by our approach. Small deviations from constant occupation show that this property is not exactly enforced by our numerical scheme and thus is a good test to validate our approach. Deviations from constant occupation are mostly a result of local Hilbert space truncation, which become more severe increasing the drive, but they are reduced by increasing the cutoff dim H . In this section we write down the Keldysh fieldtheory associated to the Lindblad master equation for our driven-dissipative Bose-Hubbard model in Eq. (3)(4)(5). This is done by writing a coherent path-integral representation of the trace over the density matrix [52,126] Z = Trρ(t) = i Db icl Db icl Db iq Db iq exp (iS) in terms of bosonic fields in the classical-quantum basis b jcl/q = (b j+ ± b j− ) / √ 2, where the Keldysh action reads where the first line describes the local non-interacting contribution, including the dissipative couplings χ ± = (γ p ± γ l ) /2 which in the case considered in the main text, where γ l = 0 and γ p = f = rη, reduce to χ ± = f /2. The second line includes non-linearities due to interaction and two-body losses while the last term account for the hopping between neighboring sites. The semiclassical limit is obtained by disregarding terms with more than two quantum fields, which leaves the action at most quadratic in the quantum fields b jq [52]. Those quadratic terms can be decoupled using a classical stochastic field [126] which allow to write the action as In this appendix we show that a direct consequence of the NDoS effect is that the average power absorbed from a weak coherent drive becomes negative, indicating the onset of gain and energy emission. We consider our lattice model in presence of a time dependent perturbation of the Hamiltonian describing a weak coherent drive with frequency ω where H is the Bose Hubbard Hamiltonian defined in Eq. (3)(4). The evolution of the density matrix of the system is described by the Lindblad master equation in presence of H(t) (D2) We notice that the time derivative of the internal energy E(t) = H(t) in an open quantum system has two contributions, respectively coming from the derivative of the system Hamiltonian and from the time derivative of the density matriẋ E(t) = Trρ(t)Ḣ(t) + Trρ(t)H(t) (D3) the latter being identified as heat flow which would be zero in a closed system evolving with unitary dynamics. We are interested in the first term which measures the absorbed power from the external perturbation [142] and can be written aṡ Within linear response theory the average of the bosonic field can be written as . Plugging this expression into Eq. (D4) we obtain the average absorbed powerẆ written in terms of the spectral function A ij (ω) = (−1/π)ImG R ij (ω). For a localized perturbation, v i = v 0 δ i,0 we obtain the result given in the main text, showing how a NDoS implies a negative absorption rate, namely the system transfer some of the energy from the perturbation to the emitted radiation.
where we have performed a change of variable in the convolution integral. The NCA self-energyŜ(t, t − τ ) in the equation above depends on the hybridization function of the bath ∆ (see Eq. 27), which is itself related to the local Green's function of the system by the DMFT self-consistency. Therefore the assumption that for t → ∞ the system reaches a unique steady-state, forgetting about its initial condition, translates within DMFT in a requirement on the behavior of the NCA self-energy for large values of its argument. In particu-larŜ(t, t − τ ) must vanish for τ ∼ t → ∞, restricting the time convolution in the Dyson equation to times of the order τ ∼ τ * t. We can therefore safely take the long time limit t → ∞ in Eq. (F1), use the fact that in this limit the self-energy becomes translational invariant, S(t, t − τ ) →Ŝ(τ ), and the stateV(t − τ, 0)ρ 0 approaches the stationary state ρ s , and finally send the scale τ * to infinity, yielding the eigenvalue condition This condition depends only on the stationary state prop-agatorV(τ ) through the NCA self-energyŜ(τ ), allowing to compute the stationary density matrix ρ s in the steady-state DMFT procedure.