Scale-free non-Hermitian skin effect in a boundary-dissipated spin chain

We study the open XXZ spin chain with a PT-symmetric non-Hermitian boundary field. We find an interaction-induced scale-free non-Hermitian skin effect by using the coordinate Bethe ansatz. The steady state and the ground state in the PT broken phase are constructed, and the formulas of their eigen-energies in the thermodynamic limit are obtained. The differences between the many-body scale-free states and the boundary string states are explored, and the transition between the two at isotropic point is investigated. We also discuss an experimental scheme to verify our results.

In this article, we investigate a non-Hermitian open XXZ chain using the coordinate Bethe Ansatz.The chain is subjected to opposite imaginary magnetic field on two ends, pointing to a prescribed direction called z.Here, the non-Hermiticity naturally stems from the ubiquitous coupling with the environment.In the previous literature, exact solutions for spin chains with arbitrary boundary fields (diagonal or off-diagonal) have been studied by systematic methods [10,11,27,28,[31][32][33]36].However, the analysis of the ground state properties has been mainly conducted for the cases of Hermitian boundary fields.Therefore, a general method to identify the steady-state (with the largest imaginary part of energy) and the ground state (with the lowest real part of energy) in the presence of non-Hermiticity is still lacking hitherto.A special case is that the strength of the boundary field takes a specific value depending on the anisotropic interaction strength between adjacent spins.The Hamiltonian then respects the q−deformed SU (2) symmetry [37] with |q| = 1, and serves as a representation of the Temperley-Lieb algebra [38,39].The spectrum of the model is purely real, though the Hamiltonian is non-Hermitian.Furthermore, when q is a root of unity (i.e., q n = 1 for a certain integer n) for some values of the boundary field, the representation of the symmetry group enjoys richer structures, such that an exact duality exists between the spin model and free-end quantum Potts model [1].The duality leads to the same conformal field theory (CFT) structure of two models, i.e. the non-unitary CFT [40,41].Recently, the entanglement criterion has been developed for the ground state of the non-Hermitian XXZ chain to show the agreement with the non-unitary CFT [42,43].
A significant consequence of q being a root of unity is that the spin Hamiltonian develops Jordan blocks, which feature exceptional points.The number of Jordan blocks for given q and size N has been counted [44], followed by the constructions of the corresponding generalized eigenstates [45].Given the existence of many-body exceptional points, it is natural to identify different phases around them.Our article exhausts the parameter space of boundary imaginary field and anisotropic interaction.For a small boundary field, the spectrum remains real, and the Bethe roots of the ground state only shift slightly compared with the Hermitian case.The spectrum becomes complex, however, when the boundary field exceeds the q−deformed SU(2) symmetric value.We show that, despite the breaking of the quantum group invariance, the model possesses a novel behaviour of scale-free localization.We figure out the structure of the steady-state and the ground state as the combination of a boundary Bethe root and a set of continuous Bethe roots in the thermodynamic limit.The continuous Bethe roots all have an imaginary part inversely proportional to the system size N , corresponding to a small imaginary wavevector (or momentum) κ ∼ α/N .In the singleparticle context, when the localization length of wavefunction is proportional to the system size N , the density |ψ N (x)| 2 ∼ exp(2αx/N ) is invariant under re-scaling transformation with a factor s: |ψ N (x)| 2 = |ψ sN (s x)| 2 , and therefore called scale-free non-Hermitian skin effect (NHSE) or critical NHSE [46][47][48].The original NHSE means the exponential localization of most eigenstates near the boundary, with localization lengthes independent of the system size [49][50][51][52][53][54][55]; the scale-free NHSE is therefore a weaker version of NHSE.Scale-free NHSE has also been found in Hermitian systems with non-Hermitian boundary field, though the mechanism is different [56].In the present work, the imaginary part of wavevector is attributed to the scattering between the boundary mode and magnons traveling in the bulk, and these Bethe roots contribute a non-negligible imaginary part to the energy.Thus, unlike previous works, our scale-free behaviour has a many-body origin.More precisely, it originates from the interplay between boundary dissipation and many-body interactions.On the other hand, compared to recent studies reporting on the many-body NHSE in exactly solvable models with non-Hermitian hopping terms in the bulk [57,58], in this article we highlight the non-perturbative effect of the boundary non-Hermiticity on the Hermitian spin chain.We derive an integral equation, dubbed the imaginary Fredholm equation, 1 to solve the scale-free localization length in the thermodynamic limit.We then give an exact formula for the imaginary part of the steady-state energy, which are then compared to finite-size numerical results.We also explain how to measure these physical quantities in cold-atom experiments.
Before proceeding, we compare our results to earlier studies on boundary-driven spin chains as open quantum systems.The evolution of those open quantum systems is generated by the Linbladian operator, composed of an integrable Hamiltonian and quantum jump operators on the boundaries.A typical example relevant to our work is [18] Here, H eff is the non-Hermitian Hamiltonian we shall focus on below (see Eq. ( 1)).Although the Lindbladian breaks integrability, the density matrix of non-equilibrium steady-state (NESS) has been established by the matrix product operator (MPO) Ansatz exactly.Furthermore, it has been found that the local matrix of MPO Ansatz is indeed the infinite-dimensional solution of Yang-Baxter relations, and thus exterior integrability emerges in the NESS [59,60].However, the dynamics towards NESS is unknown yet.Our work about the non-Hermitian effective Hamiltonian is complementary to the NESS solution because H eff governs the time evolution of the open quantum system under post-selection, which is relevant to numerous experiments [61].Our solution is enabled by the Yang-Baxter integrability of the model.Another related system is the XXZ model with only one jump operator L 1 on the left boundary [24].Since the dissipator is purely lossy, the Lindbladian becomes upper-triangular under an appropriate basis choice, so that the Liouvillian spectrum can be completely determined by the effective non-Hermitian Hamiltonian.The Hamiltonian has scale-free eigenstates even in the single-magnon sector, but P T symmetry is absent due to that the dissipator occurs only on one of the two ends.By contrast, our Hamiltonian preserves P T symmetry, and in single-magnon sector there are only Bloch-wave modes and exponentially localized states.Scale-free modes originate from many-body interactions in our model.
The rest part is organized as follows.In the next section, we introduce the model Hamiltonian, its general Bethe equations, and the phase diagram.In Sections 3.1 and 3.2, we consider the single-magnon and two-magnon state as a warm-up.We then generalize the results to the 1 The (second kind) Fredholm integral equation refers to an inhomogeneous linear integral equation of the form φ(x) + b a d y K(x, y)φ( y) = f (x), and solving for φ(x).In another exactly sovable model: Lieb-Linger model of repulsive bosons, the corresponding integral equation for the momentum distribution on the ground state is specifically named the Lieb equation.many-body cases to obtain the steady-state with scale-free NHSE in Section 3.3.Section 3.4 is devoted to another type of steady-state solution, the boundary string states, which emerges for the highly anisotropic case.In Section 4, we apply the Ansatz of scale-free solutions to the ground state for different parameters.A possible experimental setup for the non-Hermitian model is discussed in Section 5. We give some concluding remarks in Section 6 .

Non-Hermitian XXZ model and the phase diagram
The Hamiltonian reads: where S α = 1 2 σ α (α = x, y, z) is the spin-1/2 operator; the anisotropic interaction strength ∆ and boundary field strength g are purely real, with g > 0. The unit circle ∆ 2 + g 2 = 1 has been the focus in previous literature [1,[37][38][39]42] because it enjoys the q-deformed SU (2) symmetry which greatly simplifies the problem.Here, we explore the entire (∆, g) parameter space beyond this circle.We find new phenomena, including the interaction-induced scale-free non-Hermitian skin effect, that exist outside the unit circle.
The model respects the P T symmetry with T i T = −i and PS α j P = S α N − j , and therefore the eigenvalues are either real or form complex conjugate pairs.When the whole spectrum is purely real, the model is said to be in the P T exact (or P T -symmetric) phase, otherwise it is in the P T broken phase [62,63].The steady-state, which has the largest imaginary part of eigen-energy in the P T broken phase, is of great importance because it captures the long-time behaviour of the system.A generic initial wavefunction evolving for a sufficiently long time under exp(−iH t) will converge to the steady-state; we shall study the phase diagram of this steady-state.The Hamiltonian also commutes with total z magnetization m = N j=1 S z j , so that it can be block diagonalized in each sector with a definite total magnetization.Furthermore, there is another symmetry operator P X with X = N j=1 σ x j which sends m to −m, and therefore it suffices to study non-positive magnetization (m ≤ 0) states.For the odd length chain, P X symmetry leads to the two-fold degeneracy of the steady-states.Thus, we only take even site number N throughout the paper.
Ferromagnetic (∆ > 0) and anti-ferromagnetic (∆ < 0) models can be related by the transformation As such, an eigenstate of ferromagnetic Hamiltonian with H(∆, g)|ψ〉 = E|ψ〉 can be transformed to an eigenstate of the anti-ferromagnetic one: . If E has the largest imaginary part among the spectrum, so does −E * .Thus, the steady-state properties remain the same for ±∆, and it suffices to study the ferromagnetic case.
Eigenstates of the model can be solved by coordinate Bethe Ansatz [1].Taking all spin down state as the reference state with energy E 0 = − 1 4 (N − 1)∆, we can excite M magnons by flipping M spins up (M ≤ N /2).We can construct the Ansatz state |ψ〉 , whose wavefunction in the onsite magnon number basis is given by:  1) in the zero magnetization sector.The critical curve ∆ 2 + g 2 = 1 separates P T exact and broken phases.When |∆| < 1 and g > g c = 1 − ∆ 2 , the steady-state is the many-body scale-free state; when |∆| > 1, the steady-state is the boundary string state.
the momentum of magnon by β j = e ik j [49,64], we have the equivalent form: The coefficient A({β j }; P, {η j }) is a function of magnon momenta {β j }, permutation P, and chirality {η j }.Imposing the condition that |ψ〉 is an eigenstate of H with energy and the boundary condition, these coefficients can be found as [1]: where Here, the M momenta have to satisfy the so-called Bethe equations: Intuitively, the left hand side is the phase accumulated by a magnon when it travels freely from one end of the chain to the other and then gets reflected back; each term of the right hand side represents the scattering phase between magnon j and l.Note that for the periodic XXZ model there is only the β j , β l term in the scattering phase, while for the open boundary condition β j also scatters with β −1 l .The solutions are inversion-symmetric in the sense that β j and β −1 j always appear in pairs in the solutions.After solving a set of Bethe roots {β j }, the energy of M magnon state is Eq. ( 4).Notably, if e., all β j 's belong to the unit circle, one must have Im(E M ({β j })) = 0. Thus, PT symmetry breaking requires |β j | ̸ = 1 for certain β j 's, which implies the presence of NHSE (which turns out to be scale-free here).It has been proved that the Bethe Ansatz technique can produce a complete set of solutions for finite-length models with such boundary terms [65].Here, we specifically obtain the closedform expression for the complex momentum distribution on the steady-state and the ground state in the thermodynamic limit.
Our results on the steady-states for different parameters are summarized in Fig. 1.We focus on the zero magnetization sector (m = 0) since for the Hermitian X X Z model (with ∆ < 1) the ground state lies in this sector.Our numerical results support that the steady-state belongs to the zero magnetization sector.The phase boundary between P T -exact and broken phase is ∆ 2 + g 2 = 1.Given the transformation between H(∆, g) and H(−∆, g), the steadystate phase diagram is symmetric about the g axis.Apart from the many-body scale-free modes which will be investigated in detail in Section 3, we identify the "boundary string" state as the steady-state when |∆| > 1.A boundary string corresponds to a multi-magnon bound state localized at the boundary, which will be explained in Section 3.4.

Bethe Ansatz solutions for scale-free skin modes 3.1 Single-magnon state
In the single-magnon sector (M = 1), the Bethe equation ( 5) is simplified to When |∆ ± | < 1, or equivalently g < g c , solutions of the equation have been found to be on the unit circle [66].We briefly review the proof.We define complex variable function Eq. ( 6) is then transformed to h(β) = h(β −1 ).For g < g c , the image of the disk |β| < 1 under h is still inside the disk, that is, |h(β)| < 1, and vice versa.The statement can be verified by writing β = ρ exp(iφ), ∆ + = ρ 0 exp(iφ 0 ) with ρ 0 < 1, then A similar argument works for |β| > 1.Thus, possible solutions of h(β) = h(β −1 ) must be on the unit circle, corresponding to purely real momentum.When g > g c , no theorem prohibits the existence of non-unitary solutions, and one can notice that a pair of isolated boundary modes with β ± ≈ ∆ ± is possible.For such a solution, |β ± | > 1 leads to the divergence of the term β 2(N −1) in Eq. ( 6) in the thermodynamic limit, but this can be compensated by the factor (β − ∆ + )(β − ∆ − ), which is close to zero.Most of the energy levels remain real, and the scale-free localization behavior is absent in this sector.We define the boundary imaginary energy contributed by one boundary mode as

Two-magnon state
In the M = 2 sector, scale-free modes appear and contribute to the 1/N scaling behaviour of energy.Fig. 2(a1) shows a typical finite-size two-magnon spectrum.We color the eigenvalues by many-body participation entropy [67][68][69], a measure of localization generalized from single-particle inverse participation ratio (IPR): where the index i sums over all the basis functions in the relevant Hilbert space, |ψ〉 is a manybody eigenstate.We choose |i〉 to be local magnon number basis for the following calculations.The participation entropy gets smaller when the eigenstate is more localized in the Hilbert space, as we observed in Fig. 2(a1): Two dark points correspond to isolated states bounded to the boundaries, and both two magnons are localized to the same side; around Im(E) = ±E b there is a continuum of states, which are the combination of a boundary mode with β 1 ≈ ∆ ± and a scale-free mode with β 2 ≈ e ik ; the continuum on the real axis is brighter, though there is one localized mode, corresponding to the state with two magnons localized on different ends.We apply Bethe equation to the state with one magnon bounded to the boundary while the other has momentum β 2 : The right hand side, which will be denoted by exp(2κ), is of order O(1).Taking logarithm of both sides, we have 2N ln|β 2 | = 2κ + O(1/N ).Thus, lnβ 2 acquires a first order correction to the real part: Traveling along the chain, the corresponding magnon accumulates an amplitude change A = |β 2 | N ∼ exp (κ).This magnon localization is distinguished from the disorder-induced Anderson localization and the original NHSE, both of which have exponential eigenstate decay ψ(x) ∼ exp(κ ′ x) with size-independent localization lengthes, so that the wavefunction amplitude change (∼ exp(κ ′ N )) diverges as the system size N grows to infinity.In contrast, the amplitude change A in our case saturates as the size N grows.In terms of the generalized Brillouin zone (GBZ) of the non-Bloch band theory [49,64], ln|β 2 | ≈ κ/N means that the GBZ (more precisely, the finite-size GBZ [70]) deviates from the unit circle by amount O(1/N ).Notably, the scale-free localization in our model has an intrinsic many-body origin because in the non-interacting limit ∆ = 0, the right hand side of Eq. ( 8) equals 1, and the imaginary momentum vanishes.Therefore, the phenomenon in our model is dramatically different from that of free-particle models [46][47][48].
The energy of such two-magnon state is The statement is verified by finite-size scaling of the average of those complex energy.As shown in Fig. 2(b1), the imaginary part scales linearly with 1/N .In the thermodynamic limit, it will converge to ±E b .Moreover, we add a next-nearest-neighbor zz interaction to break integrability, and the numerical results are shown in Fig. 2(b1)(b2).After adding the integrability breaking terms, the many-body scattering process can never be factorized into consecutive two-body scattering.However, on the two-body level the Bethe equations (for β 1 and β 2 ) are still valid with minor modifications on the exact form.Specifically, the continuum of states out of the real axis in Fig. 2(b1) are all the combinations of a boundary mode β 1 ≈ ∆ + ∆ ′ ± i g and a scale-free mode β 2 ≈ e ik .Therefore, Eq. ( 9) and Eq. ( 10) can still capture the qualitative behaviour of the scale-free modes after adding H nnn .

Imaginary Fredholm equation at 0 < ∆ < 1
After the warm-up on two-body scale-free modes, we now generalize it to the many-body cases.We will study the parameter space 0 < ∆ < 1, where it is the Luttinger liquid phase in the Hermitian limit and the ground state lies in M = N /2 sector (with zero magnetization).
We assume that the steady-state is composed of a boundary mode and a set of continuous scale-free Bethe roots, and then derive the Bethe equations in the thermodynamic limit.We adopt a conventional parametrization of magnon momentum [71,72]: where γ = arccos(−∆) such that 0 < γ < π.The kinetic energy of the magnon is Taking the logarithm of Bethe equations ( 5), we have where the function θ n is defined as The second term φ b (x j ) on the left of Eq. ( 14) is the scattering phase between the magnon j and the boundary, which has the form: where θ m ± is obtained by taking n in Eq. ( 15) as complex numbers This involved boundary term will not have significance in the rest part of solution.For the right hand side, I j is an integer and the set of {I j } determines the set of Bethe roots.We take the occupation of one boundary mode and I j = j on the steady-state, and the corresponding boundary Bethe root β 0 = ∆ − has the parametrization The steady-state Bethe equations become We rewrite x j = λ j + iσ j /N with purely real λ j , σ j .We note that |σ j /N | ≪ λ j , and therefore any real function of x j can be expanded as The real and imaginary part of Bethe equations are: Eq. ( 19) and Eq. ( 20) is accurate only to the O(N ) and O(1) order, respectively, and their leading terms are: Eq. (22) indicates that each scattering between the ( j, l) magnon pair generates a O(1/N ) contribution to σ j , and therefore the sum over l is of order O(1).Since a nonzero σ j ∼ O(1) characterizes the scale-free NHSE (with an imaginary part of momentum ∼ σ j /N ), Eq. ( 22) clearly demonstrates the many-body origin of the scale-free NHSE in the present model.Eq. ( 21) is the same as the ground state Bethe equations of the Hermitian open XXZ model [1].It is standard to calculate the difference between ( j +1)−th and the j −th equation, taking f (λ j+1 ) − f (λ j ) = f ′ (λ j )(λ j+1 − λ j ): The thermodynamic limit is taken by sending lim then the integral equation of ρ(λ) is where K n (λ) is the derivative of θ n (λ) with respect to λ: Notably, the energy function is proportional to K 1 up to a constant: Note that only when the steady-state belongs to the zero magnetization sector, the integral interval of λ can be taken as (−∞, ∞).This is the case here, and it corresponds to filling the Fermi sea k ∈ (−π + γ, π − γ).The integral equation ( 25) is commonly solved by Fourier transformation: Applying the convolution formula on Eq. ( 25), we have a linear equation then the distribution function is solved: ρ(ω) = 1/(2π cosh ω), ρ(λ) = 1/(2 cosh(πλ/2)).Eq. ( 22) counts two kinds of mechanisms of scale-free localization.On the right hand side, the first term sums interactions between magnons in the bulk, while the second term is the scattering with the boundary mode.The last one is much smaller than the first in a few-body state, but becomes comparable when the magnon number is the same order of system size N , e.g. in the zero magnetization sector.Define σ(λ) as a function of λ, the continuous version of this equation is Denoting the above expression by Θ(ω), we can solve σρ (ω): The summation of energy of all those scale-free modes becomes an integral over λ in the thermodynamic limit.The imaginary part of energy formula of a single magnon is Here, Eq. ( 26) has been used.While each one contributes O(1/N ) to the total imaginary part, the sum of contributions from all scale-free magnons is comparable to the boundary mode contribution E b : The imaginary part of the steady-state eigen-energy is then given by adding E b : In Fig. 3, we compare our formula with the exact diagonalization (ED) results in the M = N /2 sector (zero magnetization sector), which agrees excellently.The boundary field g controls the imaginary part of the energy totally via λ 0 = 1 γ ln g−g c g+g c .As g crosses g c , the steady-state energy becomes complex.

Boundary bound state at ∆ > 1 and phase transition
Anisotropic interaction ∆ > 1 prefers bounding all magnons together, and therefore the magnons in the steady-state tend to localize to the boundary.Specifically, the first magnon is bounded to the boundary, and the next one is bounded to the previous one recursively.In the context of integrable spin models, the bound state is named a "string"; we shall follow this terminology and call our bound state near the boundary a "boundary string".We note that similar states have been identified in the spin model subjected to a non-Hermitian magnetic field at only one end [24].In the thermodynamic limit, Bethe roots {β j } satisfy a recursive relation: The imaginary part of energy is given by For large N , β N /2 approaches the fixed point of recursive relations Eq. ( 35): It is clear that the structure of Bethe roots of the steady-state is different for 0 < ∆ < 1 and ∆ > 1.We may take the isotropic limit ∆ = 1 from the two sides to understand the phase transition point.In the scale-free phase, we have to deal with the γ → π limit of Eq. ( 33) carefully.Note that the Fermi sea ranging from −π + γ to π − γ shrinks to a Fermi point when γ → π.We take the limit by substituting ω by ωγ/ sin γ with sin γ ≪ 1, and the integral can be simplified to In the boundary-string phase, the imaginary part is a constant γ 0 = g/2.Bethe roots can be determined by the recursive equation Eq. ( 35) at ∆ = 1, which results in an explicit solution: For small n, the magnon localizes exponentially at the boundary.However, for n sufficiently large (n/N ∼ O(1)), it behaves in a scale-free fashion.This result is also confirmed by comparing the imaginary part with γ 0 for different system sizes, and the differences scale linearly with 1/N (see Fig. 4).We emphasize that the solution Eq. ( 37), though qualitatively valid, is not exact because the scattering between the large n scale-free modes have been neglected.This approximation has been implied in Eq. ( 35).

Ground state phase diagram
For the ferromagnetic case ∆ > 0, the steady-state and the ground state coincide, and the phase boundary ∆ = 1 separates the scale-free phase and the boundary string.However, the transformation Z T relating the ferromagnetic and the anti-ferromagnetic steady-state is not applicable to the ground state, because it changes to the highest-energy (real part) state when reversing the sign of ∆.Therefore, one cannot borrow the ground-state phase diagram from that of the steady-state.Moreover, the comparison between ground states at zero and finite non-Hermiticity provides another perspective on the effect of boundary dissipation.In Fig. 5, we summarize the results of the ground state.The Ansatz of many-body scale-free state in the region ∆ < 0 is studied in the following subsections.1).In the P T exact phase (g < g c ≡ 1 − ∆ 2 ), the ground state is Luttinger liquid (LL); for g > g c , gapless excitations become scale-free.When ∆ < −1, the gap between the ground state and excited states opens, yet our Ansatz of scale-free skin modes remains valid.

Scale-free solutions for −1 < ∆ < 0
For −1 < ∆ < 0, we find that Eq. ( 33) can still be applied to the ground state, though it does not coincide with the steady-state anymore.The formula is compared with exact diagonalization results in Fig. 6.It seems that the data does not agree as well as in the ferromagnetic case (see Fig. 3).This is due to the larger finite-size error.In fact, we observe that as size N increases, the ED results become closer to the analytical ones.We also show the numerical results obtained by solving the Bethe equations (5) directly for N = 40 (see Sec.A for the detailed process of calculations).The results match the analytical solutions better than ED, thus supporting our conjecture of the finite size effect.
The single-magnon kinetic energy is  We also adopt here a new definition of the function θ n : The imaginary Fredholm equation is then given by Since functions of λ are periodic functions with periodicity 2π/φ, we can expand them as The imaginary part of energy is: As illustrated in Fig. 7, there are finite-size errors between the numerical results and our Bethe Ansatz formula, which is similar to the gapless antiferromagnetic case.

Experimental scheme
The onsite non-Hermiticity can be realized in cold atom systems by coupling the spin down (up) degrees of freedom of the first (last) site to an auxiliary state by optical pumping [73,74].Each atom in the bulk has two effective energy levels to mimic a 1/2 spin.XXZ interaction can be implemented by the anisotropic spin-exchange interactions adjusted by Feshbach resonances [75,76], or by the dipole-dipole interactions between different parity states in Rydberg atoms [77][78][79][80].We may introduce the third energy levels on the two ends so that the effective spin down (up) state on the left (right) end can decay to it spontaneously.The effective loss is described by a non-Hermitian term To evolve the open system under the non-Hermitian Hamiltonian without quantum jump, i.e., by post-selection, the population of auxiliary energy levels should be monitored by exciting the states with laser.The absence of fluorescence signals the absence of the quantum jump from the magnetization-conserved quantum trajectory.Starting from an initial state in the zero magnetization sector, the system will relax to the steady-state after sufficiently long time.The imaginary part of the corresponding eigenvalue can be obtained by measuring the expectation of boundary spin polarization: where |ψ s 〉 is the steady-state.The measured boundary spin polarization can be compared to our analytical result of Im(E s ).
Numerical simulations for post-selection evolution under H eff are conducted on N = 14 chain to back up the above proposal.We consider two kinds of initial states.The first one is a "local quench", in which the spin chain is prepared in the ground state of H XXZ , and boundary coupling to the auxiliary energy levels is turned on at certain moment.The other initial state is a domain-wall configuration, in which the spins of the left half chain point down while those of the right half point up.We discretize the continuous time evolution by fourth order Runge-Kutta method, and obtain the spin polarizations in Fig. 8(a-d).It is clear that for both local quench and domain wall initial states, the edge spin polarization converges to the predictions of Bethe Ansatz solution.We also show the steady-state expectation value of local spin polarizations through exact diagonalization in Fig. 8(e)(f).
Here, we compare the relaxation dynamics and steady-state behaviour between different anisotropy ∆.It is shown clearly that the steady-state expectation of S z N is well below 1/2 for |∆| < 1 , while it saturates to 1/2 for |∆| > 1.Furthermore, the polarization increases smoothly from the left end to the right in Fig. 8(e), since the steady-state Bethe roots include only one exponentially localized magnon together with scale-free modes, which contribute to a mildly unbalanced spatial profile of the polarization.Conversely, for |∆| > 1 [Fig.8(f)], N /2 exponentially localized modes on the steady-state lead to the approximate domain-wall configuration, with sharp transition of the polarization in the middle of the chain.All of the above confirm our theoretical predictions from the Bethe Ansatz.The relaxation time also depends on ∆ significantly: while in the scale-free phase [Fig.8(a)(c)] it typically takes 60 ∼ 80 units of time (1/g) for the boundary polarization to approach its steady-state value due to the gapless spectrum near the steady-state, the relaxation time can be comparable to 1/g for |∆| > 1, as shown in Fig. 8(b).For the later case, although the boundary-string steady-state degenerates exponentially with other shorter string states, those states share the same spatial profile of the polarization, so that the relaxation time is proportional inversely to the imaginary gap between those boundary-string states and other eigenstates without the boundary mode, i.e. g/2.

Conclusion
In this work, we applied coordinate Bethe Ansatz to solve the steady-state and the ground state of a P T -symmetric one-dimensional boundary-dissipated spin chain, focusing on the P T broken phase.We found the many-body scale-free state, which is composed of one boundary mode and a continuum of scale-free modes in our particular model.We then derived the Bethe equations of the scale-free Bethe roots, and obtained a compact formula for the eigenenergy in the thermodynamic limit.We then proposed an experimental scheme to measure the dissipative part of the energy, and discussed how to compare it with our analytical results.
Our findings shed a light on exceptional points and P T transition in many-body physics.Particularly, our solution is a generalization of the concept of scale-free NHSE from free-particle to many-body systems.Although we focused on the scale-free behaviour in the XXZ spin chain, it is expected that this feature is universal in a family of non-Hermitian models with interactions in the bulk and dissipation-induced defect mode at the boundary.For example, non-integrable models also exhibit scale-free properties, as demonstrated in Fig. 2.Moreover, in integrable models solvable by nested Bethe Ansatz (Fermi-Hubbard model, higher spin XXX chain, etc.), boundary-operator-induced P T symmetry transition and the corresponding steady-states may have richer structures to uncover.
We have demonstrated the scale-free skin effect by the difference between the imaginary part of the steady-state energy and that of a single boundary mode.Intuitively, the scale-free skin effect may also manifest itself in the excitation spectra near the steady-state.A thorough analysis of those eigenstates requires preserving the O(1/N 2 )-order terms in the imaginary Fredholm equations, which is left for future studies.Algebraic Bethe Ansatz is another possible approach to the solutions of excitations, though it remains a question how scale-free Bethe roots emerge in the monodromy and transfer matrices.

Figure 1 :
Figure 1: Steady-state phase diagram of model Eq.(1) in the zero magnetization sector.The critical curve ∆ 2 + g 2 = 1 separates P T exact and broken phases.When |∆| < 1 and g > g c = 1 − ∆ 2 , the steady-state is the many-body scale-free state; when |∆| > 1, the steady-state is the boundary string state.

Figure 3 :
Figure3: Imaginary part of the steady-state energy E s as a function of the non-Hermitian parameter g.We take ∆ = 0.8, so that g c = 0.6.Red curve is from the Bethe Ansatz (BA).As a comparison, blue line is imaginary energy of a boundary mode (E b ).The black dots are computed from exact diagonalization (ED) in the zero-magnetization subspace, with system size N = 16.

6 Figure 4 :
Figure 4: (a) The imaginary part of the steady-state energy in the zero magnetization sector, which are obtained from exact diagonalization with system size N = 16.(b) Finite-size scaling of the steady-state energy at the isotropic point ∆ = 1.γ 0 = g/2 [see the paragraph above Eq.(37)].

Figure 5 :
Figure 5: Ground state phase diagram in the zero magnetization sector.The right half of the diagram, ∆ > 0, coincides with the counterpart of steady-state (see Fig.1).In the P T exact phase (g < g c ≡ 1 − ∆ 2 ), the ground state is Luttinger liquid (LL); for g > g c , gapless excitations become scale-free.When ∆ < −1, the gap between the ground state and excited states opens, yet our Ansatz of scale-free skin modes remains valid.

Figure 6 :
Figure6: Imaginary part of the ground state energy E g as a function of the non-Hermitian parameter g.We take ∆ = −0.8,so that g c = 0.6.Red curve is from the Bethe Ansatz (BA).As a comparison, blue curve is imaginary part of the energy of a boundary mode (E b ).The diamonds are obtained by solving the Bethe equations numerically (Num BA).The hollow dots are obtained from exact diagonalization (ED) with N = 10, and solid dots with N = 16.

Figure 7 :
Figure 7: Imaginary part of the ground state energy as a function of the non-Hermitian parameter g.Left panel: ∆ = −1.2;Right panel: ∆ = −2.0.Red curve is obtained from the Bethe Ansatz (BA).Green triangles and blue squares are numerical results from exact diagonalization of N = 10 and N = 16, respectively.

Figure 8 :
Figure 8: (a)-(d) Time evolution of spin polarization under post-selection dynamics.The dashed lines indicate the values of Im(E s ) g obtained from the Bethe Ansatz.For (a)(b), the time evolution starts from a local quench; in (c)(d), it starts from the "domain-wall" initial state (see text).The time is measured in unit of 1/g.(e)(f) Steady-state spin polarization profile obtained from exact diagonalization.Parameter values: N = 14 and g = 0.8 are fixed.For (a)(c)(e), ∆ = 0.8; for (b)(d)(f), ∆ = 1.2.