Bethe ansatz approach for dissipation: exact solutions of quantum many-body dynamics under loss

We develop a Bethe ansatz based approach to study dissipative systems experiencing loss. The method allows us to exactly calculate the spectra of interacting, many-body Liouvillians. We discuss how the dissipative Bethe ansatz opens the possibility of analytically calculating the dynamics of a wide range of experimentally relevant models including cold atoms subjected to one and two body losses, coupled cavity arrays with bosons escaping the cavity, and cavity quantum electrodynamics. As an example of our approach we study the relaxation properties in a boundary driven XXZ spin chain. We exactly calculate the Liouvillian gap and find different relaxation rates with a novel type of dynamical dissipative phase transition. This physically translates into the formation of a stable domain wall in the easy-axis regime despite the presence of loss. Such analytic results have previously been inaccessible for systems of this type.


Introduction
Particle loss is an important mechanism of environmental dissipation. It strongly affects the dynamics of many particle systems in the classical and quantum regimes and has an immense impact on technological applications. It is present in numerous experimental platforms including cold atoms [1,2], non-linear waveguides [3], coupled cavity arrays [4][5][6][7], THz cavities [8][9][10][11], quantum wires [12,13], condensed matter systems [14], and solid-state devices [15]. Indeed, the primary source of dissipation in these settings is a consequence of particles escaping from the system either through one-or two-body processes [2], or due to the coupling to an external electromagnetic field (e.g. [4,8,15]). The platforms underlie future quantum technologies, which will require efficient manipulation of many constituents.
Understanding the behaviour of such systems is of paramount importance, and sheds light on properties that are robust to dissipation and could therefore allow for more efficient methods of information storage and the development of novel error correction mechanisms. However, due to the exponential complexity, numerical simulations of these systems are challenging. Thus, gaining a better understanding of their properties through uncovering their analytical structure is highly desirable. Thus far, exact solutions of such systems have been limited only to the stationary states of boundary driven systems [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35] and to those with non-interacting Hamiltonians . Beyond this only certain approximate methods [58][59][60][61][62], e.g. introducing dissipation on hydrodynamical scales, are available.
In this letter we go beyond these results and develop an analytic approach to describing the dynamics of a wide class of fully interacting dissipative systems. Our approach opens a novel avenue for the analytical study of experimentally relevant many-body models experiencing loss, provided that the system's effective non-Hermitian Hamiltonian is integrable. In experimental settings, examples of systems treatable by our method can be found in cold atom quantum simulators subjected to single and two body losses [1,2,8], and driven-dissipative cavity arrays of bosons [4].
As an example of the power of our method we study the instructive and paradigmatic XXZ spin chain, often used to describe limiting cases of the aforementioned experimental setups [2], which we subject to boundary spin loss. We find that our model exhibits intriguing physical phenomena. Additionally, these types of localized loss processes recently attracted a lot of theoretical and experimental interest due to their importance for understanding transport properties and as an experimentally realistic venue for preparing interesting quantum states, see e.g. [3,12,13,16,30,34,[63][64][65][66][67][68][69][70][71][72].
Using our method we first characterize the relaxation dynamics, uncovering a dynamical dissipative phase transition [73,74], by calculating the closure of the Liouvillian gap. Next, we analytically show the presence of a novel type of dynamical dissipative phase transition that corresponds to non-analyticity in many relaxation rates beyond the leading decay mode. Physically this implies a transition in the dynamics on both short and asymptotic time scales. This should be contrasted with phase transitions in the stationary state [75][76][77][78] or the leading decay mode [73]. In our case the stationary state is always the same and a phase transition occurs in the leading decay mode and in other parts of the spectrum. Related to this, we show that a stable domain wall state is formed in the easy-axis regime. Interestingly, the domain wall formation occurs spontaneously if the system is initialized in the maximally polarized state. It arises as a consequence of boundary bound states that we solve for. Formation of domain walls in both integrable and non-integrable closed systems has also recently attracted considerable interest [79][80][81][82][83], but is currently still analytically unsolved.

Solving lossy models
We will focus on systems described by the Lindblad master equation which characterizes open quantum systems in the weak system-bath coupling limit. The dynamics of the density matrix ρ is provided by the generator L as [84,85], where H is the system Hamiltonian and L μ are the Lindblad jump operators modeling the influence of the environment on the system. The time evolution of an observable O can be computed by diagonalizing the generatorL, where λ μ are the eigenvalues and ρ μ ,ρ μ the right and left eigenvectors. The general setup that we consider comprises an integrable Hamiltonian H with a conservation law M and Lindblad operators L μ that change the eigenvalue of M by well-defined amounts m μ > 0, [M, L μ ] = −m μ L μ , inducing the loss of the quantity M in the system. For instance M can be the total particle number and L μ particle annihilation operators. In the following we will discuss the integrability requirements for applying our technique. The Liouvillian superoperator can be represented on the vector space with doubled degrees of freedom by the channel-state transformation |ψ φ| → |ψ ⊗ |φ , yieldinĝ We will show that in order to obtain eigenvalues of the Liouvillian it suffices to obtain eigenvalues E j of the non-Hermitian HamiltonianH ≡ −iH − μ L † μ L μ [86,87]. Since [H, M] = 0, we can assume that the eigenvectors |ψ j ofH are also eigenvectors of M. The generator (3) can now be decomposed into two partŝ with H ≡H ⊗ + ⊗H * and D = 2 L μ ⊗ L * μ . Since H is a sum of two operators acting on the factors in a tensor product independently, its eigenvalues read H|ψ i ⊗ |ψ j = (E i + E * j )|ψ i ⊗ |ψ j . Let us now order the eigenvectors |ψ i ⊗ |ψ j of H by the corresponding eigenvalues m i,j of M ≡ M ⊗ + ⊗ M. Due to the purely lossy dynamics the nondiagonal matrix elements of D lie strictly above the diagonal. This immediately implies thatL takes the upper triangular form in the basis |ψ i ⊗ |ψ j and that the eigenvalues of the Liouvillian λ i,j = E i + E * j coincide with those of H. Thus, provided that the non-Hermitian HamiltonianH = −iH − μ L † μ L μ is exactly solvable, we have found the full spectrum of the Liouvillian. However, the structure of Liouvillian eigenvectors corresponding to the eigenvalue λ i,j is more complicated and includes the states |ψ k ⊗ |ψ l , with m k,l m i,j as demonstrated in appendix A. Note that the dynamics describing pure gain can be treated on the same footing.
For general interacting, many-body systems, solvingH is still intractable due to the exponential complexity. However, we will now outline how Bethe ansatz techniques [88,89] may be employed to solvẽ H in a wide range of models. In many physically relevant situations the dissipative contribution, μ L † μ L μ , modifying the system's integrable Hamiltonian, H, will leaveH integrable. For instance, single particle bulk loss throughout the system in any integrable model with particle number conservation, the dissipative contribution corresponds simply to the particle number operator, which clearly implies integrability ofH. For two-level systems with conserved magnetization the L μ would correspond to on-site spin lowering operators. This is known to be a primary dissipative loss mechanism in numerous experimental setups such as optical lattices (due to interactions with the background vacuum), wave guides, solid state contacts, and coupled cavity arrays [2][3][4]15]. Other examples of dissipative mechanisms that preserve integrability ofH are provided by nearest-neighbour dissipation [90,91] and two-body loss processes [2,92]. Further details are provided in the following section which demonstrate the wide utility and applicability of our technique.
It is instructive to contrast this situation with cases where the full Liouvillian can be mapped to a non-Hermitian integrable Hamiltonian [49-51, 53, 54, 56, 57] (where the physical system's Hamiltonian is quadratic). In our case the system's Hamiltonian is interacting and the full Liouvillian does not necessarily correspond to some non-Hermitian integrable Hamiltonian. Rather here it is onlyH (and hence H in equation (4)) that is integrable.
We now provide further explicit details of the possible models where our technique can be applied before demonstrating the utility of our method we by focusing on the example of the Heisenberg XXZ chain in the presence of a spin sink at a single boundary.

Examples of dissipative quantum models solvable by Bethe ansatz
Every integrable system Hamiltonian offers at least one corresponding lossy dissipative process that renders it solvable according to our approach. Here we detail such examples.

Localized loss on two sides
A direct generalization of the loss studied in the main text is putting a loss process on the other side of the chain. This could be done for instance by studying the 1D Hubbard model, We then study four Lindblad operators which is integrable [93]. Of course, the Hubbard model is just a concrete example-any other boundary integrable model would work here, too.

Single body loss
Take any integrable U(1) symmetric Hamiltonian with symmetry N = j a † j a j , and corresponding loss which is trivially integrable.

Two-body loss
Consider the Lieb-Liniger model of cold bosons, a very experimentally relevant loss process is two-body loss [1,2,94], This may be understood as a continuous set of Lindblad operator L(x) = √ γΨ 2 (x) inducing pure loss. This leads to the triangular form the Liouvillian required. Furthermore, which is a modification of the interaction leavingH integrable. Another concrete physical example of two-body loss treatable with our method is the 1D Hubbard model with two-body recombination of fermions of spin-down and spin-up [92].

Dissipative lossy nearest-neighbour hopping
A very related loss process is the following. Consider again the XXZ spin chain described by H from the main text with In that caseL is again triangular and, which is again integrable because it corresponds up to an irrelevant constant to an addition of constant magnetic field and boundary field. It also renders the Δ complex. This is interesting from the point of quantum phases because it allows for a physically motivated reason to extend the already rich quantum phase diagram of the XXZ spin chain to complex Δ. This is similar to two-level systems coupled by dipolar interactions and subject to nonlocal dissipation, i.e. decay through optical emissions [90,91].
Other possible examples of models with collective loss may include Richardson-Gaudin models [108].

Boundary driven XXZ spin chain dynamics
We now study the example of an XXZ spin-1/2 chain with a single boundary spin sink. Our method allows us to analytically find the scaling of the Liouvillian gap and to understand the formation of domain walls in the easy-axis regime via boundary bound states which we can solve for.

Dissipative Bethe ansatz solution
The Heisenberg Hamiltonian reads where the Pauli spin− 1 2 operators are σ x,y,z , Δ is the anisotropy, and N is the number of spin sites. We study the setup with an arbitrary loss rate on the first site i.e. iH describes an XXZ spin chain Hamiltonian in the presence of an imaginary magnetic field at the boundary. The Hamiltonian has a U(1) symmetry M = j σ z j and [M, L 1 ] = −L 1 , i.e. m 1 = 1. In contrast to boundary driven spin chains [16-18, 36, 109-113] the stationary state,Lρ ∞ = 0, is not of interest in our system since it is a trivial vacuum state. However, we obtain the full spectrum of the non-Hermitian Hamiltonian equation (14), and therefore of the Liouvillian 3 by the Bethe ansatz. The Bethe equations were obtained using Sklyanin's reflection algebra [114], and equivalently the coordinate Bethe ansatz [89,[115][116][117] with imaginary boundary magnetic field 2iΓσ z 1 . The complex energies ofH corresponding to m magnons read where the momenta of the magnons, {k j }, are obtained by solving the Bethe equations The scattering matrix of two magnons takes the form From the triangular form of the Liouvillian, which couples different magnetization sectors, we can express its eigenstates in terms of Bethe states ofH by simplified Gaussian elimination (for details see appendix A). For later convenience we introduce a pair of labels (m L , m R ) re the eigenstferring toates ofL comprised of tensor product of two Bethe states with m L and m R magnons as well as tensor products of Bethe states with less q L < m L and q R < m R magnons. In what follows we will utilise our results to address two physically interesting questions.

Eigenvalue structure and the Liouvillian gap
The first problem that we consider is the Liouvillian gap, R, ofL. It corresponds to the maximum real part of the eigenvalues different from 0, which is the inverse relaxation time of the longest-lived eigenmodes. We plot it in figure 1 for different values of the anisotropy parameter Δ. The scaling of the gap with the system size N is one of the primary features of open quantum systems, governing the late time dynamics. We observe numerically that the gap corresponds to the eigenstates ofL with m L + m R = 1. By examining the single magnon m = 1 solutions of (16) on top of the steady state we find that the gap closes at different rates depending on the value of Δ. In particular, we demonstrate that for Δ 1 the longest lived excitations correspond to the solutions with lim N→∞ (k j ) = 0, while for Δ > 1 they correspond to solutions with lim N→∞ (k j ) = − log(Δ).
We find that the real part of the eigenvalue with the smallest (non-zero) real part scales as 4 , The result for Δ 1 matches the scaling of the gap for free fermions [39]. 4 See appendix B for details. Strictly speaking this is an upper bound on the gap because there could possibly be slower decay modes outside the single magnon sector, but extensive numerical calculations confirm that we find the correct value.

Boundary bound states and domain wall formation in the easy-axis regime
In the second setup we consider the case where the system is initialized in a highly excited, i.e. maximally polarized (all spins-up) state. In this case, due to the structure of the eigenstates ofL, we need only consider eigenstates with m L = m R = m. In order to study the dynamics we now focus on the most stable (maximum real part) eigenvalues in the m top-magnon sector, corresponding to spin-down excitations on top of the all spins-up state. The Bethe equations for top-magnons can be obtained from equation (16) by replacing Γ → −Γ in the sector with m magnons. Focusing on the easy-axis, Δ > 1, regime, we show that in the m top-magnon sector states with lim N→∞ (k j ) > 0, which are localized at the boundary appear and are the most stable. For these bound states the m top-magnon Bethe equations can be easily solved in the N → ∞ limit, since e ik j N → 0 . A recursive solution of gives an appealingly simple result for the leading Liouvillian eigenvalues in the m top-magnon sector, Physically this means that the first top-magnon with momentum k 1 is localized near the loss site, while the jth top-magnon is recursively bound to the (j − 1)st. Importantly, we can show that as the number of top-magnons is increased they become exponentially stabilized, i.e. lim m→∞ Re λ m = 0. In turn, this implies that exponentially large times (in m) are needed for the loss site to dissipate the state with m top-magnons. More details can be found in appendix D or alternatively [118]. The existence of these boundary bound states has intriguing physical consequences. It results in domain wall formation if the system is initialized in the maximally polarized state. Naively one might think that such a state is the most unstable, however tDMRG simulations, as shown in figure 2(a), in the Δ > 1 regime reveal that the total spin leaking out of the system increases only logarithmically with time (see figure 2(b)). This can be understood as a consequence of the exponential stability of the boundary bound states. Namely, that exponentially long times (in m) are required for the loss site to remove all the states with m down-turned spins. Moreover, in figure 2(c) we show that the dynamics of magnetization close to the spin loss site is well described by the decay rates of boundary top-magnons. On the other hand the decay of the maximally polarized state in the Δ < 1 regime is very rapid ( figure 2(a)), and the total loss of magnetization increases linearly with time. There have recently been a number of studies addressing the dynamics of domain walls in integrable [79][80][81]83] and nonintegrable [82] systems without dissipation. While the ballistic expansion in the Δ < 1 regime is well understood, the domain wall freezing was analytically unresolved.
The existence of boundary bound states also has profound consequences on the spectral properties ofL. It results in a dissipative phase transition (shown in figure 3). In contrast to standard dissipative phase transitions [75], the stationary state remains the same (all spins-down). The phase transition rather happens in all sectors of the relaxation spectrum. At leading order in Γ as N → ∞ the transition occurs at Δ = 1. This is similar to dynamical dissipative phase transitions [73], but the discontinuous eigenvalues that are relevant for the dynamics are not only the Liouvillian gap. This is reflected in the fact that already the short time dynamics well inside the easy-plane and easy-axis regimes are qualitatively different (see figure 2(a)). The discontinuity is shown in figure 3 where, by taking small Γ, we can see non-analyticity in eigenvalues in three different top-magnon sectors at Δ = 1. This is characteristic of all sectors. The non-analyticity shown corresponds to the non-existence of boundary bound state solutions for Δ < 1 for all values of Γ. More specifically, we prove the existence of {k j } such that lim Δ→1 dk j dΔ → ∞ for large N and small Γ, which implies the divergence in the corresponding eigenvalues. See appendix C for further details.

Conclusion
We have developed a framework for diagonalizing quantum Liouvillians with integrable system Hamiltonians and dissipative loss. We demonstrate the utility of our method in an example of the Heisenberg XXZ spin chain with boundary loss. The method allows us to directly identify phase transitions in the Liouvillian spectrum and calculate the Liouvillian gap. This led us to observe two intriguing physical phenomena, namely domain wall formation, and a dissipative phase transition, which we link to the existence of boundary bound top-magnons. Such remarkable phenomena could occur in other models with localized loss, e.g. 1D Hubbard and interacting bosons in 1D [2], which can be studied analytically with our method.
A number of questions remain open. The first natural extension of our results is directly calculating the full eigenstates of the quantum Liouvillian. We also envisage using the thermodynamic Bethe ansatz [119, 120] to explore the decay of states with a finite density of excitations, and the connection with boundary states [121,122] and strong edge modes [123] in closed systems. Additionally, the Liouvillian spectrum exhibits a multi-band structure in the Ising limit, Δ → ∞, (see figure 4 of appendices), the ramifications of which remain unexplored.
More generally our method can be applied to a number of physically relevant systems that are quantitatively very different from the example studied here. We discuss these further in section 3. Here the interest is two-fold. On one hand, judging by our example, such systems hide a plethora of interesting physical phenomena, which are yet to be uncovered. On the other hand they describe realistic experimental setups and therefore provide an indispensable tool for understanding future experiments.

Note Added
While nearing the completion of this manuscript a related preprint appeared [92] studying exact solutions in the Hubbard model with two-body loss.

Appendices
In these appendices we provide the discussion of technical details that were omitted from the main text. In the first section we construct the eigenstates of the Liouvillian. In the second and third sections we provide details on the calculation of the Liouvillian gap, and the results related to the boundary magnons, respectively. We also plot in figure 4 the full Liouvillian spectrum for different values of anisotropy Δ which shows the formation of intriguing band structure in the Ising, Δ → ∞ limit.

Appendix A. Eigenstates of the boundary loss XXZ Liouvillian
In this section we will construct the right eigenstates of the Liouvillian In the case of the XXZ chain with a single loss at the first site. First of all, we will remind the reader of the basic structure of Bethe eigenstates, which will then serve to construct the eigenstates of the Liouvillian. The eigenstate of the HamiltonianH pertaining to the energy E φ a η reads  where the summation is performed over all permutations and negations of {k j }, and ε P changes sign with each such mutation. We then use the triangular form for the Liouvillian, i.e. that where M a is the subspace with magnetization a, to make the ansatz that the eigenstate of L with eigenvalue Substituting this gives the recurrence relation for μ = 1, . . . , min{a, b}, where we defined In general, this recursion relation is significantly more complex than computing the Bethe states ofH. It does however provide an insight into the structure of the eigenstates. We also note that one of the main powers of Bethe ansatz lies in the thermodynamics and that efficient calculations might still be possible in such a limit.
Appendix B. Thermodynamic limit of the leading decay rates B.1. Easy plane, Δ 1, regime Here we will study the solutions of the Bethe equations that have purely real momenta in the thermodynamic limit. These solutions can be physically understood as free magnons that live in the bulk of the system and only experience the effects of other magnons and the boundary in sub-leading order 1/N. To do this we start with the logarithmic form of the Bethe equations, In order to simplify discussion we focus on the single magnon case, though the solutions for m magnons are also straightforward in the above discussed limit, We denote by and expand the momenta k as the power series in 1/N, k = k (0) + 1/Nk (1) + · · · , which we truncate at the order O(1/N 3 ). It is important to distinguish the cases when the integer I 1 is finite and when it is of the order of the system size N or close to N. Since the leading decay mode corresponds to the latter case, we make the transformation I 1 → N − I 1 and focus on finite I 1 .
Expanding (B.2) is straightforward, as is solving it order by order. We arrive at, which we may solve for a 1 . Using (B.5) with k i = 0, k j = i log 1 Δ − a 1 Δ −2N gives the result in the main text.
Strictly speaking these results for both the easy axis and easy plane regime represent an upper bound on the gap, although extensive numerical evidence also shows that they are also a lower bound. We hope to prove that this is the gap by extending techniques for the isolated XXZ spin chain [130].
Appendix C. Calculation of the phase transition in the highly excited eigenstates We will now study the one top-magnon (spin-down in a background of all spins up) sector. The single top-magnon cases correspond to setting Ω(a) = log (aΔ−1)(−1+a(Δ+2iΓ)) (a−Δ)(a−2iΓ−Δ) in (B.2). The corresponding energies ofH are, while the eigenvalues ofL read, We now look at the case when I 1 is finite. Numerically we observe that this corresponds to the leading decay rates in the single top-magnon sector for Δ < 1 in the small Γ limit. Performing the same expansion as before, relabeling I 1 → p, we obtain (now it is sufficient to look at only up to order 1/N 2 ), We take the derivative w.r.t. to Δ of (C.2), Using (C.3) we obtain, which diverges as Δ → 1 at leading order in Γ, signalling the phase transition in these highly excited states.

Appendix D. Boundary bound modes and stability of the domain wall in the easy-axis regime
In the easy-axis regime, Δ > 1, an infinite number of solutions that have non-zero imaginary part in the thermodynamic limit appear, lim N→∞ (k j ) = 0. In contrast, in the easy-plane regime we observe that only a finite number of such solutions can appear at finite Γ = 0. We study the solutions with lim N→∞ (k j ) > 0. These correspond to top-magnons localized at the boundary loss site, which we refer to as boundary bound modes. More specifically, we may solve the m top-magnon Bethe equations in the N → ∞ limit by observing that e ik j N → 0. Focusing on the top-magnon boundary bound modes, we arrive at the following simple form of the Bethe equations in the N → ∞ limit, 1 + e ik j (−Δ − 2iΓ) i =j −2Δe ik j + e i(k j −k i ) + 1 e i(k j +k i ) − 2Δe ik j + 1 = 0. Physically, this means that the k 1 top-magnon is localized at the loss site, whereas the jth top-magnon is bound to the (j − 1)st one. This characterises a domain wall state.
We will now show that the real part of the eigenvalues decay exponentially with top-magnon number m. First, we decompose the equations for energies (C.1) as Regrouping the terms, we simply get Solving for the stationary value of recursion, z = exp(ik 2j+1 ) = exp(ik 2j−1 ), we obtain two real solutions, z 1 = Δ − −1 + Δ 2 , z 2 = Δ + −1 + Δ 2 , ( D . 6 ) with the stable point being z 1 . We numerically observe that this fixed point is converged to for any initial value of Δ > 1 and Γ. The decay of the most stable eigenvalue λ m = 2 Re(E m ) is thus exponential in top-magnon number m, demonstrating the stability of the domain wall.