Liouvillian skin effects and fragmented condensates in an integrable dissipative Bose-Hubbard model

Strongly interacting non-equilibrium systems are of great fundamental interest, yet their inherent complexity make then notoriously hard to analyze. We demonstrate that the dynamics of the Bose-Hubbard model, which by itself evades solvability, can be solved exactly at any interaction strength in the presence of loss tuned to a rate matching the hopping amplitude. Remarkably, the full solvability of the corresponding Liouvillian, and the integrability of the pertinent effective non-Hermitian Hamiltonian, survives the addition of disorder and generic boundary conditions. By analyzing the Bethe ansatz solutions we find that even weak interactions change the qualitative features of the system, leading to an intricate dynamical phase diagram featuring non-Hermitian Mott-skin effects, disorder induced localization, highly degenerate exceptional points, and a Bose glass-like phase of fragmented condensates. We discuss realistic implementations of this model with cold atoms.

A seminal exact solution strategy of a strongly interacting many body system was provided by Bethe [60] in his famous exact solution of the Heisenberg model.Much more recently, third quantization [61] was invented as a framework to solve quadratic (non-interacting) Liouvillians, describing driven-dissipative systems, and in the last few years the notion of integrability was extended to more general systems including certain interacting Liouvillians [62][63][64][65][66].The application of integrability methods to the solution of Liovillians has thus made it possible to construct exactly solvable models that combine both dissipation and strong interactions, but exactly solvable models of that kind are rare, and exactly solvable models that incorporate disorder, dissipation and interactions simultaneously remain elusive.
In this Letter, we present and analytically solve a disordered Bose-Hubbard chain with environmental coupling that is tailored to render it exactly solvable while remaining realistic in cold atom systems (see Fig. 1).In the process we derive and solve the non-Hermitian Hamiltonian of the unidirectional Bose-Hubbard chain that was first shown to be Yang-Baxter integrable in Ref. 40.Here we extend this analysis in several ways that reveal qualitatively new physics.We show how the model emerges as an effective description of a realistic setup, and that it remains integrable in the presence of arbitrary on-site potentials, including random on-site disorder, and with open boundary conditions.The open boundary conditions enable non-Hermitian Mott-skin effects and highly degenerate exceptional points, and adding disorder leads to a novel phase that is reminiscent of a Bose glass.Furthermore, we stress that our analysis goes beyond the effective non-Hermitian Hamiltonian approach by providing a complete solution of the Lindblad Liouvillian.
Model and physical realization.-Dissipativecold atom systems may be realistically modeled using the Lindblad master equation [67,68], with γ j ≥ 0. We consider the coherent dynamics to be described by the disordered Bose-Hubbard Hamiltonian The upper sites have dissipation with rate Γ, while the lower ones do not.A fictitious magnetic flux ϕ leads to a relative phase in the hopping terms between the different sites.Upon integrating out the green sites we obtain (b) with dissipation acting on neighbouring sites, where the magnetic flux has been chosen to be ϕ = π 2 .(c) Shows the system described by the integrable unidirectional effective non-Hermitian (4), where the dissipation has been absorbed into the hopping parameter and the on-site potential.
where a j are bosonic annihilation operators, n j = a † j a j are the number operators, and µ j are disorder variables.
One way to incorporate realistic dissipation so that the system becomes exactly solvable is to consider a triangular ladder [13,14], as in Fig. 1(a).Particles dissipate from the sites on the top and these sites can be integrated out [69] (see also Ref. 70), leaving a chain with jump operators as in Fig. 1(b), cf.Ref. 24.The relative phase factor may be obtained experimentally by including a magnetic field.
Remarkably, at the level of the effective non-Hermitian Hamiltonian, the hopping in one direction vanishes if γ j = t for j ̸ = N and γ N = ϵ 2 .Explicitly, the resulting effective non-Hermitian Hamiltonian (Fig. 1(c)) becomes where we defined µ ′ j = µ j −2it for j ̸ = 1, N and µ . Ref. 14 realized and analyzed the non-interacting and non-disordered limit of Eq. ( 4) in a cold atom system.Here we include both interactions and disorder, and in addition to Eq. ( 4), also consider the full Liovillian quantum dynamics as described by Eqs.
(1), ( 2) and ( 3).An immediate conclusion that can be drawn is that in the absence of interactions a state with N particles decays at a rate of −Im j µ ′ j n j ≈ 2tN .Bethe ansatz solution.-Usingthe machinery of the quantum inverse scattering method [71] it may be shown that the right eigenvectors of the Hamiltonian, and of all its constants of motion defined in the supplemental material [69], have the form given by the algebraic Bethe Ansatz for an operator B(λ) that takes a continuous argument λ.Each application of B creates a quasiparticle, so the state given above has M quasiparticles.In the Fock basis the vacuum |Ω⟩ is equal to the state without any particles in it.The states given by Eq. ( 5) are eigenvectors of the Hamiltonian if the rapidities λ n solve the Bethe equations, which in this case are given by N k=1 In terms of the rapidities the eigenvalue corresponding to the eigenstate (5) of the non-Hermitian effective Hamiltonian is E({λ n }) = M n=1 λ n .In the case without an onsite potential, meaning µ ′ k = 0 for all k, and with ϵ ̸ = 0, this was first shown in Ref. 40.
The coordinate Bethe Ansatz in terms of the rapidities is [72] λ and S M is the set of all permutations of the integers between 1 and M .This expression matches the algebraic Bethe Ansatz.
Integrability and solvability.-Theintegrability structure outlined above can be used to explicitly diagonalize the Liouvillian (1).It is important that there is only dissipation (or only gain), and that the effective non-Hermitian Hamiltonian commutes with the number operator, as this ensures that the Liouvillian is triangular in the basis given by states of the form and the eigenoperators have the form where the coefficients are given in the supplemental material [69] (see also Refs.[74][75][76][77][78].Note that the unique steady state of the system is the empty state, so there cannot be any conserved quantities.The Liouvillian therefore cannot easily be said to be integrable, in spite of the underlying integrability structure that allows its exact analytical solution. Correlation functions and phase diagram for open boundary conditions.-Wenow examine the properties of the model defined by the effective Hamiltonian (4) with random disorder and open boundary conditions (ϵ = 0), which exhibits a remarkable degree of solvability.As we detail below this leads to the phase diagram in Fig. 2. For arbitrary disorder the Bethe equations now take the factorized form U where n = 1, 2, . ... Inserting this into the expression for the energy we see that the spectrum is given by The effective Hamiltonian exhibits exceptional points, i.e. at a non-Hermitian degeneracy at which both eigenvalues and eigenstates coalesce, whenever there are distinct partitions of integers, {n i }, giving identical energies in Eq. (11).In particular this happens when any two (or more) disorder variables, µ i , are equal.Moreover, if the potential is constant (no disorder), there are exceptional points of very high order.It is worth noting that although at t = 0 the integers {n i } labeling energies and eigenstates correspond to particle densities, their interpretation is less straightforward at t ̸ = 0.
In the limit U → 0 it is sufficient to treat the singleparticle eigenstates.These are given by which are eigenstates if λ = µ ′ n for n = 1, 2, . . .N .If there is no disorder, meaning µ ′ n = 0 for all n, all the eigenstates coalesce to ψ(x) = δ x,1 .One simple observation is in that if λ = µ ′ n then the particle density is 0 for x > m where m is the smallest integer such that µ ′ m = µ ′ n .The particle is therefore sharply localized to a finite region.In addition, if t >> |µ ′ m ′ −µ ′ n ′ | for all m ′ , n ′ we see that the particle is exponentially localized in the  15) and (16).For U > 2∆ the ground state is in a Mott phase, which can either be localized at the boundary or be evenly distributed depending on t, whereas for U < 2∆ the system is in Bose-Glass-like phase of fragmented condensates, which may also be localized either at the boundary or at some site.At U = 0 there is a transition between a Bose-Einstein condensate at the boundary and a Bose-Einstein condensate at a single site in the bulk.left side of that region for all eigenstates.On the other hand, if t << |µ ′ m ′ − µ ′ n ′ | on average each eigenstate is exponentially localized to its corresponding site.Extending to the many-body ground state, we see that it is a gapless condensate of particles in the state with λ equal to the µ ′ n with smallest real part.Due to the simplicity of this state, it is possible to place strong bounds on the two-point correlation function.Let x 2 > x 1 both be in the region where the wave function is non-zero.Note that the imaginary parts of the µ n :s cancel in the eigenstates, except for the factors that include µ ′ 1 and µ ′ N .These do not have a qualitative effect on the physics, so we neglect them.Using the notation In order to proceed we must place some constraints on the disorder variables.Suppose that they are bounded and not all the same.Then we may choose the smallest interval I such that µ j ∈ I for all j to be I = [−∆, ∆] with ∆ > 0. Clearly µ n = −∆.Suppose further that the least µ m > −∆ is given by µ m = −∆ + δ.Then the two-point correlation function satisfies the bounds for x 1 , x 2 < n.If either x 1 or x 2 is greater than or equal to n then G(x 1 , x 2 ) = 0. We see that the correlation length satisfies 1/ξ ∈ ln δ 2t , ln ∆ t and that ξ diverges for some t ≤ ∆ signifying a phase transition.One of the phases is a skin-effect-phase, where all the particles are localized at one of the edges (corresponding to the red line at the left edge of Fig. 2), and the other phase is a disorder-induced localized phase, where all the particles are localized at some site (corresponding to the green line at the left edge of Fig. 2).Other states show a similar behaviour.
Now we let U > 0. This has dramatic effects for arbitrarily weak interactions in the thermodynamic limit, N, M → ∞.In the ground state the particles no longer condense into the same mode, due to the repulsive interactions.Instead, the particles spread out over the different modes, so that there are competing condensates.The Bose-Einstein condensate has fragmented due to the interactions.This is reminiscent of a Bose-Glass [79,80], where the competition between disorder and interactions causes disjoint localized pockets of superfluid condensate to form.Here, in contrast, the condensates are not spatially disjoint unless the disorder is much larger than the hopping parameter.
If U is sufficiently small one of the modes will still have the largest fraction of particles, say M 0 > M/N of them, so that there are rapidities λ n = −∆ + (n − 1)U for n = 1, 2, . . .M 0 .The leading part of the two-point correlator can then be shown to satisfy where we defined X = x 1 + x 2 and D > 0.
We see that the correlation length satisfies 1/ξ ∈ ln (M0−1)U +δ 2t , ln (M0−1)U +2∆ 2t .Clearly ξ diverges at some point, marking a transition between a boundarylocalized skin-effect phase (the purple part of Fig. 2) and a site-localized phase (the green part of Fig. 2) where M 0 particles are localized at the site where the disorder potential is minimized, while the rest are distributed in condensates elsewhere in the lattice.
If U > 2∆ the energy is minimized by equally distributing the integers in Eq. (11).Suppose for simplicity that M/N = n is an integer.In that case the gap is U − 2∆, so this ground state is a Mott insulator at integer filling.The leading order behaviour of the two-point correlator satisfies again for D > 0. The correlation length satisfies 1/ξ ∈ ln (n−1)U +δ 2t , ln (n−1)U +2∆ 2t and diverges at some point, signaling a phase transition between a boundarylocalized Mott-skin effect phase and a Mott phase where the particles are evenly distributed across the lattice, as shown in the yellow and red parts of Fig. 2, respectively.Numerical calculations corroborating these results are presented in the supplemental material [69].
Liouvillian dynamics.-Nowlet us turn to the dynamics of the full Lindblad master equation.For a generic initial state we expect that the time evolution includes every eigenoperator (9), hence the two-point correlator consists of some time-dependent linear combination of all two-point correlators of the model defined by the Hamiltonian (4).By the same reasoning as for the ground state correlators the leading order time dependent asymptotics at time τ have the form G(X; τ ) = M n=1 D n (τ )g n (X) where and the D n (τ ) are obtained from the coefficients in Eq. ( 9).Since there is only dissipation it follows that the functions D n (τ ) → 0 as τ → ∞.We see that if U is non-zero there will only be a phase transition between boundary localization and site localization once the mean particle number has decreased to a finite number.If U = 0, on the other hand, there is always such a dynamic phase transition.
Discussion.-In this work, we have shown that the paradigmatic example of a non-integrable model, namely the Bose-Hubbard model, remarkably becomes exactly Bethe Ansatz solvable in the presence of dissipation terms that may be realistically implemented with cold atoms.The effective non-Hermitian Hamiltonian as well as the full non-equilibrium quantum dynamics feature a rich phase diagram that invites ample future theoretical as well as experimental work.
As mentioned above, the realization of a momentumspace lattice version of this model in the non-interacting limit in a cold atom setup was reported in Ref. 13.There the parameters were chosen so that interactions were unimportant and attractive.We expect that similar setups can be tuned to the regime of repulsive interactions described in this paper, although the real-space interactions must be chosen to be attractive in order the be repulsive in the momentum-space lattice [81].While previous realizations of the non-interacting limit have used 87 Rb atoms, 39 K may be more suitable for the interacting system, as its relevant Feschbach resonances are wider, allowing for more precise control of the interactions [82].The existing experimental implementations of the noninteracting limit can therefore be modified to include interactions, which we expect would exhibit the phase transitions described above.
The model studied in this work may be generalized to higher dimensions and more involved interaction terms.
To this end one has to devise Lindblad operators such that the particles in the system can only hop in one direction along each dimension.Then the Fock basis can be ordered so that the kinetic term is triangular, and consequently the non-Hermitian effective Hamiltonian and the Liouvillian are triangular as long as the remaining terms are polynomial in the on-site densities.If open boundary conditions are chosen it is straightforward to calculate the eigenvalues of the effective Hamiltonian, while determining the corresponding eigenvectors is technically more complicated.
An interesting conceptual aspect of this work is that the Lindbladian does not have any conserved quantities, yet it is solvable by the Bethe Ansatz, similarly to the models discussed in [65,66].There, however, the system is integrable in the absence of dissipation in contrast to the present case [47].
The topological nature of the non-interacting non-Hermitian skin effect suggests that the quantum manybody generalizations thereof studied here is also robust to generic perturbations, which is promising in the context of its realization.
Finally, we also note that an alternative interpretation of the non-Hermitian Hamiltonian Eq. ( 4) can be gleaned from [83][84][85].There the difference in the hopping to the left and to the right, respectively, are obtained from the strength of a magnetic field applied to interacting flux lines in an array of defects in a superconductor.The unidirectional Hamiltonian is obtained when the transverse field component is large.

SUPPLEMENTAL MATERIAL FOR 'LIOUVILLIAN SKIN EFFECTS AND FRAGMENTED CONDENSATES IN AN INTEGRABLE DISSIPATIVE BOSE-HUBBARD MODEL'
The supplemental material provides derivations and details in addition to the main text.In section 1 the effective Lindbladian used in the main text is derived from a more realistic model.In section 2 the structure underlying the exact solvability of the model is shown, and in section 3 more explicit expression for the eigenoperators of the Liouvillian are derived.

DERIVATION OF THE EFFECTIVE LINDBLADIAN
An effective model describing the cold atom setup shown in figure S1 is given by the Lindblad equation The physics described by this equation may equivalently be captured by a path integral in the general Keldysh formalism, see Ref. [70] for a textbook treatment.In this framework the field operators are doubled, with one set living on the forward part of the closed time contour, and the other set living on the backwards path.Let the annihilation operators acting on the lower sites be given by a j and the annihilation operators acting on the upper sites be given by b j for sites j = 1, 2, ..., N .In addition we use a superscript +(−) to indicate that a field lives on the forward(backward) part of the time contour.Taking the system to be noninteracting for simplicity it may be described by the following path integral where the Keldysh action S is given by and the Hamiltonian is given by The dissipation occurs entirely through the upper sites, so the jump operators can be taken to be L j = b j , and we let Γ j = Γ be the same for each site.In order to simplify the notation we define Then the part of the action involving the b operators is where we defined the inverse Green's function and The action S B is quadratic, so the integral over the b fields may be carried out.This procedure corresponds to moving from (a) to (b) in figure S1.We obtain the following contribution to the effective action for the a fields Now we make the assumption that M and Γ are large, implying that the timescale of the degrees of freedom that have been integrated out is small compared to the remaining system of interest.Then to leading order where we defined ℓ n = a n + e −iϕ a n+1 .Taking ϕ = π 2 yields the same jump operators as in the main text.Note that the integration led to additional nearest-neighbor hopping terms, which renormalize hopping parameter and the on-site potential in the Hamiltonian.Using a Hubbard-Stratonovich transformation it can be shown that the addition of interactions in the b fields does not impact the above calculation at the order we are working to.

DETAILS ON THE ALGEBRAIC BETHE ANSATZ
We show this by generalizing the argument in [40].The R-matrix of this Hamiltonian is the R-matrix of the isotropic Heisenberg model [71] where Take the Lax operator

THE PARTICLE DENSITIES -NUMERICAL RESULTS
In this section we present some numerical calculations of the particle densities for the unidirectional disordered Bose-Hubbard model with open boundary conditions that is considered in the main text.The results in the figures were obtained using exact diagonalization.
Figure S2 shows the disorder-averaged density ⟨GS| a † x a x |GS⟩ /⟨GS|GS⟩ of the non-interacting ground state.The disorder was chosen so that its minimum −∆ = −1 was located at position 22 and the rest was chosen to be uniformly distributed in the interval [−∆.∆].Note that the number of particles in the system was chosen to be 1, since manybody effects are not relevant in the non-interacting case.We see that for t = 0.15∆ the particle is localized at the minimum of the potential, which is located at site 22, while for t = 0.8∆ the particle is localized at the boundary.There must therefore be a transition between these phases when t is varied in the interval [0.15∆, 0.8∆], in agreement with the main text, where it was shown that this transition must occur for some t c ≤ ∆.
Figure S3 again shows the disorder-averaged density, but now in the interacting case.The disorder was chosen so that its minimum −∆ = −1 was located at position 3 and the rest was chosen to be uniformly distributed in the interval [−∆.∆].The interaction strength U was chosen to be small, U = 0.01, to ensure that the system is in the fragmented condensates phase discussed in the main text.We see that for t = 0.15∆ the largest fraction of particles are localized at 3 while for t = 0.8∆ all the particles are localized at the boundary, so there must be a phase transition when t is varied in the interval [0.15∆, 0.8∆], which is in agreement with the main text, since U is small.
Finally, figure S4 shows the disorder-averaged density in the Mott phases described in the main text.Again, the disorder was chosen so that its minimum −∆ = −1 was located at position 3 and the rest was chosen to be uniformly distributed in the interval [−∆.∆].However, now the interaction strength was chosen to be large, U = 5, to ensure that the system is in one of the Mott phases.We see that for t = 0.15∆ the particles are evenly distributed across the lattice, while for t = 0.8∆ the particles are localized at the boundary.In the latter case they are not as sharply localized as in the fragmented condensates case.Nonetheless, since the system is at unit filling this is in agreement with the main text, since there must be a transition as t is varied in the interval [0.15∆, 0.8∆].
FIG. 1.The models under consideration in this letter.(a)A triangular ladder as implemented in Refs.13 and 14, with an additional on-site potential and interactions of strength U .The upper sites have dissipation with rate Γ, while the lower ones do not.A fictitious magnetic flux ϕ leads to a relative phase in the hopping terms between the different sites.Upon integrating out the green sites we obtain (b) with dissipation acting on neighbouring sites, where the magnetic flux has been chosen to be ϕ = π 2 .(c) Shows the system described by the integrable unidirectional effective non-Hermitian (4), where the dissipation has been absorbed into the hopping parameter and the on-site potential.
FIG. S1. Figure (a) shows the triangular chain, where there is dissipation only in the upper sites.Figure (b) shows the remaining sites after the upper sites have been integrated out.An effective non-local dissipation with decay rate γ has emerged, while the on-site potentials and the hopping parameter are renormalized.
FIG. S2.Density of single particle ground states with ∆ = 1 and N = 45.The minimum of the potential was chosen to be at position 22, and the rest were chosen to be uniformly distributed in the interval [−∆, ∆] with 500 disorder iterations.The grey lines correspond to individual disorder realizations, while the red line is the average.In the left figure the hopping parameter is t = 0.15 while in the right figure the hopping parameter is t = 0.8.
FIG. S3.Density of the many-particle ground states with ∆ = 1, U = 0.1 and unit filling, N = M = 7.The minimum of the potential was chosen to be at position 22 and the rest were chosen to be uniformly distributed in the interval [−∆, ∆] with 10 disorder iterations The grey lines correspond to individual disorder realizations, while the red line is the average.In the left figure the hopping parameter is t = 0.15 while in the right figure the hopping parameter is t = 0.8.
FIG. S4.Density of the many-particle ground states with ∆ = 1, U = 5 and unit filling, N = M = 7 .The minimum of the potential was chosen to be at position 22 and the rest were chosen to be uniformly distributed in the interval [−∆, ∆] with 10 disorder iterations The grey lines correspond to individual disorder realizations, while the red line is the average.In the left figure the hopping parameter is t = 0.15 while in the right figure the hopping parameter is t = 0.8.