# Signatures of Self-Trapping in the Driven-Dissipative Bose-Hubbard Dimer

SSignatures of Self-Trapping in the Driven-Dissipative Bose-Hubbard Dimer

Matteo Seclì, ∗ Massimo Capone,

1, 2 and Marco Schirò † International School for Advanced Studies (SISSA), Via Bonomea 265, I-34136 Trieste, Italy CNR-IOM Democritos, Via Bonomea 265, I-34136 Trieste, Italy JEIP, USR 3573 CNRS, Collége de France, PSL Research University,11 Place Marcelin Berthelot, 75321 Paris Cedex 05, France (Dated: February 08, 2021)We investigate signatures of a self-trapping transition in the driven-dissipative Bose Hubbarddimer, in presence of incoherent pump and single-particle losses. For fully symmetric couplings thestationary state density matrix is independent of any Hamiltonian parameter, and cannot thereforecapture the competition between hopping-induced delocalization and the interaction-dominated self-trapping regime. We focus instead on the exact quantum dynamics of the particle imbalance afterthe system is prepared in a variety of initial states, and on the frequency-resolved spectral propertiesof the steady state, as encoded in the single-particle Green’s functions. We ﬁnd clear signaturesof a localization-delocalization crossover as a function of hopping to interaction ratio. We furthershow that a ﬁnite a pump-loss asymmetry restores a delocalization crossover in the steady-stateimbalance and leads to a ﬁnite intra-dimer dissipation.

I. INTRODUCTION

Recent years have seen an increase of interest in openMarkovian quantum systems, which describe a numberof experimental platforms for quantum information pro-cessing and quantum simulation, both in the realm ofatomic physics and quantum optics as well as in the solidstate framework. Among these we can mention for ex-ample cavity QED experiments [1] and their analoguewith superconducting circuits [2]. Here the basic degreesof freedom, photons and qubits, are inevitably exposedto dissipative processes such as losses and decoherenceinduced by the environment. The quantum dynamicsof Markovian systems is described theoretically withinthe framework of a Lindblad master equation which en-codes the competition between coherent (Hamiltonian)evolution and dissipative processes described by a set ofjump operators [3]. Out of this competition one can ex-pect non-trivial stationary states and dynamical behav-ior to emerge, leading to novel dissipative phase transi-tions [4, 5], both in small systems made by few quantumnon-linear oscillators [6, 7] as well as in larger arrays [8–14].An intriguing question which has recently attractedlarge interest is to understand what kind of dynamicalphenomena can arise in these Markovian quantum sys-tems and their relationship with analogous phenomenain the ﬁeld of classical non-linear dynamical systems inpresence of non-linearities, noise and dissipation [15].A prototype example in this context is provided by thedriven-dissipative Bose-Hubbard dimer (BHD), whichcan be seen as a toy model of strongly correlated openMarkovian quantum systems since it encodes the basic ∗ [email protected] † On Leave from: Institut de Physique Théorique, Université ParisSaclay, CNRS, CEA, F-91191 Gif-sur-Yvette, France competition between local dissipative processes, interac-tions and non-local coherent hopping processes.Besides its paradigmatic relevance, the driven-dissipative BHD has also been realized experimentallyin a variety of quantum light-matter platforms, includ-ing superconducting circuits [16, 17] and semiconductormicrocavities [18–20] and photonic crystals [21, 22].In the closed isolated case, corresponding to a purelyconservative Hamiltonian evolution, the BHD has beenextensively studied, in particular its self-trapping, orlocalization-delocalization [23–29]. Here, an initial im-balance of particles between the two sites of the dimeris either rapidly redistributed by hopping processes lead-ing to an homogeneous conﬁguration or conserved indef-initely, leading to a self-trapped state below a criticalratio between hopping and interaction. This transitioncorresponds to a spontaneous breaking of the reﬂectionsymmetry between the two sites of the dimer. Open-Markovian extensions of the BHD have been mostly fo-cused on the coherently driven case [17, 30–33] or, in thecase of the related Jaynes-Cummings Dimer model [34],the purely dissipative case in absence of any externalpumping.In this work we theoretically study the driven-dissipative BHD in presence of single-particle losses andincoherent single-particle drive. This case is somewhatpeculiar, since it is known that for a perfectly symmetricmodel the stationary state of the problem is completelyindependent of Hamiltonian parameters and only set bythe ratio between pump and losses [35], so it cannot con-tain any signature of a putative delocalization transition.In order to explore the competition between hopping andinteractions in a dissipative setting one has therefore togo beyond the analysis of steady-state observables andfocus instead on response functions, or to consider anasymmetry between the two sites of the dimer.In particular we prepare the system in diﬀerent ini-tial states and follow the exact quantum dynamics of themodel, characterizing also the properties of the station- a r X i v : . [ qu a n t - ph ] F e b ary state reached at long times. Furthermore we focuson the spectral properties of the BHD as encoded inthe Green’s functions which for open-Markovian quan-tum system, much like their closed system counterpart,contain rich insights on the structure of the single-particleexcitations around the stationary state.The paper is organized as follows. In Sec. II we in-troduce the BHD model and brieﬂy review some of itsproperties, while in Sec. III we present details on its nu-merical solution. In Sec. IV we review the known resultsabout the semiclassical limit and the self-trapping transi-tion in the isolated and dissipative cases. Our results forthe quantum dynamics in the symmetric pumping regimeare discussed in Sec. V, while those for ﬁnite pump/lossasymmetry in Sec. VI. In Sec. VII we present results forthe Green’s functions of the BHD, while Sec. VIII is de-voted to conclusions. II. THE MODEL

We start by considering the Hamiltonian of a Bose-Hubbard dimer (BHD). The model is a paradigmatic in-teracting lattice model which can be realized in a num-ber of platforms. Our implementation including pumpingand losses is naturally realized using optical cavities. Forthis reason in the following we will refer to the two latticesites as cavities and to the bosonic degrees of freedom in-volved in the physics as photons. The Hamiltonian reads ˆ H = ω (ˆ n L + ˆ n R ) + U (cid:16) ˆ n L ˆ n L + ˆ n R ˆ n R (cid:17) + J (cid:16) ˆ a † L ˆ a R + ˆ a † R ˆ a L (cid:17) , (1)where ˆ n L = ˆ a † L ˆ a L and ˆ n R = ˆ a † R ˆ a R are the number oper-ators of the left and the right cavities, respectively. Thetwo cavities have the same resonant frequency ω andKerr non-linearity U , and photons can hop between thecavities at a rate J .We can add a simple mechanism for incoherent driv-ing and dissipation at the master-equation level, by us-ing single-particle pump and loss operators. In practice,we describe the driven-dissipative dimer by a reduceddensity matrix ˆ ρ that evolves according to the Lindbladmaster equation ˙ˆ ρ = ˆˆ L ˆ ρ = ˆˆ L H ˆ ρ + ˆˆ L D ˆ ρ (2)where ˆˆ L H ρ = − i (cid:104) ˆ H, ˆ ρ (cid:105) (3)is the Hermitian part of the evolution, while the dissipa-tive piece reads as ˆˆ L D ˆ ρ = 2 (cid:88) i = L,R (cid:40) Γ i (cid:18) ˆ a i ˆ ρ ˆ a † i − (cid:110) ˆ a † i ˆ a i , ˆ ρ (cid:111)(cid:19) + P i (cid:18) ˆ a † i ˆ ρ ˆ a i − (cid:110) ˆ a i ˆ a † i , ˆ ρ (cid:111)(cid:19) (cid:41) . (4) In this form, Γ L/R are interpreted as loss rates while P L/R as pumping rates. It is convenient to parametrize themas Γ i = Γ ± ∆Γ , ∆Γ = Γ L − Γ R (5) P i = P ± ∆ P, ∆ P = P L − P R (6)to distinguish the case in which pump/loss rates are sym-metric in the dimer, ∆Γ = ∆ P = 0 or asymmetric dueto an imbalance of pump and/or losses. In fact it isknown [35] that for a Bose-Hubbard lattice with uniformparameters and identical single-particle pump and lossrates, i.e. ∆Γ = ∆ P = 0 the structure of the stationarystate density matrix is particularly simple and reads ˆ ρ ss = (cid:88) N π N | N (cid:105) (cid:104) N | where | N (cid:105) is a Fock state with N bosons and π N ∼ ( P/ Γ) N up to a normalization factor. We note in theabove expression that ˆ ρ ss is independent of any Hamil-tonian parameter and only set by pump/loss ratio. Thisimplies in particular that the stationary state occupancy n α = Tr ( ρ ss ˆ n α ) is equal in the two cavities and given by n L = n R = P Γ − P (7)which coincides with the value of an uncoupled Kerr res-onator. Given these results, it is clear that any non-trivial dependence from J/U has to be looked for in prop-erties other than the stationary-state observables, as wewill discuss in Sec. V and VII A. The above result is how-ever no longer true in presence of a ﬁnite asymmetry inthe dissipative couplings, leading to ∆Γ , ∆ P (cid:54) = 0 , as wewill see more in detail in Sec. VI and VII B. III. METHODS

The vectorized version of equation (2) is solved by ex-act diagonalization, yielding a bi-normalized set of leftand right eigenvectors ( (cid:104) l α | and | r α (cid:105) , respectively) thatsatisfy (cid:104) l α | ˆ L = L α (cid:104) l α | and ˆ L | r α (cid:105) = L α | r α (cid:105) (8)where ˆ L is the matrix representation of the superoperator ˆˆ L . The cokernel and the kernel of ˆ L are, respectively,the left vacuum (cid:104) I | and the steady-state density matrix | ρ ss (cid:105) .The diagonalization problem can actually be simpliﬁedby realizing that both the Hamiltonian and the dissipatorconserve the total number of photons, yielding a global The left and right eigenvectors corresponding to the special eigen-value L = 0 . gauge symmetry expressed by an operator functional ˆˆ K that commutes with ˆˆ L . By exploiting this symmetry thematrix ˆ L can then be written in a block-diagonal form,where each block is labeled by the eigenvalues of ˆˆ K . A. Time Dynamics

Having solved the eigenproblem, we can then expand[36] | ρ ( t ) (cid:105) = e L t | ρ (0) (cid:105) = (cid:88) α ρ α ( t ) | r α (cid:105) , (9)where ρ α ( t ) (cid:43) e L α t (cid:104) l α | ρ (0) (cid:105) = e L α t ρ α (0) . (10)We note that the form of the Lindblad equation ensures Re L α ≤ ∀ α , which prevents the dynamics from un-bounded growth with time. Again, if we can exploit theglobal gauge symmetry, then it is suﬃcient to diagonalizejust the largest diagonal block of the Lindbladian. Theknowledge of the time-evolution of the density matrix canthen be used to calculate the time-evolution of other ob-servables, for example the occupations of the two cavities( i = { L, R } ): n i ( t ) = Tr (cid:16) ˆ n i ρ ( t ) (cid:17) = (cid:88) α (cid:104) I | ˆ n i | r α (cid:105) ρ α ( t ) . (11) B. Källén-Lehmann Spectral Representation ofGreen’s Functions

Albeit not necessary if one only wants to calculate thesteady-state density matrix | ρ ss (cid:105) , the full knowledge ofthe spectrum can be used to explore the Green’s functionsof the system. In fact, one can obtain frequency-domainexpressions for the retarded and the Keldysh componentsof the steady-state Green’s function deﬁned respectivelyas G RAB ( t ) = − iθ ( t ) (cid:104) [ A ( t ) , B (0)] (cid:105) (12) G KAB ( t ) = − i (cid:104){ A ( t ) , B (0) }(cid:105) (13)where the average is taken over the stationary stateand the operator A is evolved with the Lindbladianof the system. Upon inserting a complete set ofleft and right eigenvectors of the Lindbladian and go-ing to the frequency domain by deﬁning G R/KAB ( ω ) (cid:43) (cid:82) dt e iωt G R/KAB ( ω ) , we obtain a spectral representation ofthose functions: G RAB ( ω ) = (cid:88) α (cid:104) I | A | r α (cid:105) (cid:104) l α | B | ρ ss (cid:105) ω − i L α − (cid:32)(cid:88) α (cid:10) I (cid:12)(cid:12) A † (cid:12)(cid:12) r α (cid:11) (cid:10) l α (cid:12)(cid:12) B † (cid:12)(cid:12) ρ ss (cid:11) ω + i L α (cid:33) ∗ (14) G KAB ( ω ) = (cid:88) α (cid:104) I | A | r α (cid:105) (cid:104) l α | B | ρ ss (cid:105) ω − i L α − (cid:88) α (cid:104) I | B | r α (cid:105) (cid:104) l α | A | ρ ss (cid:105) ω + i L α + (cid:32)(cid:88) α (cid:10) I (cid:12)(cid:12) A † (cid:12)(cid:12) r α (cid:11) (cid:10) l α (cid:12)(cid:12) B † (cid:12)(cid:12) ρ ss (cid:11) ω + i L α (cid:33) ∗ − (cid:32)(cid:88) α (cid:10) I (cid:12)(cid:12) B † (cid:12)(cid:12) r α (cid:11) (cid:10) l α (cid:12)(cid:12) A † (cid:12)(cid:12) ρ ss (cid:11) ω − i L α (cid:33) ∗ (15)We see that the Green’s functions of an open Markovianquantum system can be generically written as sum ofsimple poles at complex frequencies given by the eigen-values of the Lindbladian and with weights, in generalcomplex, given by the transition matrix elements betweenthe stationary state and some excited state of the system[36, 37].From the practical point of view, if one focuses onthe single-particle Green’s functions, the calculation canbe further simpliﬁed via the block-diagonal structure ofthe Lindbladian outlined above. In fact, since the cal-culation of the single-particle Green’s functions involvesstates that diﬀer at most by one particle from the sta-tionary state, it turns out that the full knowledge of thespectrum is not necessary; it is suﬃcient to diagonalizejust the 3 largest blocks of the diagonal-block structure.Assuming that the diagonlization scales as the cube ofthe matrix linear dimension, this yielded a theoretical speedup of the diagonalization with the 20-bosonscutoﬀ we have used in both cavities, as well as a . reduction of the memory required to store the results. IV. REVIEW OF SEMICLASSICAL DYNAMICSAND SELF-TRAPPING TRANSITION

In order to have a reference point for the analysis ofour results we can start by recalling the predictions of asemiclassical treatment of the quantum dynamics for theBHD [23, 38]. This is obtained by writing the exact equa-tions of motion for the cavity ﬁeld operators ˆ a L/R andby closing them by taking ˆ a L/R ≡ α L/R , where α L/R are c -numbers. It is important to remark that this approach,which assumes a coherent state of bosons, works for largephotons number, while in the quantum treatment we aretypically interested in a few-photons treatment. The re-sulting equations of motion read ˙ α L = − i (cid:104) ( ω − U ) + 2 U | α L | (cid:105) α L − iJα R − Γ eﬀ L α L ˙ α R = − i (cid:104) ( ω − U ) + 2 U | α R | (cid:105) α R − iJα L − Γ eﬀ R α R where Γ eﬀ L/R = Γ

L/R − P L/R are the eﬀective loss rates,which for single-particle losses must always be positive.As discussed in more detail in Appendix A, it’s possibleto write semiclassical equations for the total number of

Figure 1. Evolution of the occupation imbalance for N = 3 , Z = 1 , U = 0 . and ∆ ω = 0 . Diﬀerent colors correspondto diﬀerent values of J/U around the critical value, predictedvia Eq. (16). The lines in the top panel are for a closedsystem, while the ones of the same color in the bottom panelare for an open system with Γ eﬀ L = Γ eﬀ R = 4 × − ; they areobtained respectively by numerically solving the full systemof equations discussed in Section IV (see also Eq. (A2)-(A3)in Appendix A). photons N = n L + n R and for the occupation imbalancebetween the two cavities Z = n L − n R , with n L/R = | α L/R | .In the closed-system case, corresponding to Γ eﬀ L/R = 0 ,number and energy conservation yield simpliﬁed analyt-ical results for the imbalance Z , predicting a transitionfrom a regime in which Z oscillates above the initial con-dition Z to a regime in which it oscillates around (solidlines in Fig. 1) as one increases the value of J/U abovethe critical coupling (cid:18) JU (cid:19) c = N (cid:32) (cid:112) − ( Z /N ) + 12 (cid:33) . (16)which depends on the initial total number of photons N and imbalance Z . This phase transition can be seenas a divergence of the oscillation period (Fig. 11) or asa sharp decay to zero of the time-averaged imbalance(Fig. 2, bottom panel) (see Appendix A).The open system case is not analytically solvable, butthe numerical solution of the equations for the total num-ber of photons and for the cavity occupation imbalanceshows that the closed-system picture is preserved for lowenough values of the loss coeﬃcients, with the diﬀerencethat even oscillations around a value that is diﬀerentfrom zero at initial times will eventually transition atlong enough times to an oscillation regime around zeroduring the dynamical evolution (Fig. 1, bottom panel).We can deﬁne the time at which this dynamical tran-sition happens to be some t cross for which the imbalance Figure 2. (Top) Time t cross at which Z ( t ) crosses the value Z = 0 , as a function of J/U . (Bottom) Semiclassical time-averaged imbalance. While the closed system has an analyti-cal expression, the open case requires to solve the full dynam-ics (A2) and to choose an upper time limit in the integration;in this plot, we integrate up to t = 200 . In both panels, thecritical value of J/U predicted via (16) (for the closed system)and as a numerical estimate (for the open system) is shownas a vertical dotted line. Z ( t ) crosses the value Z = 0 for the ﬁrst time. If we plotthis time as a function of J/U we expect that for theclosed system this time is divergent for values of

J/U be-low the critical value; for the open system, however, thistime assumes ﬁnite values even below the critical pointand the critical point itself is at a slightly lower valuethan its closed-system counterpart ( ( J/U ) c = 2 . vs. ( J/U ) c = 2 . ).Albeit holding in the limit of large photon number only,these semiclassical results provide a useful hint for thequantities to look at in the quantum case, as well as apoint of comparison that highlights the intrinsic diﬀer-ences between the two types of analyses. V. RESULTS: DISSIPATIVE QUANTUMDYNAMICS

We now move on to discuss the full dissipative quan-tum dynamics of the BHD introduced in Sec. II. We fo-cus in particular on the occupation imbalance Z ( t ) = n L ( t ) − n R ( t ) between the two cavities, which in thesemiclassical limit shows a clear change of behavior asa function of the parameters.In the following we set ω = 1 , U = 0 . and consider asituation of symmetric pump and loss rates, ∆ P = ∆Γ =0 , so that by construction the imbalance is zero at longtimes. We set the eﬀective losses Γ eﬀ L/R = Γ

L/R − P L/R =1 × − and the pump P L/R = 2 × − , such that the Figure 3. (Inset) Imbalance Z ( t ) = n L ( t ) − n R ( t ) for dif-ferent values of J/U = 0 . [ U = 0 . ], starting from a state | , (cid:105) at t = 0 ( Z = 2 ). The cavities have a base frequency ω = 1 . . The eﬀective loss is Γ eﬀ L = Γ eﬀ R = 1 × − and thepumping rate realizes a steady-state occupation equal to inboth cavities, so that Z ss = 0 by construction. The semiclas-sical, non-dissipative critical value of J/U for this particularconﬁguration is ( J/U ) c ≈ . . The time-averaged occupa-tion imbalance computed over the time interval [0 , isshown in the main panel. identical occupation in the two cavities is n L = n R = 2 (see Eq. (7), independently on J/U .We start discussing the imbalance dynamics as a func-tion from

J/U , at a ﬁxed initial condition which we taketo be a Fock state | , (cid:105) , corresponding to an initial im-balance Z = 2 and an initial number of photons N = 4 .At the semiclassical level, see Eq. (16), this would corre-spond to a critical coupling ( J/U ) c = 3 . for the self-trapping transition.In the inset of Fig. 3 we plot the time-dependent im-balance Z ( t ) for diﬀerent values of J/U . We ﬁnd a clearcrossover as the hopping is increased, from a pure ex-ponential decay to zero at small

J/U = 0 . , to an un-derdamped decay with fast oscillations superimposed at J/U = 0 . which evolves further into strongly anhar-monic oscillations at large values of the hopping, whosefrequency grows with J/U . We can interpret this be-havior as a signature of the self-trapping transition inthe dissipative quantum dynamics. In the small hoppingregime each site of the dimer evolves almost indepen-dently and the imbalance goes to zero, while for largervalues of the hopping there is a substantial transfer ofphotons across the dimer, resulting in coherent Rabi-likeoscillations, before the imbalance reaches the stationarystate.The

J/U dependence can also be studied from thepoint of view of the time-averaged occupation imbal-ance (cid:104) Z (cid:105) T . In contrast to the semiclassical case (Fig. 2), Figure 4. Evolution of the imbalance Z ( t ) for the samesettings of Fig. 3, shown at longer times and at log scale. Atlong times, log ( Z ( t ) /Z ) ﬁts a straight line in log scale; theinset shows the corresponding decay rate as a function of J/U . where one expects a sharp transition between (cid:104) Z (cid:105) T (cid:54) = 0 and (cid:104) Z (cid:105) T = 0 , in the quantum case we have a smoothcrossover between the two regimes. The average imbal-ance drops quickly with J/U due to the development ofdamped Rabi oscillations, reaching a minimum around

J/U (cid:39) . . Quite interestingly, though, we ﬁnd theappearance of a region in which the imbalance actu-ally increases as a function of J/U before completelydropping to at higher values of J/U . We note that,with respect to the semiclassical case, the localized (self-trapped) phase with (cid:104) Z (cid:105) T (cid:54) = 0 is strongly suppressed andthat already for J/U (cid:39) . the average imbalance iszero. This is consistent with the expectation that quan-tum ﬂuctuations, included in the exact solution and notproperly treated in the semiclassical approach, tend toreduce the broken symmetry phase.We now discuss the dynamics on longer time scales,where we expect the small dissipative couplings to domi-nate over the Hamiltonian parameters. To this extent inFig. 4 we plot the time-dependent imbalance over a broadrange of time scales and for diﬀerent values of J/U . Wesee a clear separation of dynamical regimes, from a short-time one - strongly dependent on

J/U , as we discussedabove - to a longer-time one where the imbalance expo-nentially decays to zero. While naively one could haveexpected the decay rate to be set only by the dissipa-tive couplings we see in the inset of Fig. 4 that instead itshows a monotonic increase with

J/U .Finally, we consider the dependence of the time-dependent imbalance Z ( t ) from the initial condition. To In the open case, the extent of the jump discontinuity in ∂ J/U (cid:104) Z (cid:105) T depends on the upper limit of the integration time. Figure 5. Imbalance Z ( t ) = n L ( t ) − n R ( t ) for ( J/U ) =0 . J/U ) c [ U = 0 . ], starting from diﬀerent | n L , n R (cid:105) num-ber states at t = 0 . The cavities have a base frequency ω = 1 . . The eﬀective loss is Γ eﬀ L = Γ eﬀ R = 1 × − andthe pumping rate realizes a steady-state occupation equal to in both cavities, so that Z ss = 0 by construction. The quan-tity ( J/U ) c refers to the semiclassical non-dissipative value in(16). this extent we ﬁx as initial density matrix a pure Fockstate ρ = | n L , n R (cid:105) (cid:104) n L , n R | , corresponding to an ini-tial imbalance Z = n L − n R and initial photon num-ber N = n L + n R , and change the values of n L , n R .At the semiclassical level, as we see in Eq. (16), thereis a critical value of J/U for any N , Z . In order tohighlight the diﬀerence between the exact quantum dy-namics and the semiclassical evolution we ﬁx the value ofthe hopping to interaction ratio J/U to be always below ( J/U ) c ( N , Z ) , such that at the semiclassical level thesystem should be localized (self-trapped) at short timesfor all the chosen initial conditions (see Eq. (16)) anddelocalized at longer times (see Fig. 1).We plot in Fig. 5 the quantum dynamics of the imbal-ance for diﬀerent initial conditions. We see that, quiteat the opposite of what expected from the semiclassicalanalysis, the evolution of Z ( t ) has a strong dependenceon the initial state in which the system is prepared. Inparticular we ﬁnd both regimes of slow decay to zeroof the imbalance (see for example the initial conditionscorresponding to | , (cid:105) or | , (cid:105) ), indicating localized/self-trapped behavior, as well as regimes of coherent Rabi-likeoscillations of the imbalance (see for example the initialconditions corresponding to | , (cid:105) or | , (cid:105) ) that we can in-terpret as signatures of delocalization. This is consistentwith the observation made earlier (see Fig. 3) that quan-tum ﬂuctuations renormalize the critical coupling and fa-vor the delocalized regime. We conclude therefore that, Figure 6. (Top) Steady-state cavity occupations as a func-tion of

J/U for a dimer with loss coeﬃcients (Γ L , Γ R ) =(6 × − , × − ) and pump coeﬃcients ( P L , P R ) = (4 × − , × − ) described in the main text, at diﬀerent valuesof U . The cavities have a base frequency ω = 1 . . The topcurves are the occupations of the left cavity, while the bot-tom ones are the occupations of the right cavity. (Bottom)Steady-state imbalance Z = n L − n R corresponding to theoccupations in the top panel. as in the semiclassical case, the self-trapping crossovercan be accessed by changing the initial condition, how-ever we do not explore here the precise dependence of ( J/U ) c from the initial state and whether it can be en-coded in a simple expression depending only on N and Z as in Eq. (16). VI. RESULTS: QUANTUM STEADY STATEFOR FINITE PUMP/LOSS ASYMMETRY

In the previous section we have considered the case ofa BHD with symmetric pump and loss rates, resultingin a trivial stationary state with zero imbalance for anyvalue of

J/U , but with a rich nonequilibrium dynamics.As we discussed in Sec. II, in presence of a ﬁnitepump/loss asymmetry among the two cavities the sta-tionary state becomes more interesting. We can thereforelook for signatures of a delocalization crossover, analo-gous to what we have shown in Fig. 3, directly in observ-ables such as the steady-state occupation or imbalance.As an example, we consider two cavities with loss co-eﬃcients (Γ L , Γ R ) = (6 × − , × − ) and pump co-eﬃcients ( P L , P R ) = (4 × − , × − ) , that thus re-alize steady-state occupations ( n L , n R ) = (2 , in theuncoupled limit J = 0 (see Eq. (7)). In Fig. 6 we plotthe dependence of the two cavity occupations (top panel)and imbalance (bottom panel) from the hopping to in-teraction ratio J/U , for diﬀerent values of U (keeping ω = 1 as unit). We see in the top panel that as J/U isincreased the two occupations both converge towards acommon value, which is essentially independent from U .The long-time limit of the occupations can be obtainedanalytically by considering the limit U = 0 and resultsin a weighted average of the two uncoupled occupations(see Eq. (B12)).As a consequence of the two occupations becomingequal at large J/U we see in the bottom panel that thesteady-state imbalance between the two cavities reducesand approaches zero for large enough

J/U , a signatureof delocalization. We note that increasing U pushes thecrossover J/U scale for delocalization to lower values andwe expect for U = 0 . to obtain a behavior comparablewith what obtained from the dynamics (see Fig. 3). VII. RESULTS: GREEN’S FUNCTIONS

A way to get some insights on the system even whenthe steady-state observables do not depend neither on J nor on U , as in the case of symmetric pump and losses, isto look instead at the single-particle Green’s functions.Either by seeing them as the resolvent of the Lindbla-dian or as response functions that link diﬀerent statesand thus participate in the calculation of transport quan-tities like the optical transmission, the Green’s functionsare sensitive to the details of the Lindbladian spectrum,and not only to the zero mode (stationary state), as itappears clearly from the Källén-Lehmann representationdiscussed in Sec. III B.In this section we present our results for the Green’sfunction of the BHD, that we obtained from the exact di-agonalization of the Lindbladian as discussed in Sec. III.Speciﬁcally we consider the single-particle Green’s func-tions, obtained from Eq. (15) with the choice A = a i and B = a † j with i, j = L/R , and in particular the spec-tral function A ij ( ω ) and the cavity correlation function C ij ( ω ) , deﬁned as A ij ( ω ) (cid:43) − π Im G Rij ( ω ) , C ij ( ω ) (cid:43) − πi G Kij ( ω ) (17)with i, j = L/R . The diagonal components (for i = j )contain information on the local (on-site) spectrum andoccupations of the bosonic mode and satisfy the sum rules (cid:90) + ∞−∞ dω A i ( ω ) = 1 (18) (cid:90) + ∞−∞ dω C i ( ω ) = 2 n i + 1 (19)where n i is the stationary state occupation. The oﬀ-diagonal components contain instead information on thedelocalized modes across the dimer. In particular thecorrelation function C LR ( ω ) has real and imaginary parts which satisfy the sum-rules J (cid:90) + ∞−∞ dω Re C LR ( ω ) = (cid:104) ˆ T (cid:105) (20) J (cid:90) + ∞−∞ dω Im C LR ( ω ) = (cid:104) ˆ I (cid:105) (21)where (cid:104) ˆ T (cid:105) = J (cid:104) ˆ a † L ˆ a R + ˆ a † R ˆ a L (cid:105) is the average ki-netic energy in the stationary state while (cid:104) ˆ I (cid:105) = − iJ (cid:104) ˆ a † R ˆ a L − ˆ a † L ˆ a R (cid:105) is the average current ﬂowing from L to R (see Appendix C). We now presents our resultsfor these Green’s functions, starting from the pump/losssymmetric case and then discussing the role of a ﬁnitepump/loss asymmetry. A. Symmetric Pump and Losses

We start considering the case of symmetric pump andloss rates, ∆Γ = ∆ P = 0 . As a result the system iscompletely symmetric upon reﬂection ( L ↔ R ) and assuch the diagonal spectral functions in Eq. (17) do notdepend on the index i = L/R . As an example, in the toppanel of Fig. 7 we plot the spectral function of the leftcavity for diﬀerent values of

J/U .At low

J/U the spectral function resembles much theone of a single driven-dissipative Kerr resonator, witha characteristic sequence of peaks located at frequenciesgiven by the energy diﬀerence between states with n + 1 and n photons, ∆ n = E n +1 − E n = ω + U + 2 U n , where E n = ω n + U n is the energy of the Kerr resonatorwith n photons (see the Hamiltonian in Eq. (1)). Thesepeaks, which start at ω + U and are equally spaced by U , would be inﬁnitely sharp in the closed system whileare broadened by the dissipative processes by an amountroughly given by Γ eﬀ L (it would perfectly match this valuein the non-interacting, decoupled case J = U = 0 , seeAppendix B 1.)As J/U is increased we see that the ﬁrst eﬀect is thecreation of sub-peaks within each resonance, particularlyin the low frequency ones, with the center of mass of each band remaining roughly located at the isolated Kerr ex-citation energies. Upon increasing further

J/U we seehow diﬀerent bands start to merge in a continuum andfor

J/U = 0 . a new features arises, namely a ﬁnitespectral weight appears below the resonator frequency ω = 1 , which becomes a sharp peak for large values of J/U (e.g.

J/U = 1 . ). This peak corresponds to a delo-calized photonic excitation as one can realize by lookingat the spectral function in the opposite limit of U = 0 (see Appendix B 1), which has two poles at frequenciesroughly ω ± (cid:39) ω ± J since in this regime the dissipativecouplings are very small.It is interesting to connect these spectral features to thebehavior of the time-dependent and of the time-averagedimbalance shown in Fig. 3 for similar values of J/U . Forsmall values of the hopping the imbalance is diﬀerent

Figure 7. (Top) Spectral function A L ( ω ) for diﬀerent val-ues of J/U [ U = 0 . ]. The cavities have a base frequency ω = 1 . (vertical dotted line), while the eﬀective loss is Γ eﬀ L = Γ eﬀ R = 1 × − and the pumping rate realizes a steady-state occupation equal to in both cavities. The circledpeaks mark the bonding/anti-bonding states resulting fromthe splitting of the ﬁrst excited state at ω + U for decoupledcavities. (Botttom) Real part of the oﬀ-diagonal cavity cor-relation function C LR ( ω ) . This function is negative (positive)for the bonding (anti-bonding) states marked by circles anddiscussed in the top panel.. from zero at short and intermediate times, i.e. photonsremain localized in one of the two cavities and the spec-tral function resembles the one of an isolated Kerr res-onator. Upon increasing J/U photons start to hop coher- ently within the dimer: the imbalance shows short-timeRabi oscillations with a period controlled by

J/U and itstime-average vanishes, while spectrally this translates inthe emergence of two peaks above and below the bareresonator frequency.In the bottom panel of Fig. 7 we plot the real-part ofthe oﬀ-diagonal correlation function, for diﬀerent valuesof

J/U and ∆Γ = ∆ P = 0 . We note that quite in-terestingly the imaginary part of this Green’s functionvanishes in this regime, a point onto which we will comeback in the next section. At small values of the hoppingthe real-part C LR ( ω ) is essentially zero, the cavities arealmost decoupled, except at frequencies corresponding tothe eigenmodes of the (interacting) single cavity (see toppanel at the same value of J/U ), where an anti-resonancelike contribution emerges. Upon increasing

J/U , as wediscussed for the spectral function, further peaks appearwhich start merging and shifting towards lower frequen-cies. We note that the structure of the peaks evolve as

Figure 8. Spectral functions A L ( ω ) (top) and A R ( ω ) (bot-tom) for diﬀerent values of J/U [ U = 0 . ]. The cavities havea base frequency ω = 1 . (vertical dotted line), while theeﬀective losses are Γ eﬀ L = 2 . × − and Γ eﬀ R = 1 × − and the pumping rates realizes uncoupled steady-state occu-pations equal to ∼ . in the left cavity and in the rightcavity. Figure 9. Cavity oﬀ-diagonal correlation function C LR ( ω ) fordiﬀerent values of J/U [ U = 0 . ]. The cavities have a basefrequency ω = 1 . (vertical dotted line), while the eﬀectivelosses are Γ eﬀ L = 2 . × − and Γ eﬀ R = 1 × − and the pump-ing rates realizes uncoupled steady-state occupations equal to ∼ . in the left cavity and in the right cavity. well: at small J/U they are almost perfectly asymmetricin frequency (leading to a vanishing integral, see Eq. (21))while upon increasing

J/U , when the system becomesmore delocalized, this asymmetry disappears. Further-more, also the strength of the peaks increases with

J/U (note the diﬀerent scale in the panels) in a way that ap-pears opposite to the peaks in the spectral function inthe top panel. This is again consistent with the ideathat upon entering in the delocalized regime the weightis transferred from the localized (on-site) modes to thedelocalized (oﬀ-diagonal ones).

B. Asymmetric Pump and Losses

We now move to discuss the case of asymmetric pumpand losses, ∆ P, ∆Γ (cid:54) = 0 , resulting as we know in a nontrivial stationary state density matrix (and ﬁnite imbal-ance, see Sec. VI). A natural question is whether this dif-ferent nonequilibrium protocol results in a qualitativelydiﬀerent behavior of the Green’s functions.We start from the spectral functions, that we plot inFig. 8 for a ﬁxed pump/loss asymmetry and diﬀerent val-ues of J/U . To highlight the comparison between the twocavities we plot the left and right spectral functions on acommon frequency scale. While we see a similar structureof peaks evolving with

J/U , as compared to the symmet-ric case of Fig. 7, we also note an interesting dependencefrom the pump/loss asymmetry and the hopping. In par-ticular, for small

J/U the right cavity spectral function(bottom panels) has slightly stronger peaks at low fre-quency than the left cavity one, reﬂecting the asymme-try in the pump/loss rates. As the hopping is increasedand the excitations are delocalized in the dimer we seethat this asymmetry in the left/right spectral functionsdecreases and for

J/U = 1 . the two spectra are essen-tially the same and very close in shape to the symmetricone for the same value of J/U (See Fig. 7).Then we consider the oﬀ-diagonal cavity correlationfunction, see Fig. 9, that we study as a function of

J/U .In the top panel we plot the real part, Re C LR ( ω ) , whichshows a qualitative behavior very similar to the symmet-ric case shown in Fig. 7, with anti-Lorentzian peaks whichbroaden and merge into a continuum at large J/U indi-cating the increase in kinetic energy. On the other hand,an interesting diﬀerence appears in the imaginary part ofthe oﬀ-diagonal cavity correlation function, Im C LR ( ω ) ,which is now diﬀerent from zero and shows a non-trivialdependence from J/U , with narrow peaks which broadenand merge into a continuum as

J/U is increased.We can understand the origin of a ﬁnite imaginary partof the oﬀ-diagonal cavity correlation function by usingthe sum rule that relates the integral of Im C LR ( ω ) to theaverage current ﬂowing from L to R (see Eq. (21) andAppendix C). In the stationary state the average currentis completely determined by the eﬀective pump/loss rates Γ eﬀ L/R = Γ

L/R − P L/R and the stationary occupation n L/R through the relation (cid:104) ˆ I (cid:105) = ∆ P − n L Γ eﬀ L + n R Γ eﬀ R , (22)where ∆ P is the pump asymmetry. We see that theright-hand side of this equation exactly vanishes in thesymmetric case ∆ P = 0 , Γ eﬀ L = Γ eﬀ R since as we know theoccupations of the two cavities become equal ( n L = n R ).On the other hand for ﬁnite pump/loss asymmetry thereis a ﬁnite current ﬂowing from L to R and therefore anintra-dimer dissipation. This is interesting since the twocavities are only coupled by a coherent hopping coupling.As a result of this ﬁnite current and dissipation the imag-inary part of the oﬀ-diagonal cavity correlation function0 Figure 10. Current ﬂowing in the dimer as obtained from (21)and (22). The shown values of

J/U are the same ones usedfor the panels of Figs. 8–9. has to be diﬀerent from zero, both based on the sum-rule in Eq. (21) and on physical intuition. In Fig. 10 weplot the average current versus

J/U and compare it withthe integral over Im C LR ( ω ) to conﬁrm the quantitativeagreement. We also see that the overall current, althoughvery small, increases with J/U , an eﬀect which does notappear clearly from the shape of Im C LR ( ω ) in Fig. 9 butthat is consistent with the idea that delocalization leadsto more coherent exchange of excitations between the twocavities and therefore an increased current. VIII. CONCLUSIONS

In this article we have studied an open Bose-Hubbarddimer and investigated the possible signatures of a dissi-pative localization-delocalization transition or crossover,where upon tuning the ratio of coherent hopping ver-sus local interaction an initial population imbalance iseither trapped in one of the two cavities (self-trapping)or equally distributed across the dimer.In the semiclassical limit of many photons per site, thatwe reviewed for completeness in Sec. IV, this transition isknown to occur sharply for a purely conservative (Hamil-tonian) dynamics and to remain present in the form ofa short-time dynamical transition in presence of pumpsand losses, while turning into a smooth crossover at longtimes.In the full quantum regime the situation is particu-larly interesting since it is known that in absence of anyasymmetry in the system parameters the stationary statedensity matrix is independent of any Hamiltonian cou-pling and only set by the pump and loss coeﬃcients.To address therefore possible signatures of a dissipativeself-trapping crossover one is forced to go beyond simplesteady-state observables or to explicitly break the sym-metry between the two cavities. To this extent we haveexactly solved the problem by numerical diagonalizationof the Lindbladian superoperator and obtained the sta-tionary state, the full dissipative quantum dynamics andproperties of the excitations on top of the stationarystate, as encoded in the single-particle Green’s functions,see Sec. III . In Sec. V we have shown that the short-time dissipativedynamics shows clear signatures of a crossover betweena localized behavior with ﬁnite residual imbalance andcoherent oscillations leading to a vanishing imbalance,which can be accessed by either changing the ratio

J/U or the initial condition. On the other hand the long-timesdynamics is largely controlled by the dissipative rates. InSec. VI we have shown that by breaking the symmetry ofpump-loss rates between the two cavities one can induce anon-trivial stationary state and a ﬁnite imbalance whichshows a smooth delocalization crossover upon increasing

J/U .Finally, in Sec. VII we have presented our results forthe single particle Green’s functions, in particular thespectral function and the cavity correlation function de-scribing spectrum and occupation of the bosonic modes.These turn out to be sensitive probes of the Hamilto-nian dynamics even in the fully symmetric case, wherethe delocalization crossover is signaled by the splitting ofthe lowest energy single-photon peak into bonding andanti-bonding modes as

J/U is increased. In presence ofa ﬁnite pump-loss asymmetry we have shown that a ﬁ-nite current ﬂows between the left and right cavities andthis has direct consequences in the emergence of a non-vanishing imaginary part of the oﬀ-diagonal cavity cor-relation function.The methodology discussed in this work, based on theexact diagonalization of a few-sites Lindbladian and onthe computation of Green’s functions, can be applied todiﬀerent problems. Within the BHD it would be inter-esting to study the role of two-particle losses recentlydiscussed in the context of the quantum Zeno eﬀect [39–43]. Another future direction is the development of anexact diagonalization Lindblad impurity solver for Dy-namical Mean Field Theory [44–47]; in this scheme theDMFT self-consistent bath is approximated with a lim-ited number of eﬀective sites. In this respect we note thata two-site model turns out to share many similarities [48]with a minimal, yet reasonably accurate, implementationof the DMFT using a single site in the bath [49]. Therationale is simply that, in the dimer, one of the two sitesplays the role of the self-consistent bath for the other.

ACKNOWLEDGMENTS

This work was partially supported by the ANR grant“NonEQuMat” (ANR-19-CE47-0001) (M. Schirò) and byItalian MIUR through the PRIN2017 project CEnTral(Protocol Number 20172H2SC4).

DATA AVAILABILITY

The data that support the ﬁndings of this study areavailable upon reasonable request from the authors.1

Appendix A: Semiclassical Dynamics

The driven-dissipative Bose-Hubbard dimer can bealso analyzed at a semiclassical level, by writing theHeisenberg equations for the cavity ﬁelds in (1) withpumping and losses as non-Hermitian terms and thentaking the expectation values: ˙ a L = − i (cid:104) ( ω L − U ) + 2 U a † L a L (cid:105) a L − iJa R − Γ eﬀ L a L ˙ a R = − i (cid:104) ( ω R − U ) + 2 U a † R a R (cid:105) a R − iJa L − Γ eﬀ R a R where Γ eﬀ L/R = Γ

L/R − P L/R are the eﬀective loss rates,which for single-particle losses must always be positive.By applying the transformation a L/R (cid:43) α L/R e iϑ L/R , α i , ϑ i ∈ R (A1)one can then reduce the equations for the two complexnumbers above into the following three equations for thereal quantities N = n L + n R , Z = n L − n R and φ = ϑ L − ϑ R : ˙ N = − (cid:0) Γ eﬀ L + Γ eﬀ R (cid:1) N − (cid:0) Γ eﬀ L − Γ eﬀ R (cid:1) Z ˙ Z = − (cid:0) Γ eﬀ L + Γ eﬀ R (cid:1) Z − (cid:0) Γ eﬀ L − Γ eﬀ R (cid:1) N − J (cid:112) N − Z sin( φ )˙ φ = − ∆ ω − U Z + 2

J Z √ N − Z cos φ (A2)where ∆ ω = ω L − ω R .

1. Closed System

In the Hamiltonian case, with ∆ ω = 0 for simplicity,the equations reduce to ˙ N = 0 = ⇒ N = N = const . ˙ Z = − J (cid:113) N − Z sin( φ )˙ φ = − U Z + 2

J Z (cid:112) N − Z cos φ (A3)By using the fact that in a closed system the energy isconserved, the two remaining equations can be furtherreduced to a single equation for the macroscopic occupa-tion imbalance: ˙ Z = − (cid:112) p ( Z ) (A4)where p ( Z ) is a polynomial that can be factorized as p ( Z ) = − U (cid:0) Z − Z (cid:1) (cid:0) Z − Z (cid:1) , (A5)with Z the initial imbalance and Z equal to Z = (cid:115) Z + 4 (cid:18) JU (cid:19) (cid:113) N − Z − (cid:18) JU (cid:19) . (A6) Figure 11. Oscillation period T obtained from (A9) as a func-tion of J/U ; the settings are the same as in Fig. 1.

Being under a square root, the sign of p ( Z ) is the realdiscriminant on the evolution of Z . In turn, the sign of p ( Z ) is completely determined by Z being real or imag-inary (since Z is real). If Z is real then the polynomialis positive only in the region between Z and Z and inthe region between − Z and − Z , no matter whether Z is greater or less than Z . If instead Z is imaginary thenthe polynomial is positive only in the region between − Z and Z .The nature of Z is in turn determined by the sign ofthe polynomial (cid:18) JU (cid:19) − (cid:113) N − Z (cid:18) JU (cid:19) − Z . (A7)If we assume that J/U is positive, then the polynomialabove provides a critical JU , given in the main text inEq. (16) that we rewrite here for simplicity (cid:18) JU (cid:19) c = N (cid:32) (cid:112) − ( Z /N ) + 12 (cid:33) . (A8)For JU < (cid:0) JU (cid:1) c Z is real and therefore Z ( t ) oscillates be-tween Z and Z ; for JU > (cid:0) JU (cid:1) c Z is imaginary andtherefore Z ( t ) oscillates between − Z and Z . Then (cid:0) JU (cid:1) c , in this sense, can be interpreted as a critical valuefor a transition from a localized regime (low J ) to a de-localized regime (high J ).This transition can also be seen through the divergenceof the oscillation period at the critical point (Fig. 11),which can be analytically expressed as T = K (cid:18)(cid:16) Z Z (cid:17) (cid:19) U √ − Z (cid:0) JU (cid:1) > (cid:0) JU (cid:1) c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:20) K (cid:18)(cid:16) Z Z (cid:17) (cid:19) − F (cid:18) sin − (cid:16) Z Z (cid:17) , (cid:16) Z Z (cid:17) (cid:19)(cid:21) U √ − Z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:0) JU (cid:1) < (cid:0) JU (cid:1) c (A9)where F ( ϕ, m ) = (cid:82) ϕ du √ − m sin u and K ( m ) = F ( π , m ) are respectively the incomplete and the com-plete elliptic integral of the ﬁrst kind.The divergence is logarithmic, as one can infer by ap-proximating the integral around the critical point ( Z → + ). The fact that the period diverges, making the os-cillations slower and slower, is a common signature of aphase transition and it’s called critical slowing down .

2. Open System

The question is now how much of the non-dissipativeanalysis done above survives in the presence of losses,at intermediate times. As we cannot go further with ananalytical treatment, we have to go back to (A2) andsolve the full system of equations.Intuitively we expect to see a similar oscillatory be-havior of Z ( t ) in the dissipative case, though the meanvalue approaches zero at large enough times since, semi-classically, dissipative cavities decay to vacuum at thestationary state.Indeed, you see that the presence of dissipation has thedouble eﬀect of increasing the oscillation period and pro-ducing an overall decay of the occupation imbalance withtime. But more interestingly, it stimulates a dynamicaltransition from the regime in which the imbalance oscil-lations are between Z and Z to a regime in which theimbalance oscillates around . Appendix B: Analytical Quantum Results at U = 0

1. Green’s Functions

The Green’s functions at U = 0 can be obtained ana-lytically via the Keldysh formalism. Here we start withthe single-cavity Green’s function and then extend to twocoupled cavities. a. Single Cavity The retarded, advanced and Keldysh components ofthe Green’s function are: G R/A ( ω ) = 1 ω − ω ± i (Γ − P ) , (B1) G K ( ω ) = − i (Γ + P )( ω − ω ) + (Γ − P ) . (B2)The loss/pumping rates appear in couple as Γ − P , ex-cept for the Keldysh Green’s function in which they alsoappear as Γ + P . This is a signature of the quantumnature of the system, encoded in the Keldysh Green’sfunction, in the same way that it appears, for example,when adding quantum noise to a semiclassical treatment. b. Two coupled cavities In the case of two coupled cavities, we distinguish be-tween left and right cavity with a subscript

L/R . The uncoupled Green’s functions, denoted with a subscript ,are the ones in (B2) that we’ve found before for the singlecavity, i.e. G R/A i ( ω ) = 1∆ i ± Γ − i , G K i ( ω ) = − i Γ + i ∆ i + Γ − i (B3)where ∆ i (cid:43) ω − ω i , i = L, R (B4) Γ ± i (cid:43) Γ i ± P i , i = L, R (B5)Then the Green’s function components for the left cav-ity are G RL ( ω ) = 1∆ L + i Γ − L − J ∆ R + i Γ − R (B6) G AL ( ω ) = ( G RL ( ω )) ∗ (B7) G KL ( ω ) = − i (cid:34) Γ + L + J Γ + R Γ R + Γ − R (cid:35) (cid:12)(cid:12) G RL ( ω ) (cid:12)(cid:12) (B8)and the corresponding Green’s functions for the right cav-ity are obtained by simply replacing L → R .

2. Steady-State Properties

The retarded Green’s function of the left cavity can bealso rewritten as G RL ( ω ) = ∆ R + i Γ − R ∆ + ∆ − + i (∆ L Γ − R + ∆ R Γ − L ) (B9)where ∆ ± = ω − ω ± and ω ± = ω L + ω R ± (cid:115)(cid:18) ω L − ω R (cid:19) + J + Γ − L Γ − R . (B10)Since the spectral function is proportional to the imag-inary part of the retarded Green’s function, this meansthat the frequency spectrum will be peaked around ω + and ω − , and J will just have the eﬀect of increasing ordecreasing the separation between these two peaks.As for the occupations of the two cavities, they canbe calculated via (19). Analytical expressions can beeasily obtained in some limiting cases. For example, if Γ ± R = Γ ± L , you obtain that (cid:82) + ∞−∞ dω C L ( ω ) = Γ + L / Γ − L and therefore n L ≡ n L = P L Γ L − P L (B11)(and similarly for the right cavity), i.e. the occupation ofthe cavities at the steady-state is equal to the occupationof the uncoupled cavities ( J = 0 ) and it’s completely3 Figure 12. (Top) Steady-state cavity occupations as a func-tion of J for the dimer with loss coeﬃcients (Γ L , Γ R ) =(6 × − , × − ) and pump coeﬃcients ( P L , P R ) = (4 × − , × − ) described in Section VI, at diﬀerent valuesof U . The dashed lines are the theoretical predictions forthe U = 0 case, calculated by combining the exact formulafor the Keldysh Green’s function (B8) and (19). (Bottom)Steady-state imbalance Z = n L − n R corresponding to theoccupations in the top panel. ﬁxed by the pump/loss rates, no matter what the value of J is. This is actually a special case of a result obtained in[35], showing that any number of cavities with the sameincoherent pump/loss rates have a trivial steady statethat does not depend on the details of their Hamiltonian,i.e. in this case neither on J nor on U . This means,in practice, that in order to have non-trivial physics atthe steady state we must have, if not a loss imbalancebetween the two cavities, at least a pump imbalance .A more interesting case is the one at “strong” J , where“strong” means much bigger than at least all the loss coef-ﬁcients. This time, we do not impose any prior conditionon the pump/loss rates. If ω L = ω R for simplicity, thenthe steady-state occupations become n L ≡ n R = Γ − L n L + Γ − R n R Γ − L + Γ − R , (B12)i.e., for strong enough coupling the occupation of the left Note that the quantity Γ − L/R used in the quantum treatmenthas the same value of the semiclassical Γ eﬀ L/R .In addition, below the lasing threshold we can always param-eterize Γ L/R and P L/R as Γ L/R = Γ − L/R (cid:0) n L/R + 1 (cid:1) and P L/R = Γ − L/R n L/R . and of the right cavities are equal and equal to a weightedaverage of their bare occupations.In particular, if the eﬀective losses are equal ( Γ − L =Γ − R ), then n L ≡ n R = n L + n R , (B13)i.e. the steady-state occupation of the two cavities is ex-actly the mean between the bare occupations.The J = 0 and strong J limits match our intuitiveexpectations, i.e. that the occupations of the cavities, asa function of J , start from their uncoupled values andget closer and closer to each other as J is increased, upto the point at which they match each other’s value. Appendix C: Sum-Rules and Particle Currents inthe BHD

We start deriving the sum-rules for the oﬀ-diagonalcorrelation function deﬁned in Eqs. (20)–(21). To thisextent we note that, by our deﬁnitions in Eqs. (13)–(17), (cid:90) + ∞−∞ dω e − iωt C LR ( ω ) = (cid:104) ˆ a L ( t )ˆ a † R + ˆ a † R ˆ a L ( t ) (cid:105) and that by taking the Hermitian conjugate we have (cid:90) + ∞−∞ dω e iωt C ∗ LR ( ω ) = (cid:104) ˆ a R ˆ a † L ( t ) + ˆ a † L ( t )ˆ a R (cid:105) . Taking the t → + limit and the sum/diﬀerence of theabove two equations we obtain (cid:90) + ∞−∞ dω ( C LR ( ω ) + C ∗ LR ( ω )) = 2 (cid:104) ˆ a † L ˆ a R + ˆ a † R ˆ a L (cid:105) as well as (cid:90) + ∞−∞ dω ( C LR ( ω ) − C ∗ LR ( ω )) = 2 (cid:104) ˆ a † R ˆ a L − ˆ a † L ˆ a R (cid:105) , from which the sum-rules quoted in the main text follow.We now relate the average stationary current acrossthe dimer to the pump-loss asymmetry. To this extentwe consider the BHD in Eq. (1) and we start writingdown the quantum equation of motion for the density ofbosons in each site of the dimer, n α ( t ) = Tr (cid:16) ˆ ρ ( t )ˆ n α (cid:17) ,with α = L/R , which read dn L dt = i (cid:68)(cid:104) ˆ T, ˆ n L (cid:105)(cid:69) + 2 (cid:16) P L + n L ( P L − Γ L ) (cid:17) (C1) dn R dt = i (cid:68)(cid:104) ˆ T, ˆ n R (cid:105)(cid:69) + 2 (cid:16) P R + n R ( P R − Γ R ) (cid:17) (C2)where ˆ T = J (cid:16) ˆ a † L ˆ a R + ˆ a † R ˆ a L (cid:17) is the kinetic energy oper-ator. The commutator gives (cid:104) ˆ T, ˆ n L (cid:105) = − (cid:104) ˆ T, ˆ n R (cid:105) = J (cid:16) ˆ a † R ˆ a L − ˆ a † L ˆ a R (cid:17) ≡ i ˆ I . (C3)4If we take the diﬀerence between the two equations weobtain for the dynamics of the imbalance Z = n L − n R the result dZdt = − (cid:104) ˆ I (cid:105) + 2 (cid:16) ∆ P − n L Γ eﬀ L + n R Γ eﬀ R (cid:17) (C4)In the stationary state the right hand side goes to zeroand we obtain (cid:104) ˆ I (cid:105) = ∆ P − n L Γ eﬀ L + n R Γ eﬀ R (C5) from which, using Eq. (B11), we immediately concludethat for symmetric pump and losses there is no aver-age current between the two sites of the dimer and as aconsequence, using Eq. (21), the imaginary part of theoﬀ-diagonal cavity correlation function has vanishing in-tegral. [1] J. M. Raimond, M. Brune, and S. Haroche, Reviews ofModern Physics , 565 (2001).[2] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wall-raﬀ, To appear in: Reviews of Modern Physics (2021),arXiv:2005.12667.[3] H.-P. Breuer and F. Petruccione, The Theory of OpenQuantum Systems (Oxford University Press, 2007).[4] E. M. Kessler, G. Giedke, A. Imamoglu, S. F. Yelin, M. D.Lukin, and J. I. Cirac, Physical Review A , 012116(2012).[5] F. Minganti, A. Biella, N. Bartolo, and C. Ciuti, PhysicalReview A , 042118 (2018).[6] H. J. Carmichael, Physical Review X , 031028 (2015).[7] W. Casteels, F. Storme, A. Le Boité, and C. Ciuti, Phys-ical Review A , 033824 (2016).[8] A. Le Boité, G. Orso, and C. Ciuti, Physical ReviewLetters , 233601 (2013).[9] M. Schiró, C. Joshi, M. Bordyuh, R. Fazio, J. Keeling,and H. E. Türeci, Physical Review Letters , 143603(2016).[10] F. Vicentini, F. Minganti, R. Rota, G. Orso, andC. Ciuti, Physical Review A , 013853 (2018).[11] A. Biella, F. Storme, J. Lebreuilly, D. Rossini, R. Fazio,I. Carusotto, and C. Ciuti, Physical Review A , 023839(2017).[12] O. Scarlatella, R. Fazio, and M. Schiró, Physical ReviewB , 064511 (2019).[13] H. Landa, M. Schiró, and G. Misguich, Physical ReviewLetters , 043601 (2020).[14] H. Landa, M. Schiró, and G. Misguich, Physical ReviewB , 064301 (2020).[15] M. C. Cross and P. C. Hohenberg, Reviews of ModernPhysics , 851 (1993).[16] J. Raftery, D. Sadri, S. Schmidt, H. E. Türeci, and A. A.Houck, Physical Review X , 031043 (2014).[17] C. Eichler, Y. Salathe, J. Mlynek, S. Schmidt, andA. Wallraﬀ, Physical Review Letters , 110502 (2014).[18] K. G. Lagoudakis, B. Pietka, M. Wouters, R. André,and B. Deveaud-Plédran, Physical Review Letters ,120403 (2010).[19] M. Galbiati, L. Ferrier, D. D. Solnyshkov, D. Tanese,E. Wertz, A. Amo, M. Abbarchi, P. Senellart, I. Sagnes,A. Lemaître, E. Galopin, G. Malpuech, and J. Bloch,Physical Review Letters , 126403 (2012).[20] M. Abbarchi, A. Amo, V. G. Sala, D. D. Solnyshkov,H. Flayac, L. Ferrier, I. Sagnes, E. Galopin, A. Lemaître,G. Malpuech, and J. Bloch, Nature Physics , 275(2013). [21] P. Hamel, S. Haddadi, F. Raineri, P. Monnier, G. Beau-doin, I. Sagnes, A. Levenson, and A. M. Yacomotti, Na-ture Photonics , 311 (2015).[22] M. Marconi, F. Raineri, A. Levenson, A. M. Yacomotti,J. Javaloyes, S. H. Pan, A. El Amili, and Y. Fainman,Physical Review Letters , 213602 (2020).[23] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy,Physical Review Letters , 4950 (1997).[24] L. Pitaevskii and S. Stringari, Physical Review Letters , 180402 (2001).[25] A. Polkovnikov, S. Sachdev, and S. M. Girvin, PhysicalReview A , 053607 (2002).[26] M. Albiez, R. Gati, J. Fölling, S. Hunsmann, M. Cris-tiani, and M. K. Oberthaler, Physical Review Letters , 010402 (2005).[27] M. Trujillo-Martinez, A. Posazhennikova, and J. Kroha,Physical Review Letters , 105302 (2009).[28] T. Venumadhav, M. Haque, and R. Moessner, PhysicalReview B , 054305 (2010).[29] T. Pudlik, H. Hennig, D. Witthaut, and D. K. Campbell,Physical Review A , 063606 (2013).[30] T. C. H. Liew and V. Savona, Physical Review Letters , 183601 (2010).[31] M. Bamba, A. Imamoğlu, I. Carusotto, and C. Ciuti,Physical Review A , 021802 (2011).[32] W. Casteels and C. Ciuti, Physical Review A , 013812(2017).[33] K. Seibold, R. Rota, and V. Savona, Physical Review A , 033839 (2020).[34] S. Schmidt, D. Gerace, A. A. Houck, G. Blatter, andH. E. Türeci, Physical Review B , 100507 (2010).[35] J. Lebreuilly, M. Wouters, and I. Carusotto, ComptesRendus Physique , 836 (2016).[36] E. Arrigoni and A. Dorda, in Out-of-Equilibrium Physicsof Correlated Electron Systems , Springer Series in Solid-State Sciences, edited by R. Citro and F. Mancini(Springer International Publishing, 2018) Chap. 4, pp.121–188.[37] O. Scarlatella, A. A. Clerk, and M. Schiro, New Journalof Physics , 043040 (2019).[38] D. Sarchi, I. Carusotto, M. Wouters, and V. Savona,Physical Review B , 125324 (2008).[39] B. Misra and E. C. G. Sudarshan, Journal of Mathemat-ical Physics , 756 (1977).[40] A. Peres, American Journal of Physics , 931 (1980).[41] W. M. Itano, D. J. Heinzen, J. J. Bollinger, and D. J.Wineland, Physical Review A , 2295 (1990).[42] N. Syassen, D. M. Bauer, M. Lettner, T. Volz, D. Dietze, J. J. Garcia-Ripoll, J. I. Cirac, G. Rempe, and S. Durr,Science , 1329 (2008).[43] D. Rossini, A. Ghermaoui, M. B. Aguilera, R. Vatré,R. Bouganne, J. Beugnon, F. Gerbier, and L. Mazza,arXiv preprint (2020), arXiv:2011.04318.[44] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozen-berg, Reviews of Modern Physics , 13 (1996).[45] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Reviews of Modern Physics , 779 (2014).[46] E. Arrigoni, M. Knap, and W. von der Linden, PhysicalReview Letters , 086403 (2013).[47] O. Scarlatella, A. A. Clerk, R. Fazio, and M. Schiró,arXiv preprint (2020), arXiv:2008.02563.[48] M. Capone and S. Ciuchi, Phys. Rev. B , 104409(2002).[49] M. Potthoﬀ, Physical Review B64