Dissipative cooling towards phantom Bethe states in boundary driven XXZ spin chain

A dissipative method that allows to access family of phantom Bethe-states (PBS) of boundary driven XXZ spin chains, is introduced. The method consists in coupling the ends of the open spin chain to suitable dissipative magnetic baths to force the edge spins to satisfy specific boundary conditions necessary for the PBS existence. Cumulative monotonous depopulation of the non-chiral components of the density matrix with growing dissipation amplitude is analogous to the depopulation of high-energy states in response to thermal cooling. Compared to generic states, PBS have strong chirality, nontrivial topology and carry high spin currents.

Introduction. Dissipation needs not to be always destructive for quantum protocols but it can represent a resource for manipulating quantum systems. Dissipation alone [1] or in combination with coherent dynamics [2][3][4][5][6][7][8][9][10], indeed, can be used to create quantum nonequilibrium stationary states (NESS) which are attractors of the dynamics and therefore are stable even in the presence of noise.
Most protocols, however, require a tailored set of operations used in pumping cycles to target each specific state [5,11]. If the protocols are implemented by stationary control fields, they usually require sophisticated dissipations that make the NESS targeting more complicated [12]. This is due to the fact that the targeted NESS must be an eigenstate of the coherent part of the dynamics and a dark state for all jump operators in the dissipator [13,14].
Here we demonstrate how to generate a remarkable family of NESS containing an arbitrary number of qubits, employing simple boundary-localized dissipation, and manipulating just one parameter. These are the phantom Bethe states (PBS), i.e. eigenstates of integrable XXZ spin chains on special parameter manifolds [15][16][17], possessing unusual chiral and topological features.
The phenomenon is based on a subtle mechanism that makes (within the quantum Zeno limit [18][19][20][21]), a highly selective 'phantom' invariant subspace the basin of attraction for the density matrix. Consequently, the NESS responds in a singular resonant way to an increase of the dissipation strength in the vicinity of "phantom" manifolds, restricting the density matrix to states with chirality of the same sign and thus rendering the NESS chiral. The resonances become sharper as the dissipation strength is increased and their number grows with the number of spins involved.
Dynamically, the "freezing out" of the non-chiral components of the density matrix with growing dissipation amplitude is analogous to the depopulation of highenergy states in response to thermal cooling.
The "dissipative cooling" method consists in coupling the ends of the open chain to dissipative baths of polarization constraining the first and the last spin to relax to predefined pure qubit states that satisfy boundary conditions necessary for the PBS existence. We show that by changing the control parameter, i.e. the misfit azimuthal angle of the dissipatively-targeted boundary polarizations, it is possible to thread "phantom" manifolds, passing from one chiral NESS to another one with different topology (see Figs. 2). To the best of our knowledge, a quantum protocol that allows to target the whole PBS family is presently lacking.
The simplest states belonging to the phantom family are the spin-helix states (SHS) (see top panel of Fig. 4 for an example). Recently, SHS were created and used as a sensitive tool to measure the anisotropy in experiments on XXZ chains implemented via ultracold atoms [22,23]. While in [22,23], the lifetime of the SHS is restricted by finite size effects, the stability of dissipatively created SHS (and other PBS) is guaranteed as long as the boundary dissipation strength is kept sufficiently strong.
Note that while SHS are pure states, all other PBS are mixed states corresponding, in terms of the mean variation of the magnetization along the chain, to helices with a variable radius (see middle panel of Fig. 4 for an example). Quite interestingly, we find that in comparison to generic pure and mixed states of the open XXZ arXiv:2206.06771v3 [cond-mat.stat-mech] 19 Oct 2022 chain, the PBS have strong chirality and carry relatively high spin currents (see Fig.3), a feature that can be of potential interest for future applications. Model setting and near Zeno limit relaxation. Our dissipative setup is schematically depicted in Fig. 1. An XXZ Heisenberg spin 1 2 chain of N +2 sites, numbered as 0, 1, . . . N +1, is coupled to Lindbladian baths of polarizations at the ends ( [18,24], see details in [25]. We assume dissipative targeting of generic pure qubit states ρ L,R = (I + n L,R · σ) /2 at sites 0 and N + 1, where n L , n R are unit vectors of polarization. The strength of the dissipation Γ is measured by the inverse time needed for edge spins to relax, e.g. ρ L/R (t) = I 2 + 1 2 σ 1 n L/R + O(e −Γt ), if the coherent part ( XXZ spin chain Hamiltonian) is neglected. If n L = n R , a nonequilibrium gradient is created and a spin current j z can flow, typically obeying the Fourier law j z = O(1/N ).
In (8) we fix a resonance condition under which steady currents can be increased up to maximally possible values, j z = O(1). Using the criterion of [26] one finds that the NESS is always unique.
It was shown in [27] that close to quantum Zeno limit Γ → ∞ the relaxation to NESS undergoes a threestage process, each stage occurring at different timescale. On the shortest timescale: t ≤ O(1/Γ), the boundary spins relax towards their targeted states. From this point on, the density matrix is approximately given by describes the time evolution of the the internal part of the system, which becomes approximately coherent again [2] and is governed by the dissipation projected Hamiltonian h D given by [27] h D = N −1 n=1 σ n ·Ĵ σ n+1 +Ĵ n L · σ 1 +Ĵ n R · σ N , (1) withĴ = diag(1, 1, ∆). Note that the dissipationprojected Hamiltonian (1) depends on polarizations of the baths through its boundary fields.
On the intermediate timescale: , the reduced density matrix R(t) acquires the approximate form where τ = t Γ and |α are h D eigenstates, with positive coefficients P α significantly changing on long time scales τ = O(1). Finally, on the long time scale O(1) t ≤ O(Γ) the coefficients P α (τ ) → P α (∞) ≡ P α relax to their stationary NESS values, by following an effective Markov process with rates w αβ ≡ w α→β given by [27] where g L , g R are operators acting on the first and the last spin respectively, given by where θ, ϕ are spherical coordinates of a unit vector. The final state, the NESS, in the Zeno limit has the form where P α is time-independent solution of (3), and |α are eigenstates of h D (1). Phantom Bethe states. Phantom Bethe states are eigenstates of the Hamiltonian (1) that have exceptional chirality and correspond to special manifolds of the boundary fields in (1)Ĵ n L (θ L , ϕ L ),Ĵ n R (θ R , ϕ R ), with θ, ϕ given by [15,16] where M + = 0, 1, . . . N + 1 [28] and γ parametrizing the XXZ model anisotropy ∆ ∆ = cos γ.
The substitution M + → M − ≡ N +1−M + in (8) leads to the physical setup with the opposite boundary gradient, ϕ R −ϕ L → −(ϕ R −ϕ L ), and consequently flips the steady current j z → −j z . The respective NESSs are related via the left-right reflection ϕ L ↔ ϕ R and subsequent rotation in XY -plane, see [25]. Using this property we restrict the M + range to the M + = 0, 1, . . . [ N +1 2 ]. For fixed M ≡ M + in this range all eigenstates |α of h D (1) which determine Zeno NESS (7) split into two chiral families, {|α + }, {|α − } characterized by the socalled phantom Bethe roots [15]. All eigenstates from each family share common chiral properties [16]. In- . According to the above, the reduced density matrix R in (2) on "phantom" manifolds (8) splits as The sum (10) contains projectors on states with opposite chiralities and is generically approximately neutral. The time evolution obeys the effective Markov process (3) i.e. depends on rates w βδ exclusively. Analyzing the rates [25] we find a remarkable property: all the rates vanish, while generic w β→α (1) + remain finite. Thus, the subfamily {|α (1) + } becomes an adsorbing basin of the Markov process (3) [29] resulting in depopulation of all other states with time, and leading to the NESS of the form All eigenstates |α (12) are phantom Bethe eigenstates [15] of the same chirality and the chirality gets more pronounced for small M . For M = 0, the sum in (12) contains just one term, a projector on the spin-helix state (SHS) (13) visualized in Fig. 1, characterized by a large current of magnetization A possibility to target spin-helix state (13) dissipatively was also noted in previous studies [14,30,31]. For M > 0 the ideal helix (13) gets blurred but the NESS (12) remains chiral. For M = 1 the current of magnetization α + |ĵ z |α + averaged over states |α + is of [16].
From the above result we predict the existence of chiral Zeno NESS with unusually high magnetization current at phantom Bethe manifolds (8). Numerical results. To check our predictions we study NESS magnetization current for fixed anisotropy ∆ and varying boundary gradient, for systems of size 4 ≤ N ≤ 30. We use exact numerical diagonalization for small chains [32] and Matrix product ansatz for NESS in the Zeno limit [33,34] for large chains. Already for N = 5 we find all j z peak positions at predicted points (8) and their vicinity, see Figs. 2. Namely, all points (8) with M ± = 0, 1, 2 correspond to peaks of various amplitudes (j z = 0 for M + = (N + 1)/2 = 3, since it corresponds to zero misfit angle Φ = ±(N + 1 − 2M + )γ = 0 and hence to absence of the boundary gradient n L = n R ). For larger chains (Fig. 3) the agreement and the phenomenon becomes striking: most j z peaks appear as sharp resonances centered at phantom Bethe states manifolds (8), on top of a background with j z = O(1/N ). Notice that the empty black and red circles in (Fig. 3) show the dependence of the j z peaks on M (top horizontalaxis).
The role of the parameter M determining the NESS rank in (12) via r(M ) = N +1 M deserves special discussion. The highest and sharpest of all j z peaks always corresponds to M ± = 0, see Figs. 2, 3, i.e. pure Zeno NESS, with magnetization winding in a perfect helix around the z axis, see (13) and Fig. 4. For M > 0, perfect chirality is lost: the basis of "phantom" manifold for M = 1 consists of 2 disjointed helix pieces separated by a kink at some position n; phantom Bethe states are linear combinations of basis states for all possible kink positions [16]. Resulting NESS magnetization profile is a distorted helix with variable radius, see mid-panel of Fig. 4. Basis states for arbitrary M contain M kinks [16], introducing more helix imperfections. For large M NESS becomes a It is also worth to note that, due to larger number of contributing states, peaks with larger M are less sharp and have smaller amplitude, but, on the other side, the NESS gets more stable with respect to perturbations (boundary misfits, lowering dissipation, etc.). This effect can be seen by comparing the behavior of peaks at increasing θ R − θ L misfit (upper panel of Fig. 2) and for different Γ values (the bottom panel).
In particular, as Γ decreases, the near Zeno limit description of an effectively coherent evolution (2) becomes invalid, and the peaks gradually smear out, the sharper peaks first.
In contrast, for parameters chosen away from "phantom" manifolds, variations of Γ do not lead to any drastic effects (data not shown), especially for large enough N . The reason is, beyond certain characteristic value with dissipatively created boundary gradient are due to the existence of special manifolds (8) where phantom Bethe roots solutions of the Bethe Ansatz equations exist. The mechanism by which the PBS can be physically accessed was identified to be strong dissipation, which drives the system towards a chiral invariant subspace. The gradual depopulation of the non-chiral states resembles the depopulation of highly energetic states in a quantum system coupled to a cold thermal bath. The amplitude of dissipation plays the role of inverse temperature and a chiral subset of states plays the role of a subset of low-energy states close to the ground state of the system. We showed that by varying the system size, the bulk anisotropy and the boundary dissipative driving one can manipulate a number of peaks of the magnetization current, distance between the peaks in the parameter space, and their magnitudes. All the resulting stationary states are easily distinguishable by the value of the carried spin current, and by their nontrivial topology [35].
We expect the dissipative cooling approach considered in this paper to be effective also for other open quantum many body systems that are integrable via the offdiagonal Bethe ansatz method.
VP acknowledges financial support from the European Research Council through the advanced grant No. 694544-OMNES, and from the Deutsche Forschungs-gemeinschaft through the DFG projects KL 645/20-1, KL 645/20-2. VP thanks the Department of Physics "E.R.Caianiello" for hospitality and for short visits partial supports (FARB 2018 -2019) during which the work was completed.

Dissipative cooling towards phantom Bethe states in boundary driven XXZ spin chain by Vladislav Popkov and Mario Salerno
This Supplemental Material contains four sections organized as follows. In S-I we collect some useful standard definitions. In S-II we derive symmetry properties of the Lindblad equation. In S-III we prove the existence of invariant subspaces at Bethe manifolds (8). In S-IV we analyze the effective Markov process (3).

S-I. LINDBLAD MASTER EQUATION AND SYSTEM EVOLUTION NEAR ZENO LIMIT
A density matrix of a generic quantum system with dissipation satisfies, under standard assumptions [18,24], a Lindblad Master equation (LME), where H † s = H s is a coherent part, D[ρ] is the dissipator, and L γ are Lindblad operators, describing the interaction with the environment. In our model H s is taken as the Heisenberg XXZ Hamiltonian with z anisotropy J z /J x = J z /J y = ∆: with n L , n R unit vectors of polarization, parametrized by spherical coordinates angles 0 ≤ θ L , θ R ≤ π and 0 ≤ ϕ L , ϕ R < 2π. The explicit form of L L,R can be obtained by observing that the targeting of a pure qubit state |ψ ψ| at a single site, is achieved by means of operators L of the form L = |ψ ψ ⊥ |, with ψ ⊥ |ψ = 0, so that the dissipator D L [ρ] automatically satisfies D L [|ψ ψ|] = 0. Observing that the normalization of the states |ψ , |ψ ⊥ entails tr(L † L) = 1 and that the linear map D L [ρ] has four eigenmodes with eigenvalues 0, −1, −1, −2, one can write a formal solution, valid for large times t, of the Eq. (S1) with H s → 0 and with a single operator L = |ψ ψ ⊥ |, as: Note that the a similar solution can be written also for Γ ||H s ||, i.e. when the coherent evolution can be considered as a small perturbation. From the above it follows that the targeting of the pure states at the edges can be achieved by a dissipator D in (S1) with Lindblad operators SM-2 acting at sites 0, N + 1, respectively. It can be proved that the steady state, ρ N ESS (Γ) = lim t→∞ ρ(Γ, t), of (S1) is unique for any given choice of the targeted polarizations and in the Zeno limit we have We also recall (see [27] for details) that at time scale t = O(1) and large Γ, the density matrix ρ(t) acquires the approximate form (S7) with R ≡ R(tΓ) = R(τ ), evolving in time according to where D {{ is an effective dissipator and h D = h D † is the so-called dissipation-projected Hamiltonian [2] given in Eq.(1) of the main text. From Eqs. (S7), (S8) it follows that [h D , R] = 0 and R is given by with |α denoting h D eigenvectors. Moreover, the coefficients P α (τ ) evolve adiabatically under the perturbative effective dissipator D ef f in (S8). The adiabatic evolution is given by the Markov process discussed in sec. S-IV.
Remarkably, it was shown in [36] that the Zeno NESS (S7) is robust against possible asymmetries of the dissipation between left and right boundaries, i.e. the asymmetric rescaling of the Lindblad operators L 1 → δ 1 L 1 , L 2 → δ 2 L 2 in (S1), with any finite nonzero δ 1 , δ 2 values, yields the same limiting result (S7).

S-II. SYMMETRY PROPERTIES OF THE LINDBLAD EQUATION
In the following we investigate symmetry properties of the NESS of the Lindblad equation (S1) with targeted polarizations lying in XY -plane, θ L = θ R = π/2, ϕ L = 0, and arbitrary ϕ R . For this we introduce the operator U ϕ performing a simultaneous rotation of all spins around the z axis and the parity operator P, respectively defined as Notice that the bulk Hamiltonian H s commute with both operators U ϕ and P and hence also with their product U ϕ P ≡ U ϕ , i.e. [H, U ϕ ] = 0. On the other hand the Lindblad operators L 1 , L 2 from (S5), (S6), under the action of U ϕ R transform as: i.e. one gets the same physical setup but with the opposite boundary gradient at the right boundary. Indeed, it is straightforward to verify that LME (S1) under the U transformation is mapped onto the same LME but with opposite boundary gradient, ϕ R → −ϕ R , with the corresponding ρ matrix transformed as The uniqueness of the NESS for any choice of the boundary parameters, implies that the above transformation is a map between NESS of opposite boundary gradients. On the other hand, taking into account that M + + M − = N + 1 we have, from the ϕ R angles parametrized by the integer M + in Eq. (8), that: thus the transformation (S10) allows to get the M − −NESS from the M + one as: Notably, the steady state current under the transformation ϕ R → −ϕ R changes its sign: j z (M − ) = − j z (M + ) , in agreement with what observed in Fig. 2 and Fig. 3 of the main text.

S-III. INVARIANT SUBSPACE WM OF hD AT PHANTOM BETHE MANIFOLDS
To investigate the invariant phantom Bethe manifolds it is convenient to introduce the following parametrization for an arbitrary pure qubit state where f, F are real numbers. Notice that the state |f, F describes a spin 1/2 pointing in the direction ( σ x , σ y , σ z ) = (sin θ cos F, sin θ sin F, cos θ) with tan(θ/2) = e f , thus −∞ < f < ∞, parametrizes the polar angle while 0 ≤ F < 2π parametrizes the azimuthal angle (or phase factor) of the usual polar coordinate system. The state |F instead, describes a fully polarized spin 1/2 lying in XY −plane. Here for simplicity we consider the condition (8) at θ = π/2 (dissipation baths polarizations in the XY plane) to prove the following Theorem.
Theorem The set of states Remark.-In (S15) n 1 = 1 and n M = N + 1 denote virtual kinks at outer links of the chain (0, 1) and (N, N + 1): n 1 = 1 and n 1 > 1 correspond to first qubit in (S15) being | − ϕ and |ϕ respectively. Likewise, n M = N + 1 and n M < N + 1 correspond to the last qubit of the form |(N − 2M + 2)ϕ and |(N − 2M )ϕ respectively. Introducing the notation the operatorh D = h D − (N − 1)∆I can be written as sum of local terms where h = σ x ⊗ σ x + σ y ⊗ σ y + ∆(σ z ⊗ σ z − I ⊗ I) is the local energy-density of the XXZ spin chain.
A key point in the proof is the divergence condition: Using (S20) one obtains the following useful relations: Due to (S20),(S22), (S23), any factorized state of the form: under the action of the operator N −1 n=1 h n,n+1 will transform into a sum of terms Ψ α of the same type of (S25), i.e. Ψ α = ⊗ N k=1 |F k,α , with F k+1,α − F k,α = ±γ, plus two extra terms :

SM-4
with A 1 = κ sign((F 1 − F 2 )/γ) and A N = κ sign((F N − F N −1 )/γ). All basis vector of W M are of the form (S17) with F 1 = ±γ , F N = F c (M ) ± γ. By acting withh D on the basis elements n 1 n 2 . . . n M , we generate other terms of the same basis m 1 m 2 . . . m M ≡ |F 1 , . . . F N and extra terms with "impurities" at site 1 and site N : Since A 1 = sign((F 1 − F 2 )/γ)κ and F 1 = ±γ, we have that the following four A 1 |F ⊥ 1 , F 2 . . . terms may arise: Terms of type (S26) cancel with respective terms from h l ⊗ I ⊗ N −1 |F 1 , F 2 . . . , using the relations and similarly, one get cancelation of the terms of type (S27) by means of the relations Proceedings in the same manner for the arising four terms of type A N | . . . F N −1 , F ⊥ N : one finds that terms (S30),(S31) cancel with respective terms respectively (note that relations (S33), (S34) differ from (S28), (S29) simply by the increase of the angle by Φ c (M )). From this we conclude that the action of h D on any basis vector of W M produces states that still belong to W M , thus proving the Theorem (QED). We close this section with the following remarks. Remark 1.-It can be shown that the repeated action of h D on any single basis vector generates the whole subspace W M . Moreover, it follows that all eigenvectors |α + ∈ W M of h D have nonzero overlaps with all vectors n 1 . . . n M of W M basis. The total number of states in W M in (S15) can be calculated combinatorially as: Remark 2.-The subspace W M is a part of a larger invariant subspace G + M described in [16]. G + M contains all vectors of W M plus additional linearly independent vectors, all with qualitatively the same chirality features. The total number of vectors in G + M is see [16]. The effective Markov evolution for the probabilities P α of the occupation of state α at time τ : is related to a reduced density matrix dynamics R(τ ), near the quantum Zeno limit, given by where dimH = 2 N is the dimension of the Hilbert space, and |α are eigenvectors of h D . The rates w αβ ≡ w α→β are given by [27]: where g L , g R are some locally acting operators [27].   (2) + are also orthogonal. We will be interested in cumulative probability current towards set of states in W M from other blocks, given by   Eqs(S40), (S41) show that one has strictly positive probability current (S39) at all times, leading to gradual increase of cumulative population of states with the same chirality P (α