Many-body out-of-equilibrium dynamics of hard-core lattice bosons with non-local loss

We explore the dynamics of hard-core lattice bosons in the presence of strong non-local particle loss. The evolution occurs on two distinct time-scales, first a rapid strongly correlated decay into a highly degenerate Zeno state subspace, followed by a slow almost coherent evolution. We analytically solve the fast initial dynamics of the system, where we specifically focus on an initial Mott insulator state, and perform an analysis of the particle arrangements in the Zeno subspace. We investigate the secondary slow relaxation process that follows and find an intricate regime where the competition between dissipation and coherence results in various types of interacting particle complexes. We classify them and analyse their spectral properties in the presence and absence of nearest-neighbor interactions. Under certain circumstances the dispersion relations of the complexes feature flat bands, which are a result of an effective spin-orbit coupling.

Introduction.-The out-of-equilibrium behaviour of open quantum many-body systems is currently under intense investigation [1][2][3][4][5]. This interest is rooted in the fact that often the competition between coherent and incoherent processes gives rise to seemingly counterintuitive phenomena. Examples are the creation of entanglement by dissipation [6][7][8][9][10][11][12][13][14] and the emergence of effective interparticle interactions [15,16]. In certain cases the latter may even lead to a binding mechanism [1,4], which is qualitatively different to the one resulting from coherent forces that bind constituent particles, in for example molecules or atoms [17]. In Ref. [1] the creation of dissipatively bound complexes was shown to be due to the quantum Zeno effect [7,[18][19][20][21][22], i.e. due to strong dissipation preventing the occupation of particular states by projecting the system onto a reduced state space, the Zeno subspace. While this leads to a good understanding of the few-body physics, a systematic exploration of out-of-equilibrium dynamics on the many-body level is so far lacking.
The purpose of this work is to provide insight into many-body dynamics resulting from a competition between coherent particle motion and strong non-local particle loss through primarily analytic analysis. To this end we consider the situation of a one-dimensional lattice filled with hard-core bosons in a Mott insulating state. We find that the evolution proceeds in two stages. The first stage is characterized by a purely dissipative dynamics that leads to a strongly correlated loss of bosons until the system reaches a highly degenerate Zeno subspace. The second stage is governed by the competition between the dissipation and coherent particle hopping that leads to the formation of dissipatively bound complexes. We identify two qualitatively different types which naturally occur in the Zeno subspace. Their dispersion relations depend strongly on the number of constituent bosons and we find for some configurations the emergence of so-called flat bands [23] which result from an effective spin-orbit coupling and gives rise to immobile complexes [24]. Such flat bands are of interest in the study of exotic topological states of matter e.g. in quantum Hall physics [25]. We further analyze the effect of interactions among neighboring bosons and between complexes.
System.-We consider a one-dimensional lattice with N sites filled with hard-core bosons [16], a scenario which can for example be realized with optically trapped cold atoms [26]. Bosons tunnel between adjacent sites at a rate J such that the Hamiltonian is given by . Here σ ± k = (σ x k ± iσ y k )/2, with {σ x , σ y , σ z } being the standard Pauli matrices. In addition to the Hamiltonian evolution we consider non-local dissipation which is given by distance-selective pair loss, meaning that two bosons separated by the critical distance R are ejected from the lattice at a rate γ [see Fig.  1]. This type of dissipation can be physically realised in cold atoms experiments by exploiting the properties of high-lying excited states, so-called Rydberg states, as shown in Ref. [1]. The dynamics of the density matrix ρ of the system is described by a master equation in Lindblad form, In this work we focus on the limit of strong dissipation, e.g. γ J. This leads to a separation of the two timescales on which the coherent L c and dissipative L d dynamics proceed.
Fast dissipative dynamics and the Zeno subspace.-We begin by analyzing the fast dissipative dynamics. Its stationary subspace -the Zeno subspace -is spanned by all states |s that satisfy L j |s = 0 ∀ j, i.e. they do not contain any two bosons at the critical distance R. To understand the dissipative non-equilibrium evolution into the Zeno subspace we consider our system starting in a Mott insulator state. The corresponding evolution of the average boson density p(t) = j n j (t)/N , with n j = σ + j σ − j , can be found analytically: The mean value of the density on site j evolves under the fast dynamics of L d according to˙ n j = −γ( n j n j+R + n j−R n j ), i.e. it depends on a two-point correlation function. Defining the correlators C k = k l=0 n j+lR and using translational invariance we obtain the hierarchyĊ k = −γ(kC k +2C k+1 ). This equation can be solved by introducing (see Ref.  (a) Evolution of the boson density p(t) under the dissipative dynamics L d from an initial Mott insulator. The stationary density is p(t → ∞) = e −2 ≈ 0.14. The inset shows a single trajectory with 40 bosons and periodic boundary conditions. For this particular trajectory 16 jumps (loss processes) occur, such that the final state contains solely 8 bosons. (b) Representative boson arrangements in the stationary state, where single free bosons and two types of particle complexes can emerge. The circles indicate sites, a filled circle indicates an occupied site, a cross indicates a site whose occupation is forbidden, as the resulting configuration would not lie in the Zeno subspace, and a box indicates the "size" of a complex. The type I complex -defined as having a size smaller than R -is, in this example, constituted of two bosons. These bosons are unable to tunnel away from each other without running into a forbidden site which leads to an effective binding. The type II complex has a spatial extent that is larger than R. It is qualitatively different to type I in the sense that the removal of one boson (in the center) destroys the binding for the remaining ones. (c) Probability distributions for single bosons, type I and type II complexes in the stationary state that is reached from a Mott insulator.
For a Mott insulator state we have the initial condition C k = 1 and therefore G(x, 0) = ∞ k=0 x k /k! = e x . With this, the general solution becomes G(x, t) = e (2+x)e −γt −2 and thus the density evolves as p(t) = C 0 = G(x = 0, t) = e 2(e −γt −1) . Numerical Monte Carlo simulations [see Fig. 1(a)] confirm the rapid exponential decay of the boson density on a timescale ∼ γ −1 . The inset shows a generic trajectory which displays the fast removal of boson pairs and a stationary configuration in which boson pairs at distance R are absent. This is one configuration of many that span the high dimensional stationary Zeno subspace, the projector onto which can be explicitly written as Q 0 = N j=1 (1 − n j n j+R ). The average density in the stationary state reached from a Mott insulator is given by p(t → ∞) = e −2 ≈ 0.14. Note, that this calculation is in fact independent of the value of R Effective coherent dynamics in the Zeno subspace.-Once having reached the Zeno subspace the dissipative evolution governed solely by L d comes to a halt. However, in the presence of quantum tunneling, due to L c , non-trivial coherent dynamics emerges which takes place on a timescale J −1 . As shown in Ref. [1] the effective master equation for the projected density matrix onto the Zeno subspace, µ ≡ Q 0 ρQ 0 , in the limit γ J, is obtained by means of Kato perturbation theory [16,28]: with the effective decay rate Γ = 2J 2 /γ and the operators By construction the dynamics under H Z is restricted to the Zeno subspace. Dissipation within the Zeno subspace affects boson pairs (L j,2 ) in configurations that are "one tunneling event away" from containing bosons at the critical distance R. Such configurations undergo an incoherent evolution at a rate Γ, which is strongly suppressed for fast two-body decay γ J. Therefore the evolution within the Zeno subspace becomes predominantly coherent.
Families of coherent particle complexes.-The approximately coherent evolution under H Z has interesting consequences. Due to the explicit appearance of the projector Q 0 , the simultaneous occupation of two sites at a distance of R is forbidden. This leads to strong correlations and the formation of bound complexes. These complexes can contain a variable number of bosons, but there are two qualitatively different configuration sets in which they can form. Let us start with the simplest case -referred to as type I -aspects of which were already discussed in Ref. [1]. Here m bosons are localized in a region with spatial extent smaller than R, an example of which is shown in Fig. 1(b). These bosons are effectively bound since they cannot separate by more than R − 1 sites under the evolution governed by H Z . The second class -type II -are distinguished by having a spatial extent greater than R. These complexes can form when the bosons and their associated critical distances overlap [see Fig. 1(b)]. Here, unlike for type I, not every particle binds all the others, but one can even encounter situations in which the removal of one boson destroys the entire complex, an example of which is shown in Fig.  1(b). Both type I and II complexes appear naturally in the stationary state that is reached from a Mott insulator. Their relative abundance is shown in the histogram presented in Fig. 1(c). Besides single bosons, there is a significant proportion that occupy a type I state and only a small number enter a type II state. In the following we will perform a detailed investigation of their properties.
Type I complexes.
-We limit our study to the dynamics of a single complex in the lattice, addressing the interactions among complexes later. In the following we will provide three qualitatively different examples: immobile complexes without internal structure, complexes with an internal structure and effective spin-orbit (SO) coupling, and complexes whose dispersion relations feature a flat band arising from this effective SO coupling.
We start with the simplest type I state: two bosons and a critical distance R = 2. The only possible configuration of these bosons, in a type I state, is to be adjacent. Thus, the basis states are |j, 1 = σ + j σ + j+1 |Φ , where |Φ is the vacuum state. In this notation j denotes the position of the complex and the second index labels the "internal state" of the complex. The projected Hamiltonian H Z in this basis is identically zero. Hence the basis states are trivially eigenstates and |j, 1 represents immobile type I complexes. These type I solutions emerge whenever R = m.
In order to see some non-trivial physics we require a complex with some "internal states". The simplest case of this is constituted by 2 bosons with R = 3, previously discussed in [1]. In order to calculate the spectrum of this complex, a basis of the internal states is defined as |j, 1 = σ + j σ + j+1 |Φ and |j, 2 = σ + j σ + j+2 |Φ . We may also define a creation operator |j, α ≡ b To obtain the corresponding dispersion relation ε ± (K) [see Fig. 2(a)] and eigenstates |K ± , we perform a discrete Fourier transform, using periodic boundary conditions and find: ε ± (K) = ±2J cos q K 2 , |K ± = 1 √ 2N j e ijq K [|j, 2 ± e −iq K /2 |j, 1 ], where q K = 2πK/N is the quasi-momentum. We see that the internal state of the complex is strongly linked to its motion on the lattice, namely the group velocity of the internal states is always in the opposite direction for the same quasi-momentum. This is what we term as effective SO coupling. Note that this spectrum has a degeneracy or crossing that occurs at q K = π.
Lastly we consider a complex where the effective SO coupling results in a flat band, namely the case of two bosons with R = 4. We define a basis with three internal states as: |j, 1 = σ + j σ + j+1 |Φ , |j, 2 = σ + j σ + j+2 |Φ and |j, 3 = σ + j σ + j+3 |Φ . The resulting dispersion relations [shown in Fig. 2(a)] and eigenstates are given by This complex has three branches labelled by η = Using one of these states as the initial condition and propagating it under the full master equation we find indeed that it remains immobile as shown in Fig. 2(b). Note, that the boson density is slowly decreasing on a timescale Γ −1 . This clearly shows that the flat bands are not an artifact of infinitely strong dissipation but instead that they indeed have a drastic effect on the boson dynamics in a system with competing coherent and dissipative evolution.
Let us make some general remarks on the emergence of flat bands in case of type I complexes: For complexes consisting of two bosons, flat bands exist provided that R is even. Furthermore, we find that for two, three and four bosons a flat band emerges when R/m ∈ N. Interactions among bosons also play an important role. In order to illustrate this we consider nearest-neighbor interactions of the form H nn = V j n j n j+1 which might, for instance, emerge in cases where non-local loss is engineered via Rydberg dressing (see Ref. [1]). Such interactions modify the dispersion relations as shown in Fig. 2(a) in the sense that they lift the degeneracy point observed for R = 3, Type II complexes.-We now move our study to type II complexes, i.e. complexes that are larger than R. We give two examples, one without internal structure and one with effective SO coupling.
First we consider three bosons and a critical distance R = 3. The only possible type II complexes have the basis |j, 1 = σ + j σ + j+2 σ + j+4 |Φ . They are immobile -similar to the first type I example -as each boson's movement is inhibited by its the nearest bosons. This is confirmed as well by numerical exact simulations as shown in Fig.  3(a). Such immobile states can be straight-forwardly generalized to larger boson numbers, e.g. in the given example by attaching bosons to either end of the complex keeping a separation of one site.
In the second example we consider four bosons and a critical distance R = 4.
The resulting complex has five internal states: |Φ and |j, 5 = σ + j σ + j+2 σ + j+5 σ + j+8 |Φ and the dispersion relations shown in Fig. 3(b): One is given by ε 0 (K) = 0 and the other four are ε η,δ (K) = η 3 + δ 5 + 4 cos(q K ), with η, δ = ±. Hence, this type II complex features a flat band and spatially localized states of the form |F Again let us conclude with some more general remarks: A flat band of similar structure exists for five bosons with R = 4. For R = 3 and 4, a flat band exists provided the number of bosons is equal to or greater than R. The dispersion relation of this type II complex is not modified by the presence of nearest neighbor interactions. This is due to the fact that given the arrangement of the bosons, the simultaneous occupation of neighboring sites is forbidden. Thus, the flat bands of certain type II complexes are in this case protected from interaction effects in contrast to the type I case.
Interaction between complexes.-As can be seen in the inset of Fig. 1(a) complexes are typically not isolated in the stationary subspace of L d . Hence, interactions between complexes, and complexes and single bosons occur. Given the abundance of each species [see Fig. 1(c)] the latter is the most common scenario. An example for such an interaction is given in Fig. 4. Here we display a single boson in the wave packet state |G = (1/ √ 2πσ 2 ) j e −iq0j e (j−j0) 2 /2σ 2 |j , where j 0 , q 0 , σ are the initial central position, quasi-momentum, and width of the wave packet, respectively, impinging an immobile type I complex with R = 2. In much the same way that the dissipation acts to bind the bosons, it results in a hard core exclusion interaction between isolated bosons and complexes that in this example extends over R sites. In the case at hand this leads to an elastic collision with the type I complex essentially acting as a hard boundary. Using this mechanism one could imagine a situation where two immobile complexes enclose a boson, thereby acting as a trap.
More generally the range of the exclusion interaction is dependent on the internal state of interacting complexes. For the type I complex of two bosons with R = 3 we define an effective complex-complex interaction as H and Θ(x) is the Heaviside step function (see Fig. 4 for an illustration).
Outlook.-In the future it will be interesting to study the quantum phases that emerge in systems that contain solely a single kind of complex, e.g. ones that feature state-dependent interactions and flat bands, and look at the case of a fermion system with equivalent dissipation. Such pure systems could be experimentally prepared in the ultra cold atoms lattice experiments discussed in Refs. [29,30].
Acknowledgements We gratefully acknowledge insight-ful discussions with J.P. Garrahan regarding the initial fast relaxation dynamics and also Matteo Marcuzzi, Sam Genway and James Hickey regarding the dispersion relations. We furthermore thank Robin Stevenson for critical reading and comments on the manuscript.

Derivation of the effective master equation
The effective master equation models the dynamics on Zeno subspace. We derive this effective master equation using the Kato resolvent method [16,28]. The form of our particular L d allows us to decompose it into a set of eigenvalues, k i , and eigenspaces or pseudo-projectors, P i , These projectors form a complete orthogonal set, The Zeno subspace has a corresponding zero eigenvalue, removing it from the expansion. Subbing Eq. (1) into the master equation we geṫ where ρ i = P i ρ and λ omits the steady state space. As the steady state space is the one of interest we define the projector onto the irrelevant space as Q = λ P λ . We split Eq. (4) into the evolution of the relevant and irrelevant spaces by applying the respective projectors: Qρ =QLQρ + QLρ 0 , where L = L d + L c . Formal integration of Eq. (6) gives We assume that we start entirely in the steady state space i.e. Qρ = 0. Expanding L we show that Eq. (7) becomes Which is substituted into Eq. (5) to givė Taking a Laplace transform of this equation gives We then use the fact that γ J, implying that the amplitudes of the Liouvillians compare as L d L c . This allows an expansion of (s − QL) −1 to second order: We then perform an inverse Laplace transform to givė (12) Expanding L d in terms of its projectors and expanding the exponential, we finḋ By integration by parts, this remaining integral is reexpressed aṡ Due to L d is a purely dissipative Liouvillian, the k λ 's are all negative. As we are interested in the long time limit, t 1/γ, the second term is considered negligible, as is the remaining integral due to it is of higher order in J/γ as dρ0(τ ) dτ ∝ J. Leaving an effective master equation with the formρ

Derivation of the Projected Hamiltonian and Jump Operators
Due to the form of Eq. (15), we are only interested in states which are coupled to the Zeno subspace via a single tunnelling event. This leads us to only study the cases of a single pair and a double pair, which shares the central boson, at the critical distance R. We define the forms of the pseudo-projectors, P i , of L d on this truncated space as: where: Q 0 was introduced previously and projects onto no pairs, Q 1 projects onto a single pair and Q 2 projects onto two pairs which share the central boson. The first projector P 0 is the steady state space of L d , P 0 = lim t→∞ L d , it includes only states with no pairs of bosons at the critical distance. The next two higher order projectors, P 1 and P 2 include states with a single pair and a double pair which share a central boson. It can be checked that P 0 , P 1 and P 2 project onto the eigenspaces of L d with eigenvalues 0, −γ/2 and −γ respectively. The exact derivation of the projected Hamiltonian from the first term of (15) relies on the assumption that the system starts in the steady state space, meaning that we reduce ρ 0 = P 0 ρ = Q 0 ρQ 0 , and the property of the Q's, Q i Q j = δ i,j Q i , allowing it to be found by the following method Giving the form of H Z as quoted. We then formulate the projected jump operators from the second term of (15). We first rewrite this term as: Which splits into two equations corresponding to the k λ eigenvalues − 2 γ (Q 0 HQ 1 HQ 0 ρ 0 + ρ 0 Q 0 HQ 1 HQ 0 Upon expansion of Q 0 HQ 1 and Q 0 HQ 2 we find a Lindblad form with the jump operators shown.

Derivation of dispersion relations
To demonstrate how the dispersion relations are calculated the single example of a type I state with 2 bosons for R = 4 will be shown. As stated, each site has an associated set of internal states, |j, α , where {α} = 1 → 3. We define a state of the system as |ψ(j) = [A(j)|1 + B(j)|2 + C(j)|3 ]|j and perform a Fourier transform on this state to give the external quasi-momentum states |K = (1/N ) j e ijq K [A(K)|1 + B(K)|2 + C(K)|3 ]|j We rewrite the projected Hamiltonian in this basis as [|j, 1 j, 2| + |j + 1, 1 j, 2| + |j, 2 j, 3| + |j + 1, 2 j, 3| + H.c.]. (25) Applying this to the |K state it is shown that you are it reduces to an operator on the spin structure: (26) Solving for the eigenvalues and eigenvectors of this matrix yields the results shown for the dispersion relations of this complex.