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

We investigate signatures of a self-trapping transition in the driven-dissipative Bose Hubbard dimer, in presence of incoherent pump and single-particle losses. For fully symmetric couplings the stationary state density matrix is independent of any Hamiltonian parameter, and cannot therefore capture 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 after the system is prepared in a variety of initial states, and on the frequency-resolved spectral properties of the steady state, as encoded in the single-particle Green's functions. We find clear signatures of a localization-delocalization crossover as a function of hopping to interaction ratio. We further show that a finite a pump-loss asymmetry restores a delocalization crossover in the steady-state imbalance and leads to a finite intra-dimer dissipation.


Recent years have seen an increase of interest in open
Markovian quantum systems, which describe a number of experimental platforms for quantum information processing and quantum simulation, both in the realm of atomic physics and quantum optics as well as in the solid state framework. Among these we can mention for example cavity QED experiments [1] and their analogue with superconducting circuits [2]. Here the basic degrees of freedom, photons and qubits, are inevitably exposed to dissipative processes such as losses and decoherence induced by the environment. The quantum dynamics of Markovian systems is described theoretically within the framework of a Lindblad master equation which encodes the competition between coherent (Hamiltonian) evolution and dissipative processes described by a set of jump operators [3]. Out of this competition one can expect non-trivial stationary states and dynamical behavior to emerge, leading to novel dissipative phase transitions [4,5], both in small systems made by few quantum non-linear oscillators [6,7] as well as in larger arrays [8][9][10][11][12][13][14].
An intriguing question which has recently attracted large interest is to understand what kind of dynamical phenomena can arise in these Markovian quantum systems and their relationship with analogous phenomena in the field of classical non-linear dynamical systems in presence of non-linearities, noise and dissipation [15].
A prototype example in this context is provided by the driven-dissipative Bose-Hubbard dimer (BHD), which can be seen as a toy model of strongly correlated open Markovian quantum systems since it encodes the basic competition between local dissipative processes, interactions and non-local coherent hopping processes.
In the closed isolated case, corresponding to a purely conservative Hamiltonian evolution, the BHD has been extensively studied, in particular its self-trapping, or localization-delocalization [23][24][25][26][27][28][29]. Here, an initial imbalance of particles between the two sites of the dimer is either rapidly redistributed by hopping processes leading to an homogeneous configuration or conserved indefinitely, leading to a self-trapped state below a critical ratio between hopping and interaction. This transition corresponds to a spontaneous breaking of the reflection symmetry between the two sites of the dimer. Open-Markovian extensions of the BHD have been mostly focused on the coherently driven case [17,[30][31][32][33] or, in the case of the related Jaynes-Cummings Dimer model [34], the purely dissipative case in absence of any external pumping.
In this work we theoretically study the drivendissipative BHD in presence of single-particle losses and incoherent single-particle drive. This case is somewhat peculiar, since it is known that for a perfectly symmetric model the stationary state of the problem is completely independent of Hamiltonian parameters and only set by the ratio between pump and losses [35], so it cannot contain any signature of a putative delocalization transition. In order to explore the competition between hopping and interactions in a dissipative setting one has therefore to go beyond the analysis of steady-state observables and focus instead on response functions, or to consider an asymmetry between the two sites of the dimer.
In particular we prepare the system in different initial states and follow the exact quantum dynamics of the model, characterizing also the properties of the station-ary state reached at long times. Furthermore we focus on the spectral properties of the BHD as encoded in the Green's functions which for open-Markovian quantum system, much like their closed system counterpart, contain rich insights on the structure of the single-particle excitations around the stationary state.
The paper is organized as follows. In Sec. II we introduce the BHD model and briefly review some of its properties, while in Sec. III we present details on its numerical solution. In Sec. IV we review the known results about the semiclassical limit and the self-trapping transition in the isolated and dissipative cases. Our results for the quantum dynamics in the symmetric pumping regime are discussed in Sec. V, while those for finite pump/loss asymmetry in Sec. VI. In Sec. VII we present results for the Green's functions of the BHD, while Sec. IX is devoted to conclusions.

II. THE MODEL
We start by considering the Hamiltonian of a Bose-Hubbard dimer (BHD). The model is a paradigmatic interacting lattice model which can be realized in a number of platforms. Our implementation including pumping and losses is naturally realized using optical cavities (see also Sec. VIII). For this reason in the following we will refer to the two lattice sites as cavities and to the bosonic degrees of freedom involved in the physics as photons. The Hamiltonian readŝ H = ω 0 (n L +n R ) + U n LnL +n RnR wheren L =â † Lâ L andn R =â † Râ R are the number operators of the left and the right cavities, respectively. The two cavities have the same resonant frequency ω 0 and Kerr non-linearity U , and photons can hop between the cavities at a rate J.
We can add a simple mechanism for incoherent driving and dissipation at the master-equation level, by using single-particle pump and loss operators. In practice, we describe the driven-dissipative dimer by a reduced density matrixρ that evolves according to the Lindblad master equationρ whereL is the Hermitian part of the evolution, while the dissipa-tive piece reads aŝ with the constraint that P i < Γ i ∀i, as if ∃i : P i > Γ i single-particle jump operators alone are no longer sufficient to provide a correct physical description of the system. In this form, Γ L/R are interpreted as loss rates while P L/R as pumping rates. It is convenient to parametrize them as to distinguish the case in which pump/loss rates are symmetric in the dimer, ∆Γ = ∆P = 0 or asymmetric due to an imbalance of pump and/or losses. In fact it is known [35] that for a Bose-Hubbard lattice with uniform parameters and identical single-particle pump and loss rates, i.e. ∆Γ = ∆P = 0 the structure of the stationary state density matrix is particularly simple and readŝ where |N is a Fock state with N bosons and π N ∼ (P/Γ) N up to a normalization factor. We note in the above expression thatρ ss is independent of any Hamiltonian parameter and only set by pump/loss ratio. This implies in particular that the stationary state occupancy n α = Tr (ρ ssnα ) is equal in the two cavities and given by which coincides with the value of an uncoupled Kerr resonator. Given these results, it is clear that any nontrivial dependence from J/U has to be looked for in properties other than the stationary-state observables, as we will discuss in Sec. V and VII A. The above result is however no longer true in presence of a finite asymmetry in the dissipative couplings, leading to ∆Γ, ∆P = 0, as we will see more in detail in Sec. VI and VII B.

III. METHODS
The vectorized version of equation (2) is solved by exact diagonalization, yielding a bi-normalized set of left and right eigenvectors ( l α | and |r α , respectively) that satisfy whereL is the matrix representation of the superoperator L. The cokernel and the kernel 1 ofL are, respectively, the left vacuum I| and the steady-state density matrix |ρ ss . The diagonalization problem can actually be simplified by realizing that both the Hamiltonian and the dissipator posses a global gauge symmetry, expressed by an operator functionalK that commutes withL and that acts aŝ K• = −i N , • . By exploiting this symmetry the matrix L can then be written in a block-diagonal form, where each block is labeled by the eigenvalues ofK.
The matrixL and its eigenvectors are written in a basis of Fock states, with a cutoff N cutoff on each particle number. We've fixed N cutoff = 20 throughout the work as a good compromise between accuracy and time and memory costs; this cutoff guarantees that the error on the displayed average steady-state occupations is equal or below 2%, while higher but more expensive cutoffs would not visibly change the results on the Green's functions.
We note that the form of the Lindblad equation ensures ReL α ≤ 0 ∀α, which prevents the dynamics from unbounded growth with time. Again, if we can exploit the global gauge symmetry, then it is sufficient to diagonalize just the largest diagonal block of the Lindbladian. The knowledge of the time-evolution of the density matrix can then be used to calculate the time-evolution of other observables, for example the occupations of the two cavities (i = {L, R}):

B. Källén-Lehmann Spectral Representation of Green's Functions
Albeit not necessary if one only wants to calculate the steady-state density matrix |ρ ss , the full knowledge of the spectrum can be used to explore the Green's functions of the system. In fact, one can obtain frequency-domain expressions for the retarded and the Keldysh components of the steady-state Green's function defined respectively as where the average is taken over the stationary state and the operator A is evolved with the Lindbladian of the system. Upon inserting a complete set of left and right eigenvectors of the Lindbladian and going to the frequency domain by defining G , we obtain a spectral representation of those functions: We see that the Green's functions of an open Markovian quantum system can be generically written as sum of simple poles at complex frequencies given by the eigenvalues of the Lindbladian and with weights, in general complex, given by the transition matrix elements between the stationary state and some excited state of the system [36,37].
From the practical point of view, if one focuses on the single-particle Green's functions, the calculation can be further simplified via the block-diagonal structure of the Lindbladian outlined above. In fact, since the calculation of the single-particle Green's functions involves states that differ at most by one particle from the stationary state, it turns out that the full knowledge of the spectrum is not necessary; it is sufficient to diagonalize just the 3 largest blocks of the diagonal-block structure.
Assuming that the diagonlization scales as the cube of the matrix linear dimension, this yielded a theoretical 10 4 speedup of the diagonalization with the 20-bosons cutoff we have used in both cavities, as well as a 99.7% reduction of the memory required to store the results.

IV. REVIEW OF SEMICLASSICAL DYNAMICS AND SELF-TRAPPING TRANSITION
In order to have a reference point for the analysis of our results we can start by recalling the predictions of a semiclassical treatment of the quantum dynamics for the BHD [23,38]. This is obtained by writing the exact equations of motion for the cavity field operatorsâ L/R and by closing them by takingâ 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 large photons number, while in the quantum treatment we are typically interested in a few-photons treatment. The resulting equations of motion reaḋ where Γ eff L/R = Γ L/R − P L/R are the effective loss rates, which for single-particle losses must always be positive.
As discussed in more detail in Appendix A, it's possible to write semiclassical equations for the total number of photons N = n L + n R and for the occupation imbalance between the two cavities Z = n L − n R , with n L/R = |α L/R | 2 .
In the closed-system case, corresponding to Γ eff L/R = 0, number and energy conservation yield simplified analytical results for the imbalance Z, predicting a transition from a regime in which Z oscillates above the initial condition Z 0 to a regime in which it oscillates around 0 (solid lines in Fig. 1) as one increases the value of J/U above the critical coupling which depends on the initial total number of photons N 0 and imbalance Z 0 . This phase transition can be seen as a divergence of the oscillation period ( Fig. 11) or as a sharp decay to zero of the time-averaged imbalance ). The open system case is not analytically solvable, but the numerical solution of the equations for the total number of photons and for the cavity occupation imbalance shows that the closed-system picture is preserved for low enough values of the loss coefficients, with the difference that even oscillations around a value that is different from zero at initial times will eventually transition at long enough times to an oscillation regime around zero during the dynamical evolution ( Fig. 1, bottom panel).
We can define the time at which this dynamical transition happens to be some t cross for which the imbalance Z(t) crosses the value Z = 0 for the first time. If we plot this time as a function of J/U , see top panel of Fig. 2, we expect that for the closed system this time is divergent for values of J/U below the critical value; for the open system, however, this time assumes finite values even be-low the critical point and the critical point itself is at a slightly lower value than its closed-system counterpart ((J/U ) c = 2.88 vs. (J/U ) c = 2.91). The peak structure visible below (J/U ) c for the open system is due to the commensurability between the period of the imbalance oscillations, that is a function of J/U itself, and t cross .
Albeit holding in the limit of large photon number only, these semiclassical results provide a useful hint for the quantities to look at in the quantum case, as well as a point of comparison that highlights the intrinsic differences between the two types of analyses.

V. RESULTS: DISSIPATIVE QUANTUM DYNAMICS
We now move on to discuss the full dissipative quantum dynamics of the BHD introduced in Sec. II. We focus in particular on the occupation imbalance Z(t) = n L (t) − n R (t) between the two cavities, which in the semiclassical limit shows a clear change of behavior as a function of the parameters.
In the following we set ω 0 = 1, U = 0.1 and consider a situation of symmetric pump and loss rates, ∆P = ∆Γ = 0, so that by construction the imbalance is zero at long times. We set the effective losses Γ eff L/R = Γ L/R − P L/R = 1 × 10 −4 and the pump P L/R = 2 × 10 −4 , such that the 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 function from J/U , at a fixed initial condition which we take to be a Fock state |3, 1 , corresponding to an initial imbalance Z 0 = 2 and an initial number of photons N 0 = 4. At the semiclassical level, see Eq. (16), this would correspond to a critical coupling (J/U ) c = 3.73 for the selftrapping transition.
In the inset of Fig. 3 we plot the time-dependent imbalance Z(t) for different values of J/U . We find a clear crossover as the hopping is increased, from a pure exponential decay to zero at small J/U = 0.1, to an underdamped decay with fast oscillations superimposed at J/U = 0.26 which evolves further into strongly anharmonic oscillations at large values of the hopping, whose frequency grows with J/U . We can interpret this behavior as a signature of the self-trapping transition in the dissipative quantum dynamics. In the small hopping regime each site of the dimer evolves almost independently and the imbalance goes to zero, while for larger values of the hopping there is a substantial transfer of photons across the dimer, resulting in coherent Rabi-like oscillations, before the imbalance reaches the stationary state.
The J/U dependence can also be studied from the point of view of the time-averaged occupation imbalance Z T . In contrast to the semiclassical case ( Fig. 2), . The cavities have a base frequency ω0 = 1.0. The effective loss is Γ eff L = Γ eff R = 1 × 10 −4 and the pumping rate realizes a steady-state occupation equal to 2 in both cavities, so that Zss = 0 by construction. The semiclassical, non-dissipative critical value of J/U for this particular configuration is (J/U )c ≈ 3.73. The time-averaged occupation imbalance computed over the time interval [0, 1000] is shown in the main panel.
where one expects a sharp transition 2 between Z T = 0 and Z T = 0, in the quantum case we have a smooth crossover between the two regimes. The average imbalance drops quickly with J/U due to the development of damped Rabi oscillations, reaching a minimum around J/U 0.25. Quite interestingly, though, we find the appearance of a region in which the imbalance actually increases as a function of J/U before completely dropping to 0 at higher values of J/U . We note that, with respect to the semiclassical case, the localized (selftrapped) phase with Z T = 0 is strongly suppressed and that already for J/U 1.25 the average imbalance is zero. This is consistent with the expectation that quantum fluctuations, included in the exact solution and not properly treated in the semiclassical approach, tend to reduce the broken symmetry phase.
We now discuss the dynamics on longer time scales, where we expect the small dissipative couplings to dominate over the Hamiltonian parameters. To this extent in Fig. 4 we plot the time-dependent imbalance over a broad range of time scales and for different values of J/U . We see a clear separation of dynamical regimes, from a shorttime one -strongly dependent on J/U , as we discussed above -to a longer-time one where the imbalance expo- Evolution of the imbalance Z(t) for the same settings of Fig. 3, shown at longer times and at log scale. The black, dotted line is obtained analytically at J/U = 0; it corresponds to an exponential decay at a rate 2Γ eff . At long times, ln (Z(t)/Z0) fits a straight line; the inset shows the corresponding decay rate as a function of J/U . nentially decays to zero. While naively one could have expected the decay rate to be set only by the dissipative couplings we see in the inset of Fig. 4 that instead it shows a monotonic increase with J/U . Finally, we consider the dependence of the timedependent imbalance Z(t) from the initial condition. To this extent we fix as initial density matrix a pure Fock state ρ 0 = |n 0L , n 0R n 0L , n 0R |, corresponding to an initial imbalance Z 0 = n 0L − n 0R and initial photon number N 0 = n 0L + n 0R , and change the values of n 0L , n 0R . At the semiclassical level, as we see in Eq. (16), there is a critical value of J/U for any N 0 , Z 0 . In order to highlight the difference between the exact quantum dynamics and the semiclassical evolution we fix the value of the hopping to interaction ratio J/U to be always below (J/U ) c (N 0 , Z 0 ), such that at the semiclassical level the system should be localized (self-trapped) at short times for all the chosen initial conditions (see Eq. (16)) and delocalized at longer times (see Fig. 1).
We plot in Fig. 5 the quantum dynamics of the imbalance for different initial conditions. We see that, quite at the opposite of what expected from the semiclassical analysis, the evolution of Z(t) has a strong dependence on the initial state in which the system is prepared. In particular we find both regimes of slow decay to zero of the imbalance (see for example the initial conditions corresponding to |3, 1 or |4, 1 ), indicating localized/selftrapped behavior, as well as regimes of coherent Rabi-like oscillations of the imbalance (see for example the initial conditions corresponding to |3, 2 or |4, 3 ) that we can interpret as signatures of delocalization. This is consistent with the observation made earlier (see Fig. 3) that quantum fluctuations renormalize the critical coupling and fa- , starting from different |n0L, n0R number states at t = 0. The cavities have a base frequency ω0 = 1.0. The effective loss is Γ eff L = Γ eff R = 1 × 10 −4 and the pumping rate realizes a steady-state occupation equal to 2 in both cavities, so that Zss = 0 by construction. The quantity (J/U )c refers to the semiclassical non-dissipative value in (16).
vor the delocalized regime. We conclude therefore that, as in the semiclassical case, the self-trapping crossover can be accessed by changing the initial condition, however we do not explore here the precise dependence of (J/U ) c from the initial state and whether it can be encoded in a simple expression depending only on N 0 and Z 0 as in Eq. (16).

VI. RESULTS: QUANTUM STEADY STATE FOR FINITE PUMP/LOSS ASYMMETRY
In the previous section we have considered the case of a BHD with symmetric pump and loss rates, resulting in a trivial stationary state with zero imbalance for any value of J/U , but with a rich nonequilibrium dynamics.
As we discussed in Sec. II, in presence of a finite pump/loss asymmetry among the two cavities the stationary state becomes more interesting. We can therefore look for signatures of a delocalization crossover, analogous to what we have shown in Fig. 3, directly in observables such as the steady-state occupation or imbalance.
As an example, we consider two cavities with loss coefficients (Γ L , Γ R ) = (6 × 10 −2 , 2 × 10 −2 ) and pump coefficients (P L , P R ) = (4 × 10 −2 , 1 × 10 −2 ), that thus realize steady-state occupations (n 0L , n 0R ) = (2, 1) in the uncoupled limit J = 0 (see Eq. (7)). In Fig. 6 we plot the dependence of the two cavity occupations (top panel) and imbalance (bottom panel) from the hopping to interaction ratio J/U , for different values of U (keeping ω 0 = 1 as unit). We see in the top panel that as J/U is increased the two occupations both converge towards a common value, which is essentially independent from U . The large-J/U limit of the occupations can be obtained analytically by considering the limit U = 0 and results in a weighted average of the two uncoupled occupations (see Eq. (B12)).
As a consequence of the two occupations becoming equal at large J/U we see in the bottom panel that the steady-state imbalance between the two cavities reduces and approaches zero for large enough J/U , a signature of delocalization. We note that increasing U pushes the crossover J/U scale for delocalization to lower values and we expect for U = 0.1 to obtain a behavior comparable with what obtained from the dynamics (see Fig. 3).

VII. RESULTS: GREEN'S FUNCTIONS
A way to get some insights on the system even when the steady-state observables do not depend neither on J nor on U , as in the case of symmetric pump and losses, is to look instead at the single-particle Green's functions. Either by seeing them as the resolvent of the Lindbladian or as response functions that link different states and thus participate in the calculation of transport quantities like the optical transmission, the Green's functions are sensitive to the details of the Lindbladian spectrum, and not only to the zero mode (stationary state), as it appears clearly from the Källén-Lehmann representation discussed in Sec. III B.
In this section we present our results for the Green's function of the BHD, that we obtained from the exact diagonalization of the Lindbladian as discussed in Sec. III. Specifically we consider the single-particle Green's functions, obtained from Eq. (15) with the choice A = a i and B = a † j with i, j = L/R, and in particular the spectral function A ij (ω) and the cavity correlation function C ij (ω), defined as with i, j = L/R. The diagonal components (for i = j) contain information on the local (on-site) spectrum and occupations of the bosonic mode and satisfy the sum rules where n i is the stationary state occupation. The offdiagonal components contain instead information on the delocalized modes across the dimer. In particular the correlation function C LR (ω) has real and imaginary parts which satisfy the sum-rules where is the average kinetic energy in the stationary state while Î = −iJ â † Râ L −â † Lâ R is the average current flowing from L to R (see Appendix C). We now presents our results for these Green's functions, starting from the pump/loss symmetric case and then discussing the role of a finite pump/loss asymmetry.

A. Symmetric Pump and Losses
We start considering the case of symmetric pump and loss rates, ∆Γ = ∆P = 0. As a result the system is completely symmetric upon reflection (L ↔ R) and as such the diagonal spectral functions in Eq. (17) do not depend on the index i = L/R. As an example, in the top panel of Fig. 7 we plot the spectral function of the left cavity for different values of J/U . At low J/U the spectral function resembles much the one of a single driven-dissipative Kerr resonator, with a characteristic sequence of peaks located at frequencies given by the energy difference between states with n + 1 and n photons, ∆ n = E n+1 − E n = ω 0 + U + 2U n, where E n = ω 0 n + U n 2 is the energy of the Kerr resonator with n photons (see the Hamiltonian in Eq. (1)). These peaks, which start at ω 0 + U and are equally spaced by 2U , would be infinitely sharp in the closed system while are broadened by the dissipative processes by an amount roughly given by Γ eff L (it would perfectly match this value in the non-interacting, decoupled case J = U = 0, see Appendix B 1.) As J/U is increased we see that the first effect is the creation of sub-peaks within each resonance, particularly in the low frequency ones, with the center of mass of each band remaining roughly located at the isolated Kerr excitation energies. Upon increasing further J/U we see how different bands start to merge in a continuum and for J/U = 0.64 a new features arises, namely a finite spectral weight appears below the resonator frequency ω 0 = 1, which becomes a sharp peak for large values of J/U (e.g. J/U = 1.50). This peak corresponds to a delocalized photonic excitation as one can realize by looking at the spectral function in the opposite limit of U = 0 (see Appendix B 1), which has two poles at frequencies roughly ω ± ω 0 ± J since in this regime the dissipative couplings are very small.
It is interesting to connect these spectral features to the behavior of the time-dependent and of the time-averaged imbalance shown in Fig. 3 for similar values of J/U . For small values of the hopping the imbalance is different from zero at short and intermediate times, i.e. photons remain localized in one of the two cavities and the spectral function resembles the one of an isolated Kerr resonator. Upon increasing J/U photons start to hop coherently within the dimer: the imbalance shows short-time Rabi oscillations with a period controlled by J/U and its time-average vanishes, while spectrally this translates in the emergence of two peaks above and below the bare resonator frequency.
In the bottom panel of Fig. 7 we plot the real-part of the off-diagonal correlation function, for different values of J/U and ∆Γ = ∆P = 0. We note that quite interestingly the imaginary part of this Green's function vanishes in this regime, a point onto which we will come back in the next section. At small values of the hopping the real-part C LR (ω) is essentially zero, the cavities are almost decoupled, except at frequencies corresponding to the eigenmodes of the (interacting) single cavity (see top panel at the same value of J/U ), where an anti-resonance like contribution emerges. Upon increasing J/U , as we discussed for the spectral function, further peaks appear which start merging and shifting towards lower frequencies. We note that the structure of the peaks evolve as well: at small J/U they are almost perfectly asymmetric in frequency (leading to a vanishing integral, see Eq. (20)) while upon increasing J/U , when the system becomes more delocalized, this asymmetry disappears. Furthermore, also the strength of the peaks increases with J/U (note the different scale in the panels) in a way that appears opposite to the peaks in the spectral function in the top panel. This is again consistent with the idea that upon entering in the delocalized regime the weight is transferred from the localized (on-site) modes to the delocalized (off-diagonal ones).

B. Asymmetric Pump and Losses
We now move to discuss the case of asymmetric pump and losses, ∆P, ∆Γ = 0, resulting as we know in a non trivial stationary state density matrix (and finite imbalance, see Sec. VI). A natural question is whether this different nonequilibrium protocol results in a qualitatively different behavior of the Green's functions.
We start from the spectral functions, that we plot in Fig. 8 for a fixed pump/loss asymmetry and different values of J/U . To highlight the comparison between the two cavities we plot the left and right spectral functions on a common frequency scale. While we see a similar structure of peaks evolving with J/U , as compared to the symmetric case of Fig. 7, we also note an interesting dependence  from the pump/loss asymmetry and the hopping. In particular, for small J/U the right cavity spectral function (bottom panels) has slightly stronger peaks at low frequency than the left cavity one, reflecting the asymmetry in the pump/loss rates. As the hopping is increased and the excitations are delocalized in the dimer we see that this asymmetry in the left/right spectral functions decreases and for J/U = 1.50 the two spectra are essentially the same and very close in shape to the symmetric one for the same value of J/U (See Fig. 7).
Then we consider the off-diagonal cavity correlation function, see Fig. 9, that we study as a function of J/U .
In the top panel we plot the real part, Re C LR (ω), which shows a qualitative behavior very similar to the symmetric case shown in Fig. 7, with anti-Lorentzian peaks which broaden and merge into a continuum at large J/U indicating the increase in kinetic energy. On the other hand, an interesting difference appears in the imaginary part of the off-diagonal cavity correlation function, Im C LR (ω), which is now different from zero and shows a non-trivial dependence from J/U , with narrow peaks which broaden and merge into a continuum as J/U is increased.
We can understand the origin of a finite imaginary part of the off-diagonal cavity correlation function by using the sum rule that relates the integral of Im C LR (ω) to the average current flowing from L to R (see Eq. (21) and Appendix C). In the stationary state the average current is completely determined by the effective pump/loss rates Γ eff L/R = Γ L/R −P L/R and the stationary occupation n L/R through the relation where ∆P is the pump asymmetry. We see that the right-hand side of this equation exactly vanishes in the symmetric case ∆P = 0, Γ eff L = Γ eff R since as we know the occupations of the two cavities become equal (n L = n R ). On the other hand for finite pump/loss asymmetry there is a finite current flowing from L to R and therefore an intra-dimer dissipation. This is interesting since the two cavities are only coupled by a coherent hopping coupling. As a result of this finite current and dissipation the imaginary part of the off-diagonal cavity correlation function has to be different from zero, both based on the sumrule in Eq. (21) and on physical intuition. In Fig. 10 we plot the average current versus J/U and compare it with the integral over Im C LR (ω) to confirm the quantitative agreement. We also see that the overall current, although very small, increases with J/U , an effect which does not appear clearly from the shape of Im C LR (ω) in Fig. 9 but that is consistent with the idea that delocalization leads to more coherent exchange of excitations between the two cavities and therefore an increased current.
Finally, we have also considered the case of extreme pump/loss asymmetry, corresponding to the situation in which one of the two cavities is non-dissipative, i.e. Γ eff R = P R = 0. Quite interestingly we have found that also in this case, as for perfectly symmetric rates, the current and the dissipative part of the off-diagonal cavity correlation function Im C LR (ω) are both zero, for any value of J/U . We can understand this result from a simple physical picture: in absence of a Markovian environment coupled to the right cavity the current flowing from left to right cannot be dissipated and bounces back, resulting in a zero net current. This can be also understood more formally, by looking at Eq. (22) and by noting that for Γ eff R = P R = 0 this reduces to Î = Γ eff L (n 0L − n L ). As we discuss in Appendix B 2 in the limit Γ eff R = P R = 0 the left cavity occupation reduces to the one of an isolated left site coupled to Markovian pump and losses, i.e. n L = n 0L resulting therefore in a vanishing current.  (21) and (22). The shown values of J/U are the same ones used for the panels of Figs. 8-9.

VIII. DISCUSSION
In this section we discuss our results on the BHD in the broader context of driven-dissipative phase transitons and comment more in detail on the experimental realization of our setup and our findings.
As for their closed system counterparts, dissipative phase transitions emerge sharply in the limit of thermodynamically large systems [5,13]. In the open-system context this has been shown to arise when taking the large volume limit at fixed finite-density or in the limit of large photon numbers, correspondingly to a well defined classical limit. From this point of view it is not surprising that for our BHD the localization-delocalization transition that exists at the semiclassical level turns in a crossover in presence of quantum fluctuations. These are in fact particularly strong in the present case where the system size is finite and therefore the Liouvillian gap is non-vanishing. This does not exclude of course the presence of sharp nonequilibrium phase transitions for arrays of driven-dissipative cavities with incoherent pumping, as it has been indeed recently discussed [11,37].
As we mentioned in the introduction, the drivendissipative BHD has been realized experimentally in a variety of quantum light-matter platforms. In circuit QED this can be done by considering the large detuning limit of two coupled Jaynes-Cummings (JC) units, which can be realized by capacitively coupling two resonators, each containing a transmon qubit. In this context the focus has been mostly on the case of coherently driven cavities, or of purely dissipative (lossy) dynamics, however an incoherent pump can be also engineered by weakly coupling each site of the dimer to a transmission line or to an incoherent noise [39]. In an actual experimental setting, the case of perfectly symmetric dimer is obviously more difficult to achieve due to local imperfections which introduce small disorder in the system. This however has been shown to remain controllable, particularly for small lattices [40,41]. Our results for the dynamics of the imbalance or its dependence from external parameters, as well as the Green's functions, can be directly measured experimentally. The former has been done in the context of a JC dimer through homodyne detection [16].
The latter can be naturally addressed in a transmission/reflection experiment. Finally, in other light-matter platforms, such as semiconductor microcavities and photonic crystals, incoherent pumping is even more natural to realize, especially for lasing applications [21]. We also mention the BHD is relevant for ultracold atomic gases experiments with double-well systems, and in this context controlled dissipative (incoherent) processes of pump and losses can be engineered by coupling to other bands.

IX. CONCLUSIONS
In this article we have studied an open Bose-Hubbard dimer and investigated the possible signatures of a dissipative localization-delocalization transition or crossover, where upon tuning the ratio of coherent hopping versus local interaction an initial population imbalance is either trapped in one of the two cavities (self-trapping) or equally distributed across the dimer.
In the semiclassical limit of many photons per site, that we reviewed for completeness in Sec. IV, this transition is known to occur sharply for a purely conservative (Hamiltonian) dynamics and to remain present in the form of a short-time dynamical transition in presence of pumps and losses, while turning into a smooth crossover at long times.
In the full quantum regime the situation is particularly interesting since it is known that in absence of any asymmetry in the system parameters the stationary state density matrix is independent of any Hamiltonian coupling and only set by the pump and loss coefficients. To address therefore possible signatures of a dissipative self-trapping crossover one is forced to go beyond simple steady-state observables or to explicitly break the symmetry between the two cavities. To this extent we have exactly solved the problem by numerical diagonalization of the Lindbladian superoperator and obtained the stationary state, the full dissipative quantum dynamics and properties of the excitations on top of the stationary state, as encoded in the single-particle Green's functions, see Sec. III .
In Sec. V we have shown that the short-time dissipative dynamics shows clear signatures of a crossover between a localized behavior with finite residual imbalance and coherent 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-times dynamics is largely controlled by the dissipative rates. In Sec. VI we have shown that by breaking the symmetry of pump-loss rates between the two cavities one can induce a non-trivial stationary state and a finite imbalance which shows a smooth delocalization crossover upon increasing J/U . Finally, in Sec. VII we have presented our results for the single particle Green's functions, in particular the spectral function and the cavity correlation function describing spectrum and occupation of the bosonic modes.
These turn out to be sensitive probes of the Hamiltonian dynamics even in the fully symmetric case, where the delocalization crossover is signaled by the splitting of the lowest energy single-photon peak into bonding and anti-bonding modes as J/U is increased. In presence of a finite pump-loss asymmetry we have shown that a finite current flows between the left and right cavities and this has direct consequences in the emergence of a nonvanishing imaginary part of the off-diagonal cavity correlation function.
The methodology discussed in this work, based on the exact diagonalization of a few-sites Lindbladian and on the computation of Green's functions, can be applied to different problems. Within the BHD it would be interesting to study the role of two-particle losses recently discussed in the context of the quantum Zeno effect [42][43][44][45][46]. Another future direction is the development of an exact diagonalization Lindblad impurity solver for Dynamical Mean Field Theory [47][48][49][50]; in this scheme the DMFT self-consistent bath is approximated with a limited number of effective sites. In this respect we note that a two-site model turns out to share many similarities [51] with a minimal, yet reasonably accurate, implementation of the DMFT using a single site in the bath [52]. The rationale is simply that, in the dimer, one of the two sites plays the role of the self-consistent bath for the other.

DATA AVAILABILITY
The data that support the findings of this study are available upon reasonable request from the authors.

Appendix A: Semiclassical Dynamics
The driven-dissipative Bose-Hubbard dimer can be also analyzed at a semiclassical level, by writing the Heisenberg equations for the cavity fields in (1) with pumping and losses as non-Hermitian terms and then taking the expectation values: where Γ eff L/R = Γ L/R − P L/R are the effective loss rates, which for single-particle losses must always be positive.
By applying the transformation one can then reduce the equations for the two complex numbers above into the following three equations for the real quantities N = n L + n R , Z = n L − n R and φ = ϑ L − ϑ R : where ∆ω = ω L − ω R .

Closed System
In the Hamiltonian case, with ∆ω = 0 for simplicity, the equations reduce to By using the fact that in a closed system the energy is conserved, the two remaining equations can be further reduced to a single equation for the macroscopic occupation imbalance:Ż where p(Z) is a polynomial that can be factorized as with Z 0 the initial imbalance and Z 1 equal to Being under a square root, the sign of p(Z) is the real discriminant on the evolution of Z. In turn, the sign of p(Z) is completely determined by Z 1 being real or imaginary (since Z 0 is real). If Z 1 is real then the polynomial is positive only in the region between Z 0 and Z 1 and in the region between −Z 0 and −Z 1 , no matter whether Z 1 is greater or less than Z 0 . If instead Z 1 is imaginary then the polynomial is positive only in the region between −Z 0 and Z 0 .
The nature of Z 1 is in turn determined by the sign of the polynomial If we assume that J/U is positive, then the polynomial above provides a critical J U , given in the main text in Figure 11. Oscillation period T obtained from (A9) as a function of J/U ; the settings are the same as in Fig. 1.
Eq. (16) that we rewrite here for simplicity For J U < J U c Z 1 is real and therefore Z(t) oscillates between Z 0 and Z 1 ; for J U > J U c Z 1 is imaginary and therefore Z(t) oscillates between −Z 0 and Z 0 . Then J U c , in this sense, can be interpreted as a critical value for a transition from a localized regime (low J) to a delocalized regime (high J).
This transition can also be seen through the divergence of the oscillation period at the critical point (Fig. 11), which can be analytically expressed as and K(m) = F ( π 2 , m) are respectively the incomplete and the complete elliptic integral of the first kind.
The divergence is logarithmic, as one can infer by approximating the integral around the critical point (Z 1 → 0 + ). The fact that the period diverges, making the oscillations slower and slower, is a common signature of a phase transition and it's called critical slowing down.

Open System
The question is now how much of the non-dissipative analysis done above survives in the presence of losses, at intermediate times. As we cannot go further with an analytical treatment, we have to go back to (A2) and solve the full system of equations.
Intuitively we expect to see a similar oscillatory behavior of Z(t) in the dissipative case, though the mean value approaches zero at large enough times since, semiclassically, dissipative cavities decay to vacuum at the stationary state.
Indeed, you see that the presence of dissipation has the double effect of increasing the oscillation period and producing an overall decay of the occupation imbalance with time. But more interestingly, it stimulates a dynamical transition from the regime in which the imbalance oscillations are between Z 0 and Z 1 to a regime in which the imbalance oscillates around 0.
Appendix B: Analytical Quantum Results at U = 0

Green's Functions
The Green's functions at U = 0 can be obtained analytically via the Keldysh formalism. Here we start with the single-cavity Green's function and then extend to two coupled cavities.

a. Single Cavity
The retarded, advanced and Keldysh components of the Green's function are: The loss/pumping rates appear in couple as Γ − P , except for the Keldysh Green's function in which they also appear as Γ + P . This is a signature of the quantum nature of the system, encoded in the Keldysh Green's function, 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 between left and right cavity with a subscript L/R. The uncoupled Green's functions, denoted with a subscript 0, are the ones in (B2) that we've found before for the single cavity, i.e.
Then the Green's function components for the left cavity are and the corresponding Green's functions for the right cavity are obtained by simply replacing L → R.

Steady-State Properties
The retarded Green's function of the left cavity can be also rewritten as where ∆ ± = ω − ω ± and (B10) Since the spectral function is proportional to the imaginary part of the retarded Green's function, this means that the frequency spectrum will be peaked around ω + and ω − , and J will just have the effect of increasing or decreasing the separation between these two peaks.
As for the occupations of the two cavities, they can be calculated via (19). Analytical expressions can be easily obtained in some limiting cases. For example, if Γ ±R = Γ ±L , you obtain that +∞ −∞ dω C L (ω) = Γ +L /Γ −L and therefore n L ≡ n 0L = P L Γ L − P L (B11) (and similarly for the right cavity), i.e. the occupation of the cavities at the steady-state is equal to the occupation of the uncoupled cavities (J = 0) and it's completely fixed 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 same incoherent pump/loss rates have a trivial steady state that 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 at the steady state we must have, if not a loss imbalance between 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 coefficients. This time, we do not impose any prior condition on the pump/loss rates. If ω L = ω R for simplicity, then the steady-state occupations become 3 i.e., for strong enough coupling the occupation of the left and of the right cavities are equal and equal to a weighted average of their bare occupations. In particular, if the effective losses are equal (Γ −L = Γ −R ), then n L ≡ n R = n 0L + n 0R 2 , i.e. the steady-state occupation of the two cavities is exactly the mean between the bare occupations. The J = 0 and strong J limits match our intuitive expectations, i.e. that the occupations of the cavities, as a function of J, start from their uncoupled values and get closer and closer to each other as J is increased, up to the point at which they match each other's value.
Another interesting limiting case is obtained if one of the cavities, say e.g. the right one, has Γ ±R = 0. Then, for any J, we get n L ≡ n R ≡ n 0L .
(B14) 3 Note that the quantity Γ −L/R used in the quantum treatment has the same value of the semiclassical Γ eff L/R . In addition, below the lasing threshold we can always parameterize Γ L/R and P L/R as Γ L/R = Γ −L/R n 0L/R + 1 and P L/R = Γ −L/R n 0L/R .
In this case the uncoupled occupation of the right cavity, n 0R , is formally ill-defined; however, it can be easily regularized by taking P R = 0 and Γ R = ε, with ε > 0 arbitrarily small, for which n 0R = 0. From a physical point of view, in this case the steadystate occupations in the system are fixed by the only available Markovian environments, i.e. the ones attached to the left cavity, so the occupations become equal as soon as the two cavities are connected (J > 0). For this reason, we expect this result to be valid at U = 0 as well.
Appendix C: Sum-Rules and Particle Currents in the BHD We start deriving the sum-rules for the off-diagonal correlation function defined in Eqs. (20)- (21). To this extent we note that, by our definitions in Eqs. (13)- (17), and that by taking the Hermitian conjugate we have Taking the t → 0 + limit and the sum/difference of the above two equations we obtain +∞ −∞ dω (C LR (ω) + C * LR (ω)) = 2 â † Lâ R +â † Râ L as well as +∞ −∞ dω (C LR (ω) − C * LR (ω)) = 2 â † Râ L −â † Lâ R , from which the sum-rules quoted in the main text follow. We now relate the average stationary current across the dimer to the pump-loss asymmetry. To this extent we consider the BHD in Eq. (1) and we start writing down the quantum equation of motion for the density of bosons in each site of the dimer, n α (t) = Tr ρ(t)n α , with α = L/R, which read whereT = J â † Lâ R +â † Râ L is the kinetic energy operator. The commutator gives If we take the difference between the two equations we obtain for the dynamics of the imbalance Z = n L − n R the result dZ dt = −2 Î + 2 ∆P − n L Γ eff L + n R Γ eff R (C4) In the stationary state the right hand side goes to zero and we obtain Î = ∆P − n L Γ eff L + n R Γ eff R (C5) from which, using Eq. (B11), we immediately conclude that for symmetric pump and losses there is no average current between the two sites of the dimer and as a consequence, using Eq. (21), the imaginary part of the off-diagonal cavity correlation function has vanishing integral.