Alexandre Le Boité

Giuliano Orso

Cristiano Ciuti

Steady-state phases and tunneling-induced instabilities in the driven-dissipative Bose-Hubbard model

Keywords: numbers: 42.50.Ar, 03.75.Lm, 42.50.Pq, 71.36.+c

We determine the steady-state phases of a driven-dissipative Bose-Hubbard model, describing, e.g., an array of coherently pumped nonlinear cavities with a finite photon lifetime. Within a meanfield master equation approach using exact quantum solutions for the one-site problem, we show that the system exhibits a tunneling-induced transition between monostable and bistable phases. We characterize the corresponding quantum correlations, highlighting the essential differences with respect to the equilibrium case. We also find collective excitations with a flat energy-momentum dispersion over the entire Brillouin zone that trigger modulational instabilities at specific wavevectors.

In recent years, the interest in the physics of quantum fluids of light in systems with effective photon-photon interactions has triggered many exciting investigations [1]. Some of the most remarkable features of quantum fluids, such as superfluid propagation [2,3] or generation of topological excitations [4][5][6][7] have been observed in experiments with solid-state microcavities. With the dramatic experimental advances in solid-state cavity and circuit quantum electrodynamics (QED), a considerable interest is growing on the physics of controlled arrays of nonlinear cavity resonators, which can be now explored in state-of-art systems [8,[START_REF]The Physics of Semiconductor Microcavities[END_REF]. This opens the way to the implementation of non-equilibrium lattice models of interacting bosons, particularly when effective on-site photonphoton interactions are large enough to enter the strongly correlated regime [START_REF] Imamoglǔ | [END_REF][11][12][13]. In this kind of systems, it is possible to realize the celebrated Bose-Hubbard model [14] for photons or polaritons. Since the first theoretical proposal for implementing this model in optical systems [15][16][17], early works have been focused on phenomena close to the equilibrium Mott insulator-Superfluid quantum phase transition [18,19]. Strongly-non equilibrium effects have been addressed only more recently [20][21][22][23] particularly in the interesting driven-dissipative regime where the cavity resonators are excited by a coherent pump which competes with the cavity dissipation processes. In such non-equilibrium conditions, these open systems are driven into steady-state phases whose collective excitations can be extremely different from the equilibrium case. However, to the best of our knowledge, very little is known so far on these important properties for the non-equilibrium Bose-Hubbard model.

In this Letter, we present comprehensive results for the steady-state phases and excitations of the drivendissipative Bose-Hubbard model in the case of homogeneous coherent pumping. The steady-state density matrix and expectation values of the relevant observables have been calculated with an efficient mean-field approach, based on exact analytical quantum optical so-lutions of the single-cavity problem. A rich diagram is shown with multiple steady-state phases, whose stability and complex energy excitation spectrum have been studied through a linearization of the Lindblad master equation around the stationary solutions. We unveil the existence of a purely imaginary excitation branch which can trigger modulational instabilities at specific wavevectors.

We consider a driven-dissipative Bose-Hubbard model under homogeneous coherent pumping describing a bidimensional square lattice of cavity resonators. In a frame rotating at the pump frequency ω p , the system is described by the following Hamiltonian [1]:

H = - J z <i,j> b † i b j - N i ∆ωb † i b i + U 2 b † i b † i b i b i + F b † i + F * b i
(1) where b † i creates a boson on site i, J > 0 is the tunneling strength, and z = 4 is the coordination number. < i, j > indicates that tunneling is possible only between first-neighbors. U > 0 represents the effective on-site repulsion, F is the amplitude of the incident laser field and ∆ω = ω p -ω c is the frequency detuning of the pump with respect to the cavity mode. The dynamics of the many-body density matrix ρ(t) is described in terms of the Lindblad master equation:

i∂ t ρ = [H, ρ] + iγ 2 N i 2b i ρb † i -b † i b i ρ -ρb † i b i ( 2
)
where γ is the dissipation rate. While for equilibrium quantum gases, the chemical potential µ is a key quantity, in this non-equilibrium model the steady-state phases depend instead on the pump parameters F and ∆ω, which compete with γ.

Given the success of mean-field theories in the investigation of the equilibrium Bose-Hubbard physics, it is a legitimate starting point for the study of its nonequilibrium version. The mean-field approximation is ob-

H mf = -∆ωb † b + U 2 b † b † bb + (F -J b )b † + (F * -J b * )b
(3) where the value of b has to be determined selfconsistently. The problem of a single cavity has been studied by Drummond and Walls [24], who obtained analytical expressions for the bosonic coherence b and the photons distribution functions via a generalized Prepresentation for the density matrix. By replacing F with F -J b in these exact expressions we find the self-

consistent formula b = (F -J b ) ∆ω + iγ/2 × F (1 + c, c * , 8| F -J b U | 2 ) F (c, c * , 8| F -J b U | 2 ) (4)
for the bosonic coherence. The mean photon density and the other diagonal correlation functions can then be easily extracted from the general expression:

(b † ) j (b) j = 2(F -J b ) U 2j × Γ(c)Γ(c * ) Γ(c + j)Γ(c * + j) × F (j + c, j + c * , 8|F/U | 2 ) F (c, c * , 8|F/U | 2 ) ( 5
)
with c = 2(-∆ω -iγ/2)/U and the hypergeometric function

F (c, d, z) = ∞ n Γ(c)Γ(d) Γ(c+n)Γ(d+n) z n
n! , Γ being the gamma special function.

All the properties of the steady-states are therefore determined by the self-consistent solutions of Eq.( 4) which we have calculated numerically. Due to the presence of the tunneling term J, multiple solutions do appear in certain region of parameters space. We investigate their stability through a linearization of the Lindblad master equation around each steady-state solution as described later in the Letter. In Fig. 1, we present a diagram showing the number of stable steady-state solutions as a function of the tunneling and the on-site interaction in units of the detuning ∆ω > 0 and for a representative set of parameters (see caption). We see that there are regions with 1 or 2 stable solutions, but also regions with no stable homogeneous solution (shown in red). Notice however that, within our mean field approach, we have direct access only to spatially uniform solutions, where all the cavity sites are equivalent.

Interestingly, we find that the bistability induced by the coupling between the cavities also appears when the pump frequency is red-detuned with respect to the cavity mode (not shown), in sharp contrast with the case of an isolated cavity. We see in Fig. 1 that the boundary between monostable and bistable phases is reminiscent of the lobe structure characteristic of the equilibrium model. But despite similar shapes, equilibrium and non-equilibrium lobes are very different in nature. In particular, as shown in Fig. 2, the mean photon density is not constant within the lobes. Moreover, in the bistable region the two phases have very different photon density: one of them, hereafter called low-density phase, has b † b ∼ 10 -2 (left panel of Fig. 2) whereas in the high density phase b † b 1. (right panel of Fig. 2). To gain further insight, let us consider quantum correlations of the considered phases. The on-site second order correlation function g 2 (0) = b † b † bb / b † b 2 is plotted in Fig. 3. It is important to keep in mind that at equilibrium the value of g 2 (0) inside the lobes is equal to 1 -1/n for a pure Mott insulator state where n is the constant integer number of particles on each site, Right panel: the same quantity but for the high-density phase.

For sake of clarity, the maximal value of the colorscale in the high-density phase has been set to 3, but the density is higher than 10 at high J and low U . Notice that the monostable region is the same for both panels. whereas it goes to 1 when J ≫ U in the thermodynamical limit. The left panel of Fig. 3 shows that the light is antibunched inside the lobes. The lowest value, 0.6, is observed for the upper lobe, reminiscent of the n = 2lobe at equilibrium. But the analogy with equilibrium stops here : in the low density phase, for U/∆ω < 3 the emitted light shows strong bunching, that is g 2 (0) ≫ 1. Furthermore, there is no one-to-one correspondence between photon density and second order correlations. For example, in the low-density phase, for U/∆ω = 2 and J/∆ω = 3.25 the density is 0.026 and g 2 (0) = 13. For U/∆ω = 4.5 and J/∆ω = 1.6, the density remains the same but g 2 (0) = 0.64.

Let us now investigate the collective excitations and the dynamical stability of the two phases. This can be done by linearizing the master equation in Fock space for small fluctuations around the steady-state assuming a Gutzwiller factorization of the density matrix:

ρ = i (ρ + δρ i ), ( 6
)
where ρ is the steady-state density matrix whose representation in Fock space can be extracted from Eq.( 5) [START_REF]a steadystate solution of the master equation (2) assuming the Gutzwiller ansatz[END_REF].

-10

0 10 Re(ω) [γ] M Γ Γ X -1 0 1 Im(ω) [γ] Γ X M Γ -10 0 10 Re(ω) [γ] Γ X Γ M -1 0 1 Im(ω) [γ] Γ X Γ M
FIG. 4: Color online. Energy-momentum dispersion of elementary excitations for points A (upper panel) and A' (lower panel), indicated in Fig. 1. Real and imaginary part of the low-energy branches (in units of γ) are plotted vs k. Γ = (0, 0), M = (π/a, π/a), X = (π/a, 0) are special points in the Brillouin zone of the squared photonic lattice. Thick blue lines depict branches with a flat real part over the entire Brillouin zone, while the imaginary part is strongly dispersive with a resonance around specific wavevectors.

-10

0 10 Re(ω) [γ] Γ X Γ M -1 0 1 Im(ω) [γ] Γ Γ M X -10 0 10 Re(ω) [γ] Γ X Γ M -1 0 1 Im(ω) [γ] Γ Γ M X
FIG. 5: Excitation dispersions for points B (upper panel) and B' (lower panel). Same conditions as in Fig. 4 The coefficients of δρ i then obey linear differential equations, which are coupled in real space but decoupled in reciprocal space due to the translational symmetry. For each k-vector in the Brillouin zone, we have:

i∂ t δρ k = L k .δρ k , (7)
where

δρ k n,m = 1 √ N N i=1 e -ik•r i δρ i n,m
and L k is the matrix associated to the linearization. The energy spectrum is given by the eigenvalues of L k and the system is dynamically stable if all eigenvalues have negative imaginary part. Stability studies revealed the onset of modulational instabilities in the low density phase (regions in yellow and red in Fig. 1). Dispersion relations at the edge of the unstable region (points A, A', B and B' on Fig. 1), are plotted on Fig. 4 and5. Remarkably, the real part of the unstable branch is zero for every k inside the Brillouin zone while the imaginary part is strongly dispersive. The existence of this purely imaginary branch can be seen analytically in the low-density regime where there is at most one photon per site so that we can approximate our description by working in a truncated Hilbert space. The vector δρ k has then only four coefficients, (δρ k 00 , δρ k 01 , δρ k 11 , δρ k 10 ) T and L k is given by:

L k =     0 (A * + F * ) + t k ρ * 10 iγ -(A + F ) -t k ρ 10 A + F -∆ω -t k (ρ 00 -ρ 11 ) -i γ 2 -(A + F ) 0 0 -(A * + F * ) -t k ρ * 10 -iγ (A + F ) + t k ρ 10 -(A * + F * ) 0 A * + F * ∆ω + t k (ρ 00 -ρ 11 ) -i γ 2     (8)
where A = -J b is the mean-field parameter, (ρ 00 , ρ 10 , ρ 11 ) the coefficients of the steady-state density matrix, t k = -J/2(cos(k x a) + cos(k y a)) is the term responsible for the dispersion in a square lattice and a is the lattice parameter.

Simple algebra shows that such 4 × 4-matrix has always a purely imaginary eigenvalue. We have checked numerically that in the low density regime such truncated matrix agrees with the results obtained by including the full linearization matrix. The dispersive nature of purely imaginary branches is a consequence of interactions. Indeed, in the low-density phase, when J ≫ U , all eigenvalues have completely flat imaginary parts, while there are only two (anti-conjugate) branches whose dispersive real parts originate from the bare boson dispersion on a square lattice. The k-dependence of L k is enclosed in the coefficient t k which itself is a function of cos(k x a) + cos(k y a). We see on Fig. 4 and Fig. 5 that on the left side of the unstable region, instabilities arise at k = (π/a, π/a). On the other edge however, the unstable k-vectors are smaller and located well inside the Brillouin zone. The region marked in red on Fig. 1 is of particular interest as there is no stable homogenous solution inside. This shows that an inhomogeneous density-wave steady-state occurs in this region of the phase diagram.

In summary, we have explored the mean-field phase diagram of a driven-dissipative Bose-Hubbard model for a wide range of parameters. In the case of spatially homogenous coherent pumping, depending on the values of the on-site repulsion and tunneling coupling, it is possible to have bistable or monostable homogeoneous solutions with peculiar quantum correlation properties. A collective excitation mode with a flat dispersion over the entire Brillouin zone and a dispersive imaginary part can occur, leading to tunneling-induced instabilities at specific wavevectors and thus a breaking of the translational invariance. Our results shows that driven-dissipative arrays of cavities can lead to very rich manybody physics, which is very different from its equilibrium counterpart. C.C. is member of Institut Universitaire de France.

1 FIG. 1 :

 11 FIG.1: Color online. Top: Sketch of a square photonic lattice made of nonlinear cavities coupled by tunneling. The system is pumped coherently by an homogeneous laser field at normal incidence. Bottom: Number of mean-field solutions and their stability plotted as a function of J/∆ω and U/∆ω, for F/∆ω = 0.4 ; γ/∆ω = 0.2 and ∆ω > 0. Light blue (top-left and bottom-right part labeled with a '1'): monostable region, only one solution to Eq.(4). Dark blue part (label '2'): bistable region, two solutions to Eq.(4). Yellow part (central region labeled with '1') has only one stable phase out of two existing solutions. Red part (label '0'): only one solution, which is unstable.

3 FIG. 2 :

 32 FIG.2: Photon occupation number as a function of J/∆ω and U/∆ω for F/∆ω = 0.4 and γ/∆ω = 0.2. Left panel: results for the boson occupation number in the low-density phase. Right panel: the same quantity but for the high-density phase. For sake of clarity, the maximal value of the colorscale in the high-density phase has been set to 3, but the density is higher than 10 at high J and low U . Notice that the monostable region is the same for both panels.

3 FIG. 3 :

 33 FIG.3: Second-order correlation function g 2 (0) as a function of J/∆ω and U/∆ω. Same parameters as in Fig.2. Left panel for the low-density phase. Right panel for the highdensity phase.