Open XXZ chain and boundary modes at zero temperature

We study the open XXZ spin chain in the regime Delta>1 and for generic longitudinal magnetic fields at the two boundaries. We discuss the ground state via the Bethe ansatz and we show that, in a regime where both boundary magnetic fields are equal and bounded by a critical field, the spectrum is gapped and the ground state is doubly degenerate up to exponentially small corrections in the length L of the chain. We connect this degeneracy to the presence of a boundary root, namely an excitation localized at one of the two boundaries. We compute the local magnetization at the left edge of the chain and we show that, due to the existence of a boundary root, this depends also on the value of the field at the opposite edge, even in the half-infinite chain limit. Moreover, we give an exact expression for the large time limit of the spin autocorrelation at the boundary, which we explicitly compute in terms of the form factor between the two quasi-degenerate ground states. This, as we show, turns out to be equal to the contribution of the boundary root to the local magnetization.


Introduction and main results
The study of condensed matter theory involves understanding many-body systems starting from their elementary constituents. This protocol, which is in general notoriously hard, can be sometimes carried out in systems of one-dimensional spin chains. These constitute one the main theoretical playgrounds for the emergent physics of strongly correlated quantum systems, see for example the seminal work of Haldane [1]. In particular, in the past years, a class of interacting spin chains which can be exactly solved by the so-called Bethe Ansatz [2,3] have been successfully applied to understand the dynamical response of real compounds [4,5] or to develop better numerical techniques [6].
While the bulk physics of spin chains can be usually studied by considering the large-size limit of systems with periodic boundary conditions, a richer phenomenology can be observed in the presence of open boundaries. By tuning different parameters at the boundaries one can explore different phase transitions (also experimentally [7]), as well as probing the existence of boundary modes. Notoriously, in topological superconducting systems, the Majorana zero modes [8] are boundary modes and they consist into two decoupled Majoranas fermions localized at the two edges of the system that can be combined to form a zero-energy regular fermion. As a consequence of their existence, all many-particle states are degenerate. While Majorana zero modes are present in the so-called Kitaev chain, which becomes the XY chain with a transverse field after a Jordan-Wigner transformation, it was recently shown by Fendley [9] that also the gapped (massive) XYZ chain contains strong zero modes, namely operators defined at the two edges of the chain that commute with the Hamiltonian up to exponentially small corrections with the size of the chain. These operators, instead of being exactly localized at the two edges, are characterized by exponential tails that decay away from the edges and are related to the Z 2 symmetry of the model. Their existence also implies an extensive number of degeneracies between the different many-body states in the spectrum.
From the physical point of view it is interesting to study the spin autocorrelation at the edge of the chain. Due to the presence of the aforementioned boundary modes, the latter should not decay to zero even at finite temperature T , in the thermodynamic limit L → ∞. Namely, given the Pauli spin operator σ z 1 at the left edge of the chain and its time evolution σ z 1 (t) = e iHt σ z 1 e −iHt with the Hamiltonian H containing a strong zero mode, one should find that, for any temperature T , T denotes the connected correlator and O T = Tr e −βH O /Tr e −βH the thermal expectation value. This prediction constituted a starting point of an active research field focused on the study of coherence time of edge spins in the open XXZ chain. When the Hamiltonian is perturbed by additional terms that do not preserve the symmetry of the zero mode, namely when the system is perturbed away from the integrable limit, it was shown [10][11][12][13] that the dephasing time can still get very large and that the spin autocorrelation remains on a long-living plateau at large intermediate times.
We here consider the open XXZ Hamiltonian with anisotropy parameter ∆ and boundary longitudinal magnetic fields h − and h + , in the massive anti-ferromagnetic regime with ∆ = cosh ζ > 1 (ζ > 0), which is indeed the regime of existence of the strong zero modes [9]. We particularly focus on the case of a chain with an even number of sites L. There are two critical values of the magnetic field at the boundary, h cr = ∆ − 1 and h (1) cr = ∆ + 1, where different crossings between eigenstates occur, see Fig. 3. In the regime where |h ± | < h (1) cr , as we shall see, the spectrum is gapped and the ground state is doubly degenerate in the large L limit whenever h + = h − = h, so that the zero-temperature spin-spin boundary autocorrelation function is expected to converge for large time to the form factor of the spin operator between these two degenerate ground states. Namely, by denoting with GS i , h |, i = 1, 2, the two normalized ground states of the open chain with boundary magnetic fields h − = h + = h, we expect that In this paper, we explicitly compute the thermodynamic and large-time limit (1.3) of the boundary auto-correlation function at zero-temperature from the study of the open chain (1.2) in the algebraic Bethe ansatz (ABA) framework [14]. By considering the large L limit of the solutions of the Bethe equations and controlling the finite-size corrections up to exponentially small order in L, we show that the difference of energy between the ground state and the first excited state becomes exponentially small in L when h + = h − ( |h ± | < h (1) cr ). Each of these two states is characterized by a Fermi sea of L 2 − 1 real Bethe roots and an isolated complex Bethe root which corresponds to a boundary mode and that we call boundary root. The latter represents a collective magnonic excitation pinned at one of the two edges of the chain, whose wave function has exponentially decreasing tails away from the boundary [15]. We show that this boundary mode is responsible for the ground state degeneracy, which in particular has two main physical consequences: 1. The boundary magnetization in the ground state for even size L depends on the value of both boundary fields, even in the infinite chain limit L → ∞ (thermodynamic limit). This is due to the fact that the presence of the boundary root in the Bethe solution for the ground state and its localization at one or the other edge of the chain depends on the values of both boundary fields. Moreover, when one of the fields is inside the interval (−h (1) cr , h (1) cr ), the boundary magnetization becomes a discontinuous function of the other field at h − = h + , point at which the boundary root jumps from one edge of the chain to the other. We here provide an analytical derivation for the boundary magnetization at the left edge of the chain, and notably for the value of its discontinuity at h − = h + = h (|h| < h (1) cr ). The latter is given by σ z 1 BR , the (thermodynamic limit of the) contribution to the boundary magnetization carried by the boundary root at the left edge, which is non-zero only when h − ≥ h + : see eq. (5.16) and (5.17) for an exact expression in terms of the parameters of the model. At exactly h − = h + , the boundary root becomes delocalized between the two edges of the chain and contributes equally to the left or the right boundary magnetization, hence the factor 2 in (1.4). In the particular case h + = 0, we recover that [16,17] lim is the bulk magnetization [18].
2. The degeneracy of the ground state implies that the matrix element of the spin operator σ z 1 in the first site of the chain between the two degenerate ground states GS 1 , h | and | GS 2 , h with fields h − = h + = h (|h| < h (1) cr ), is non-zero in the thermodynamic limit L → ∞ and its absolute value squared provides the infinite time limit of the boundary spin-spin autocorrelation function, see equation (1.3). We here exactly compute this matrix element in the ABA framework and explain how to derive its thermodynamic limit. We show that, in this limit, it is directly related to the contribution to the boundary magnetization carried by the boundary root at the left edge as for any |h| < h (1) cr , so that it is given by half of the boundary magnetization discontinuity (1.4) at h + = h − . For h + = h − = h and |h| < h (1) cr the quantity σ z 1 BR , and so the matrix element (1.6), is non-zero, see (5.44). When both fields are zero (h = 0) this reduces to the value (1.5): In Fig. 1 and Fig. 2 the boundary spin autocorrelation function σ z 1 (t) σ z 1 c T =0 is computed numerically by tDMRG as a function of time for a chain with finite size L: at large times the correlation attains the value given by (1.3)-(1.6). T =0 vs. time t obtained by tDMRG. In red: h ± = 0, in blue, h ± = 1. System size L = 100, anisotropy ∆ = 3. Corresponding exact thermodynamic limit values from eq. (1.6) and eq. (1.7) are shown dashed.
Note that a similar relation to (1.7) exists also in the bulk of the chain. There, as shown first by Baxter [18], the local magnetization of the ground state is staggered, namely we have where h is the global magnetic field. The value of s 0 is also equal to the form factor between the two degenerate ground states of the chain with periodic boundary conditions, see [19].
This article is organized as follows. In section 2, we recall the diagonalization of the Hamiltonian (1.2) in the framework of the boundary algebraic Bethe ansatz introduced by Sklyanin in [14]. In section 3, we explain how to derive, in this framework, compact determinant representations for the finite-size matrix elements (form factors) of the σ z 1 operator between two Bethe eigenstates. In section 4 we study the solutions of the Bethe equations in the thermodynamic limit L → +∞ and explain how to control their finite-size corrections up to exponentially small order in L. We identify the solution corresponding to the ground state for the different values of the boundary magnetic fields h + and h − . In the regime where both fields are between −h (1) cr and h (1) cr , with h (1) cr = ∆ − 1, we show that the two states of lowest energy are given by a particular solution of the Bethe equations with L 2 − 1 real Bethe roots and one complex Bethe root which has to be chosen between the two possible boundary roots given in terms of the boundary parameter at the left or the right end of the chain. We moreover show that, when h + = h − , the deviation between the two boundary roots becomes exponentially small in L, and so does the difference of energy between the two corresponding states. In section 5, we compute the thermodynamic limit of the determinant representation that we obtained in section 3 in two particular cases: the mean value of σ z 1 in the ground state, which gives the boundary magnetization, and the σ z 1 form factor between the two degenerate ground states identified in section 4, which gives the infinite time limit of the boundary autocorrelation function.

The integrable open XXZ spin chain
The Hamiltonian (1.2) is integrable and can be diagonalized in the framework of the representation theory of the reflection algebra [20], by means of the boundary version of algebraic Bethe ansatz introduced by Sklyanin in [14].
The key object in this approach is the boundary monodromy matrix U(λ) ∈ End(C 2 ⊗ H) where H is the space of states of the system. It is such that V(λ) ≡ U t (−λ) satisfies the reflection equation 1 , The relation (2.1) has to be understood C 2 ⊗ C 2 ⊗ H, and the subscripts parameterize the subspaces of C 2 ⊗ C 2 on which the corresponding operators act non-trivially. The parameter ζ is related to the anisotropy parameter ∆ of (1.2) as ∆ = cosh ζ.
In the case of the spin chain (1.2) with longitudinal boundary fields, the boundary monodromy matrix solution of (2.1) can be constructed from the bulk monodromy matrix T (λ) and a diagonal scalar solution of the reflection equation (2.1), More precisely, we introduce two such boundary scalar matrices, where ξ ± are some complex parameters which parameterize the left and right boundary fields h ± as h ± = − sinh ζ coth ξ ± . The boundary monodromy matrix U(λ) is then constructed as where the bulk monodromy matrix T (λ) is itself constructed as a product of Rmatrices (2.2) as Here the index a denotes the so-called auxilliary space V a C 2 , and ξ 1 , . . . , ξ L are a set of inhomogeneity parameters which may be introduced for technical convenience.
One then define a one-parameter family of commuting transfer matrices as In the homogeneous limit in which ξ = −iζ/2, = 1, . . . , L, the Hamiltonian (1.2) of the spin-1/2 open chain can be obtained as In the algebraic Bethe ansatz framework, the common eigenstates of the transfer matrices can be constructed in the form where | 0 (respectively 0 |) is the reference state (respectively the dual reference state) with all spins up. By using the commutation relations issued from (2.1), it can be shown that states of the form (2.10) are eigenstates of the transfer matrix (2.8) provided the set of spectral parameters {λ} ≡ {λ 1 , . . . , λ N } satisfies the system of Bethe equations (2.14) Here and in the following, we use the shortcut notations: The corresponding transfer matrix eigenvalue is From (2.9), in the homogeneous limit in which ξ = −iζ/2, = 1, . . . , L the transfer matrix eigenstates (2.10) become eigenstates of the Hamiltonian (1.2) with energy where the bare energy ε 0 (λ) is defined as .
Eigenstates of the form (2.10) are called on-shell Bethe states. States of the form (2.10) for which the parameters {λ} do not satisfy the Bethe equations are instead called off-shell Bethe states. The study of the solutions of Bethe equations, and in particular of the ground state of the Hamiltonian (1.2) in the thermodynamic limit, has been performed in [15,21,22].
Building on this ABA description of the spectrum and eigenstates, it is possible to compute the zero-temperature correlation functions of the open spin chain [17,23]. However, this program has not yet reached the level of achievement as what has been done for the bulk correlation functions [24][25][26][27][28][29][30][31][32][33][34][35]. In the latter case, it was indeed possible to derive the large distance and long time asymptotic behavior of the two-point (or even multi-point) correlation functions in the thermodynamic limit from their exact representations on the lattice. At the root of this approach was the fact that there exist some compact and simple determinant formulas for the form factors of local operators in the finite periodic chain [24]. Such determinant representations were also of uttermost importance for the numerical studies of the correlation functions [36][37][38][39]. They were obtained thanks to two main ingredients: a determinant representation of the scalar product of an off-shell and an on-shell Bethe states [40], and the fact that the local spin operators could be expressed as a simple element of the monodromy matrix dressed by a product of transfer matrices (solution of the quantum inverse problem) [24,41,42].
In the open case, however, such nice determinant representations for the form factors do not exist in general. It is still possible to expressed the scalar product of an off-shell and an on-shell Bethe states of the form (2.10) as a generalized version of the Slavnov determinant [17,43], but a convenient expression of the local spin operators in terms of the boundary monodromy matrix elements dressed by a product of boundary transfer matrices is presently not known, except at the first (or last) site of the chain [44]. In fact, the formulas obtained in [17,23] relied on a cumbersome use of the bulk inverse problem, which resulted into multiple integral formulas for the zerotemperature correlation functions in the thermodynamic limit (half-infinite chain) similar to the one that were previously obtained in [16] from a different approach.
At the first (or last) site of the chain, however, the situation is different. Indeed, the solution of the quantum inverse problem proposed in [44] is in that case sufficient, together with the determinant representation for the scalar products, to obtain determinant representations for the form factors of local operators at site 1 which are very similar to the bulk ones. Hence, we may expect to be able to study their thermodynamic limit similarly as what has been done in [19,28,29,34]. In particular, we are in position to compute and study the thermodynamic limit of the form factors which are relevant for the long-time limit of the boundary autocorrelation (1.1). This is the purpose of the next sections. . , µ N } and arbitrary set of parameters, the latter is given by [17,43] where the elements of the N × N matrix H(λ, µ) are for λ ≡ (λ 1 , . . . , λ N ) and µ ≡ (µ 1 , . . . µ N ). The reconstruction of the σ z 1 operator in terms of the boundary monodromy matrix elements reads [44] where ξ 1 is a generic inhomogeneity parameter that should be sent to −iζ/2 at the end of the computation. We also recall the action of the boundary monodromy matrix element A(ξ 1 ) on an off-Bethe state (2.10), which follows from the commutations relations issued from (2.1): with It follows from (3.4), (3.5) and (3.1) that the matrix element of the σ z 1 operator between two eigenstates {λ} | and | {µ} is where H(λ, µ) is the matrix (3.2) and P (λ, µ) is a rank one matrix with elements . (3.9) So as to express the determinant in a more convenient form for taking the thermodynamic limit, let us introduce, as in [19], an N × N matrix X with elements . (3.10) Multiplying and dividing (3.8) by det X , computing the matrices X H and X P , and factorizing the quantity outside of the determinant, we obtain , (3.14) , (3.15) in which we have defined Using the Bethe equations for {µ} and taking the limit ξ 1 → −iζ/2, we can rewrite (3.14) and (3.15) as (3.18) in which we have set . (3.20) It remains to take into account the normalization of a Bethe state, which is given by the formula (3.21) The matrix M(λ, λ) reads explicitly in which ξ (µ|{λ}) is the following meromorphic function: with p and K given respectively by (3.20) and (3.19), and with The ground state(s) in the thermodynamic limit In this section, we explain how to characterize the configuration of Bethe roots for the ground state(s) of the open XXZ Hamiltonian (1.2) in the regime ∆ > 1. As we shall see, the total number of these Bethe roots and their pattern in the complex plane for large L depend non-trivially on the values of the magnetic fields at the boundaries, and so does the presence of an energy gap and of an exponential double degeneracy at h + = h − , see Fig. 3. Hence, we now focus on the regime ∆ > 1. We use the following parametrization for the anisotropy parameter ∆ and the boundary fields h σ (σ ∈ {+, −}) in this regime: whereξ σ ∈ R, and We also suppose that the number of sites L of the chain is even 2 . The Bethe equations (2.11) can be conveniently rewritten 3 as in terms of the function (3.16). In the homogeneous limit ξ n → −iζ/2, n = 1, . . . , L, the latter reads explicitely Due to the parity and periodicity properties of the problem, we can in fact restrict ourselves to the roots which are contained in the following subspace of the complex plane: The ground state for the open XXZ chain in the regime ∆ > 1 was studied in [15]. It is given by a solution of the Bethe equations where all Bethe roots are real, except a possible isolated complex root. In the thermodynamic limit L → +∞, the real roots α j of the Bethe equations for the ground state form a dense distribution on the interval (0, π 2 ) (which can be extended by parity on the interval (− π 2 , π 2 )), with density ρ(α) solution of the integral equation: The latter can be solved explicitly as where ϑ i (α, q), i ∈ {1, 2, 3, 4}, are the Theta functions of nome q defined as in [45], with ϑ 1 ≡ ϑ 1 (0, q), ϑ 2 ≡ ϑ 2 (0, q). It was moreover argued in [15] that the possible additional complex root was issued from the presence of the boundary factors in (4.5). 2 The degeneracies of the ground state in the case L odd are different: in that case, we do not have any more quasi-degenerate ground states for h− = h+, but an exact degeneracy for h− = −h+ due to the Z2 symmetry of the model. 3 When doing this, we have to exclude the possible roots 0 and π 2 which are always solutions of (4.4) but should actually correspond to a zero of order 2 in the numerator of (2.11). By treating them apart, it is in fact easy to see that low-energy states do not contain these roots for large L.
More precisely, according to the study of [15], the latter should correspond to a root which approaches, in the large L limit, one of the two zeroes of the boundary factors of the Bethe equations (4.5) : with exponentially small correction σ = O(L −∞ ) so as to compensate the exponentially large factor in L in the first line of (4.5). Such a complex root α σ BR , which in the following will be called boundary complex root or more simply boundary root 4 , was predicted [15] to exist only ifξ σ < ζ/2, i.e. when the corresponding boundary field h σ is not in the interval delimited by the two boundary critical fields h (1) cr and h (2) cr defined as [15,16] The presence of this kind of boundary root in the ground state was also discussed in [15], in particular in the regime h − > 0, h + < 0. It is however not completely clear, from [15], what is the accurate configuration of the Bethe roots for the ground state and the first low-energy states according to the values of the two boundary fields h ± , notably in our case of interest h + = h − (see Fig. 3) for which we may a priori expect a degeneracy. In the remaining part of the present section, we therefore perform a more detailed study of these configurations, so as to make more precise (and sometimes slightly correct) the predictions of [15]. We in particular show how to control the finite-size corrections up to exponentially small order in L, which enables us to discuss the degeneracy at h + = h − .

Properties of low-energy states for large L
Low-energy states are given in the thermodynamic limit L → ∞ by an infinite number of real roots (i.e of order L/2) and a finite number of complex roots. Using the same argument as in [46], we can show that, if a set of solutions {λ} ≡ {λ 1 , . . . , λ N } of (4.4) contains a complex root λ such that (λ ) = 0, π 2 , then it also contains the conjugate rootλ . Hence complex roots appear by pairs λ ,λ , except possible isolated imaginary roots λ such that (λ ) = 0, π 2 .

Bethe equations for real roots and counting function
Let us consider a real root λ j ∈ {λ}. It is convenient to rewrite the corresponding Bethe equation in logarithmic form, where n j is an integer and ξ(α|{λ}) is the counting function. The latter is defined, for the given set of Bethe roots {λ}, as cr and the gray states are other states with different values of total magnetization. Here for a chain with L = 12, ∆ = 4. The spectrum is gapped and there are two exponentially degenerate ground states in the region −h (1) cr < h < h (1) cr , while in all other regions the spectrum is gapless and there is no exponential degeneracy of the ground state.
in terms of the functions where we have set . (4.17) Here we have used the fact that the complex roots λ k always appear in pairs λ k ,λ k , except if (λ k ) ∈ {0, π 2 }. Note that the functions p (3.20) and g (3.24) that appeared in the expression of the form factor correspond indeed to the derivatives of p and g, and that the function K (3.19) is related to the derivative of θ as K(α) = − θ (α) 2π , so that the function ξ (3.23) is indeed the derivative of (4.12).
It is possible to determine the range of allowed quantum numbers n j for the real roots in (4.11) by continuity arguments from the Ising limit ζ → +∞. This is done in appendix A. We obtain that 1 ≤ n j ≤ M − 1, where M is given by (A.7). Hence we can rewrite the logarithmic Bethe equations (4.11) for the real roots as where M is given by (A.7) and where h 1 , . . . , h n are the positions of the holes in the adjacent set of quantum numbers for the real roots. It is also convenient to define the rapiditiesλ h k of the holes from the relation In the thermodynamic limit L → +∞, the derivative (3.23) of the counting function (4.12) tends to the density function (4.8) solution of (4.7) multiplied by π. This comes from the fact that, in the thermodynamic limit, the sums over real Bethe roots turn into integrals with measure given by the density function (4.8). As explained in appendix B, it is possible to control more precisely the finite-size corrections to this transformation sum-integral, in the spirit of what was done in [19,47] (see Proposition B.1 and Corollary B.1), and to decompose the counting function (4.12) in the large L limit according to the different contributions of the real roots, complex roots and holes up to exponentially small corrections in L: (4.20) in which the first sum runs over the set Z of indices corresponding to the complex roots (i.e. (λ k ) = 0 if k ∈ Z), whereas the second sum runs over the positions of the holes. In (4.20), the term stands for the contribution of the "Fermi sea" of real roots, taking into account the finite-size corrections

Bethe equations for complex roots, boundary roots and wide roots
We now investigate the large L behaviour of the complex solutions of the Bethe equations, and more particularly of the isolated complex roots which, according to [15] and to the study of appendix A, may appear in the set of solutions corresponding to the ground state. Hence, let us now suppose that λ j ∈ {λ} is a complex root. One can still use Corollary B.1 to rewrite the sum over real Bethe roots as integrals in the corresponding Bethe equation for large L, which gives In (4.22) we have set 23) and the functions p and θ are defined such that and such that they coincide with the definitions (4.13) and (4.14) for α real.
It is interesting to investigate the behaviour of (4.23) so as to see how the first line of (4.22) behaves with L. Using the terminology of [46], we find that if λ j is a close root, i.e. if | (λ j )| < ζ. The real part of (4.25) is moreover positive if −ζ < (λ j ) < 0 (see Fig. 4), which means that, in that case, the first factor in (4.22) is exponentially diverging in L. Hence, for (4.22) to be satisfied, λ j has to approach a zero of the expression with exponentially small corrections in L. If we suppose moreover that λ j is the only complex root of the set {λ}, i.e. that all other roots are real, this means that λ j has to approach one of the two zeros of the boundary factor in the last line of (4.22), i.e. that λ j is indeed a boundary root of the form (4.9). This can of course only be possible If instead λ j is a wide root, i.e. if | (λ j )| > ζ, then, using similar arguments as in [46], we find that (F (λ j )) = 0, (4.27) so that the first factor in (4.22) remains finite. If we suppose moreover that λ j is the only complex root of the set {λ}, i.e. that all other roots are real, this means that λ j does no longer converge exponentially fast towards one of the zeros (or poles) of the boundary factor in the last line of (4.22), and therefore is not strictly speaking a boundary root as defined in (4.9). Let us finally remark that we have here found a domain of existence of the boundary root, given by (4.26), which is more narrow than the one cr ]) found in [15]. This comes from the fact that the reasoning of the authors of [15] did not take into account the full exponential factor given by (4.23), but only the part given by p(z) 5 .

Expression of the energy
Proposition B.1 can also be applied, together with (4.20), to obtain an asymptotic expansion of the energy (2.17) of the corresponding Bethe state in the large L limit up to exponentially small order in L: where E 0 is the contribution of the real roots taking into account the finite-size corrections which are common to all low-energy states, see (B.25), whereas ε(µ) is the dressed energy of an excitation with rapidity µ, defined as Explicitly, in terms of the meromorphic elliptic function ρ(α) given by the ratio of Theta functions (4.8) if µ stands for the rapidity of a hole or of a close root (i.e. if | (µ)| < ζ), whereas in the case of a wide root (i.e. if | (µ)| > ζ), see (B.27).
In particular, the dressed energy of a hole with rapidityλ h ∈ (0, π 2 ) is as found in [15,46]. Note that (4.32) is a positive and decreasing function ofλ h on the interval [0, π 2 ], see The dressed energy of the boundary root (4.9) is (4.33) in its domain of validity (4.26), as found in [15]. We recall that δ σ is given by (4.3). Note that the expression (4.33) is an odd function ofξ σ (and therefore of h σ ). It is moreover a decreasing function of cr ], and we have so that the dressed energy (4.33) of the boundary root can be compared to the dressed energy (4.32) of a hole with rapidityλ h as The dressed energy of the boundary root α σ BR is plotted as a function of the field h σ for a specific value of ∆ in Fig. 6.
We finally recall that, according to [46], the bulk close complex roots are arranged either in 2-strings or in quartets whose dressed energy vanishes.
cr , and it tends to the one of a hole with rapidityλ h = 0 when h → −h (2) cr , h < −h (2) cr .

Configuration of Bethe roots in the ground state
We now discuss the configuration of the Bethe roots for the ground state according to the values of the two boundary fields h + and h − .
Let us first suppose that the boundary field of maximal absolute value is nonpositive, so as to ensure that the magnetization of the ground state is non-negative 6 . It follows from the study of the previous subsection that the ground state should be given by a configuration of Bethe roots that minimizes the number of holes, except if it is at the cost of containing a boundary root with higher energy, see (4.35). The allowed configurations of Bethe roots according to the values of the boundary fields in the different magnetization sectors can be deduced from the study of appendix A. Combining the results of appendix A with those of subsection 4.1, we can therefore distinguish seven different cases 7 .
The ground state is in the sector with magnetization 0 (i.e. with number of Bethe roots N = L 2 ). It corresponds to the state with L 2 −1 real roots with adjacent quantum numbers n j = 1, . . . , L 2 − 1 (no hole) and the boundary root cr , it follows from the study of appendix A that all other 6 The number of Bethe roots of a given Bethe state constructed as in (2.10) is related to its total magnetization m = L n=1 S z n as m = L 2 − N . In this framework, the states with negative magnetization would correspond to "going beyond the equator", with a number of Bethe roots exceeding the value L/2. To avoid this, one has to construct the corresponding Bethe states from the multiple action of C on the reference state | 0 with all spins down (or in the dual space from the multiple action of B on 0 |). It is also possible to reach these sectors with negative magnetization by simply using the invariance of the model under the reversal of all spins together with a change of sign of the boundary fields h±. 7 Notice that a similar classification has been recently proposed in [49] in a slightly different context. configurations with N ≤ L 2 contain either the boundary root α σ 2 BR with higher dressed energy than α σ 1 BR , or one or more hole(s), and therefore from (4.36) have higher energy. Note that the conclusion still holds even if the boundary field of maximal absolute eigenvalue is positive (but less that h (1) cr ), since by symmetry of the model under the reversal of all spins together with a change of sign of the boundary fields h ± we know that the ground state is in this case still in the sector with magnetization 0.
The consideration of the Ising limit indicates that the ground state is in the sector with magnetization 1 (i.e. with number of Bethe roots N = L 2 − 1). It therefore corresponds to the state with L 2 − 1 real roots with adjacent quantum numbers n j = 1, . . . , L 2 − 1 (one hole at position h = L 2 ). Indeed, in that case, the dressed energy of the hole is smaller than the dressed energy of the boundary root, see (4.35), so that the aforementioned state has indeed lower energy than a state with L 2 − 1 real roots and a boundary root in the sector L 2 . One can moreover notice that the latter is not even the lowest energy state in its sector, since any state with one hole (and a bulk 2-string or possibly a wide root) would also have a lower energy.
. The ground state has to be found within the states with minimal number of holes n h = 1, which may have the following configurations: (i) the state, in the sector N = L 2 − 1, with L 2 − 1 real roots and a hole at position h = L 2 ; (ii) a state, in the sector N = L 2 , with L 2 − 1 real roots with adjacent quantum numbers n j = 1, . . . , L 2 − 1, a wide root, and a hole at position h = L 2 ; (iii) a state, in the sector N = L 2 , with L 2 − 2 real roots with adjacent quantum numbers n j = 1, . . . , L 2 −2, a pair of close roots (2-string), and a hole at position h = L 2 − 1. Since the dressed energy of a wide root or of a 2-string vanishes, the difference of energy between these states is only given at leading order by the small shift between the hole rapidities. Considering the large ζ limit (see appendix A), we find that the rapidity of the hole for the configuration (iii) is given at leading order in ζ by whereas the rapidities of the hole for the configurations (i) or (ii) are given at leading order in ζ by which seems to indicate that the ground state has to be found within configurations (i) or (ii) only. To conclude further would require a more advanced study of the solutions of the Bethe equations, which is anyway unnecessary for the purpose of the present paper. Let us just mention here that the numerical results from exact diagonalization suggest that, when h σ 1 > −∆ the ground state remains in the sector m = 0 whatever the value of h σ 2 , whereas when h σ 1 approaches −h (2) cr , there exists a certain value h(h σ 1 ) < −h (2) cr at which a transition from the m = 0 (for h σ 2 > h(h σ 1 )) to the m = +1 (for h σ 2 < h(h σ 1 )) sector occurs, see Fig. 7. In red: m = +1. Notice that, as we pick h σ 1 closer to −h (2) cr , a crossing of levels appears when h σ 2 becomes small enough. Here ∆ = 3, L = 14 .
The ground state is in the sector with magnetization 0 (i.e. with number of Bethe roots N = L 2 ), even if |h σ 1 | > |h σ 2 | (this follows by symmetry from the study of previous cases). It therefore corresponds to the state with L 2 real roots with adjacent quantum numbers n j = 1, . . . , L 2 (no hole).
. This case can be obtained by symmetry from Case 2. The ground state is in the sector with magnetization −1 and hence is beyond the equator.
This case can be obtained by symmetry from Case 3. Depending on the values of the magnetic fields, it may be: (i) the state, in the sector N = L 2 of magnetization 0, with L 2 real roots and either no hole cr ]), or one hole at position cr ]); (ii) a state with magnetization −1, which is beyond the equator.

The ground state degeneracy at h + = h −
We now focus on the regime where both boundary fields h ± are such that |h ± | < h , we have seen that, in this regime, the ground state is the state in the sector N = L 2 with L 2 − 1 real roots with adjacent quantum numbers n j = 1, . . . , L 2 − 1 (no hole) and the boundary root α σ 1 BR . Moreover, it is easy to see from similar arguments that the excited state with lowest energy is the state in the sector N = L 2 with L 2 − 1 real roots with adjacent quantum numbers n j = 1, . . . , L 2 − 1 (no hole) and the boundary root α σ 2 BR . Hence, at h − = h + = h, the ground state becomes degenerated in the thermodynamic limit. Let us investigate more thoroughly this degeneracy.

e R x F O I F T O A c P r q A O d 9 C A J j A Y w j O 8 w p s j n R f n 3 f l Y t B a c f O Y Y / s j 5 / A H t M 4 2 M < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Z T c m + y z m b x 6 R N F B X v g Y 0 j b l d j 4 4 = " > A A A B 6 n i c b Z B N S 8 N A E I Y n 9 a v W r 6 p H L 4 t F 8 G J J p K D H o h e P F e 0 H t K F s t p t 2 6 W Y T d i d C C f 0 J X j w o 4 t V f 5 M 1 / 4 7 b N Q V t f W H h 4 Z 4 a d e Y N E C o O u
+ + 0 U 1 t Y 3 N r e K 2 6 W d 3 b 3 9 g / L h U c v E q W a 8 y W I Z 6 0 5 A D Z d C 8 S Y K l L y T a E 6 j Q P J 2 M L 6 d 1 d t P X B s R q 0 e c J N y P 6 F C J U D C K 1 n o Y 9 S / 6 5 Y p b d e c i q + D l U I F c j X 7 5 q z e I W R p x h U x S Y 7 q e m 6 C f U Y 2 C S T 4 t 9 V L D E 8 r G d M i 7 F h W N u P G z + a p T c m a d A Q l j b Z 9 C M n d / T 2 Q 0 M m Y S B b Y z o j g y y 7 W Z + V + t m 2 J 4 7 W d C J S l y x R Y f h a k k G J P Z 3 W Q g N G c o J x Y o 0 8 L u S t i I a s r Q p l O y I X j L J 6 9 C 6 7 L q W b 6

v V e o 3 e R x F O I F T O A c P r q A O d 9 C A J j A Y w j O 8 w p s j n R f n 3 f l Y t B a c f O Y Y / s j 5 / A H t M 4 2 M < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Z T c m + y z m b x 6 R N F B X v g Y 0 j b l d j 4 4 = " > A A A B 6 n i c b Z B N S 8 N A E I Y n 9 a v W r 6 p H L 4 t F 8 G J J p K D H o h e P F e 0 H t K F s t p t 2 6 W Y T d i d C C f 0 J X j w o 4 t V f 5 M 1 / 4 7 b N Q V t f W H h 4 Z 4 a d e Y N E C o O u
+ + 0 U 1 t Y 3 N r e K 2 6 W d 3 b 3 9 g / L h U c v E q W a 8 y W I Z 6 0 5 A D Z d C 8 S Y K l L y T a E 6 j Q P J 2 M L 6 d 1 d t P X B s R q 0 e c J N y P 6 F C J U D C K 1 n o Y 9 S / 6 5 Y p b d e c i q + D l U I F c j X 7 5 q z e I W R p x h U x S Y 7 q e m 6 C f U Y 2 C S T 4 t 9 V L D E 8 r G d M i 7 F h W N u P G z + a p T c m a d A Q l j b Z 9 C M n d / T 2 Q 0 M m Y S B b Y z o j g y y 7 W Z + V + t m 2 J 4 7 W d C J S l y x R Y f h a k k G J P Z 3 W Q g N G c o J x Y o 0 8 L u S t i I a s r Q p l O y I X j L J 6 9 C 6 7 L q W b 6      E E i u D a u + + 0 U 1 t Y 3 N r e K 2 6 W d 3 b 3 9 g / L h U U v H q W L Y Z L G I V S e g G g W X 2 D T c C O w k C m k U C G w H 4 9 t Z v f 2 E S v N Y P p h J g n 5 E h 5 K H n F F j r c Z d v 1 x x q + 5 c Z B W 8 H C q Q q 9 4 v f / U G M U s j l I Y J q n X X c x P j Z 1 Q Z z g R O S 7 1 U Y 0 L Z m A 6 x a 1 H S C L W f z R e d k j P r D E g Y K / u k I X P 3 9 0 R G I 6 0 n U W A 7 I 2 p G e r k 2 M / + r d V M T X v s Z l 0 l q U L L F R 2 E q i I n J 7 G o y 4 A q Z E R M L l C l u d y V s R B V l x m Z T s i F 4 y y e v Q u u i 6 l l u X F Z q N 3 k c R T i B U z g H D 6 6 g B v d Q h y Y w Q H i G V 3 h z H p 0 X 5 9 3 5 W L Q W n H z m G P 7 I + f w B m O + M y Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " W r b 4 3 V 9 n 9 r H g + t O I p / p D Y o W t q b s = " > A A A B 6 H i c b Z B N S 8 N A E I Y n 9 a v W r 6 p H L 4 t F 8 F Q S E f R Y F M F j C / Y D 2 l A 2 2 0 m 7 d r M J u x u h h P 4 C L x 4 U 8 e p P 8 u a / c d v m o K 0 v L D y 8 M 8 P O v E E i u D a u + + 0 U 1 t Y 3 N r e K 2 6 W d 3 b 3 9 g / L h U U v H q W L Y Z L G I V S e g G g W X 2 D T c C O w k C m k U C G w H 4 9 t Z v f 2 E S v N Y P p h J g n 5 E h 5 K H n F F j r c Z d v 1 x x q + 5 c Z B W 8 H C q Q q 9 4 v f / U G M U s j l I Y J q n X X c x P j Z 1 Q Z z g R O S 7 1 U Y 0 L Z m A 6 x a 1 H S C L W f z R e d k j P r D E g Y K / u k I X P 3 9 0 R G I 6 0 n U W A 7 I 2 p G e r k 2 M / + r d V M T X v s Z l 0 l q U L L F R 2 E q i I n J 7 G o y 4 A q Z E R M L l C l u d y V s R B V l x m Z T s i F 4 y y e v Q u u i 6 l l u X F Z q N 3 k c R T i B U z g H D 6 6 g B v d Q h y Y w Q H i G V 3 h z H p 0 X 5 9 3 5 W L Q W n H z m G P 7 I + f w B m O + M y Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " W r b 4 3 V 9 n 9 r H g + t O I p / p D Y o W t q b s = " > A A A B 6 H i c b Z B N S 8 N A E I Y n 9 a v W r 6 p H L 4 t F 8 F Q S E f R Y F M F j C / Y D 2 l A 2 2 0 m 7 d r M J u x u h h P 4 C L x 4 U 8 e p P 8 u a / c d v m o K 0 v L D y 8 M 8 P O v E E i u D a u + + 0 U 1 t Y 3 N r e K 2 6 W d 3 b 3 9 g / L h U U v H q W L Y Z L G I V S e g G g W X 2 D T c C O w k C m k U C G w H 4 9 t Z v f 2 E S v N Y P p h J g n 5 E h 5 K H n F F j r c Z d v 1 x x q + 5 c Z B W 8 H C q Q q 9 4 v f / U G M U s j l I Y J q n X X c x P j Z 1 Q Z z g R O S 7 1 U Y 0 L Z m A 6 x a 1 H S C L W f z R e d k j P r D E g Y K / u k I X P 3 9 0 R G I 6 0 n U W A 7 I 2 p G e r k 2 M / + r d V M T X v s Z l 0 l q U L L F R 2 E q i I n J 7 G o y 4 A q Z E R M L l C l u d y V s R B V l x m Z T s i F 4 y y e v Q u u i 6 l l u X F Z q N 3 k c R T i B U z g H D 6 6 g B v d Q h y Y w Q H i G V 3 h z H p 0 X 5 9 3 5 W L Q W n H z m G P 7 I + f w B m O + M y Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " W r b 4 3 V 9 n 9 r H g + t O I p / p D Y o W t q b s = " > A A A B 6 H i c b Z B N S 8 N A E I Y n 9 a v W r 6 p H L 4 t F 8 F Q S E f R Y F M F j C / Y D 2 l A 2 2 0 m 7 d r M J u x u h h P 4 C L x 4 U 8 e p P 8 u a / c d v m o K 0 v L D y 8 M 8 P O v E E i u D a u + + 0 U 1 t Y 3 N r e K 2 6 W d 3 b 3 9 g / L h U U v H q W L Y Z L G I V S e g G g W X 2 D T c C O w k C m k U C G w H 4 9 t Z v f 2 E S v N Y P p h J g n 5 E h 5 K H n F F j r c Z d v 1 x x q + 5 c Z B W 8 H C q Q q 9 4 v f / U G M U s j l I Y J q n X X c x P j Z 1 Q Z z g R O S 7 1 U Y 0 L Z m A 6 x a 1 H S C L W f z R e d k j P r D E g Y K / u k I X P 3 9 0 R G I 6 0 n U W A 7 I 2 p G e r k 2 M / + r d V M T X v s Z l 0 l q U L L F R 2 E q i I n J 7 G o y 4 A q Z E R M L l C l u d y V s R B V l x m Z T s i F 4 y y e v Q u u i 6 l l u X F Z q N 3 k c R T i B U z g H D 6 6 g B v d Q h y Y w Q H i G V 3 h z H p 0 X 5 9 3 5 W L Q W n H z m G P 7 I + f w B m O + M y Q = = < / l a t e x i t >  When h + = h − = h, namely ξ − = ξ + = ξ, the Bethe equations (3.16) contain a zero of second order which is given by the product of the two field-dependent factors:

l A a M 5 Q T i x Q p o X d l b A R 1 Z S h T a d k Q / C W T 1 6 F 1 m X V s 3 x f q 9 R v 8 j i K c A K n c A 4 e X E E d 7 q A B T W A w h G d 4 h T d H O i / O u / O x a C 0 4 + c w x / J H z + Q P q K 4 2 K < / l a t e x i t >
Let us consider a state, in the sector N = L 2 , with N − 1 real roots α 1 , . . . , α N −1 with adjacent quantum numbers n j = 1, . . . , N − 1 and a complex root α BR at and let us evaluate more precisely the deviation of this complex root with respect to the position of the double zero in the large L limit. The Bethe equation (3.16) for the complex root is Hence, using (4.40) and keeping the leading order terms in , we obtain We can now use Corollary B.1 so as to replace the sum over the real roots in (4.41) by an integral in the large L limit by means of (B.13). It leads to where F (−iζ/2 − iξ) is given by (4.23) and G(ξ) is a term of order 1 when L → +∞. We recall that F (−iζ/2 − iξ) is negative for |h| < h (1) cr (see Fig. 4, or more explicitly Fig. 9 in which this term is plotted as a function of the boundary magnetic field h), so that the deviation is exponentially decreasing in L in this regime, as expected 8 . Moreover, we see from (4.43) that there are two possible choices ± for the deviation , corresponding to the two possible choices of the sign in (4.43). Hence there are two different states with N − 1 real roots α ± 1 . . . , α ± N −1 and one boundary complex root α ± BR , that we shall denote by superscripts + or − according to the sign of the leading correction of the complex root in (4.43) (note that the + or − denomination is here not related to the left or right boundary, but only to the fact that there are two different solutions for the complex root position corresponding to the two different signs in (4.43)).
From (4.43), the boundary roots for these two states are exponentially close in L. If we denote by ξ ± the corresponding counting function, it follows from (4.20) and Appendix B.2 that Hence, using the fact that ξ + (α + j ) = ξ − (α − j ) for j = 1, . . . , N − 1, we can deduce from (4.44) that (4.45) so that the real roots of these two states are also exponentially close in L. It moreover follows from (4.44) and (4.28) that the difference of energy between these two states is also exponentially small in L: Furthermore, since from appendix A other types of states are given by solutions of the Bethe equations with at least one hole, there is a gap of energy between these two quasi-degenerate ground states and the other excited states. Let us finally remark that the exponential degeneracy at h + = h − and the gap in the spectrum are no longer present in the other regimes. Indeed, in the regimes h ∈ (−h (2) cr , −h (1) cr ) and h < −h (2) cr , it follows from our previous study that the lowest energy states contain one hole, and that their difference of energy is a direct consequence of the difference of rapidities of the hole. This is in agreement with numerical results obtained by exact diagonalization (see Fig. 3).

Form factors in the thermodynamic limit
In this section, we compute the thermodynamic limit L → ∞ of the expression for the boundary form factor of σ z 1 obtained in section 3 in two particular cases. We first consider the case of the boundary magnetization (i.e. both Bethe states coincide with each other and with the ground state) for generic values of the boundary magnetic fields. We then compute the form factor between the two quasi-degenerate ground states in the regime h + = h − = h and |h| < h (1) cr , which gives the large time limit of the boundary spin-spin autocorrelation function (1.3).

Boundary magnetization in the ground state
Let us first explain how to obtain from (3.13) the value of the boundary magnetization in the thermodynamic limit, namely the mean value σ z 1 in the ground state. This quantity has already been computed by different methods for T = 0 and h + = 0 in [16,17,23] 9 , and for finite T in [51], together with [50] where the boundary free energy was obtained for generic boundary conditions at one edge of the chain. It is relevant to see how one can derive it directly from the finite-size form factor by taking into account the precise large-L structure of the Bethe roots for the ground state that we have obtained in the previous section. We shall see in particular that, since this structure depends on both boundary fields (and therefore also on the right boundary field h + at infinity), so does the large-L limit of the boundary magnetization.
From the expressions (3.13) and (3.21), the mean value of the operator σ z 1 in an eigenstate | {λ} is in which M(λ, λ) is given by (3.22), and P(λ, λ) by (see (3.18)) where p is the derivative of the function (3.20) and .
Let us now particularise the state | {λ} in (5.1) to be the ground state of the open XXZ spin chain. We denote by α 1 , . . . , α N the corresponding Bethe roots, by ξ(µ) ≡ ξ(µ|{α}) the corresponding counting function, and set α ≡ (α 1 , . . . , α N ). From the results of section 4.2, either all N Bethe roots are real, or N − 1 of them are real whereas one of them, say α N , is an isolated complex root. We need then to compute the following trace in the thermodynamic limit: in which the vector (u(α 1 ), . . . , u(α N )) is obtained as the result of the action of the matrix M(α, α) −1 on the vector (p (α 1 ), . . . , p (α N )), i.e. is such that Let us suppose that this vector can be obtained from an odd π-periodic function u (so that in particular u(0) = u( π 2 ) = 0) which is moreover C ∞ on the real axis. Then we can use Corollary B.1 to change the sum over real roots into an integral in the left hand side of (5.5). It gives Note that, in the case of the ground state that we consider here, the set of complex roots is either empty or equal to α N , and the number n of holes is either 0 or 1. It is easy to solve (5.5) at leading order in L, by noticing that the function p can be obtained as in terms of the derivative ρ of the function (4.8), see (4.7). Therefore u solving (5.5) is of the form where u 1 (α) is a correction of order O( 1 L ) (or even of order O(L −∞ ) if the ground state does neither contain a complex root nor a hole). Note that the leading term in (5.8) is indeed an odd π-periodic meromorphic function with no pole on the real axis. Hence, combining this result with (5.4), we obtain that in which we have again replaced the sum over real roots by integrals. Note that the contributions of the complex root and/or hole vanish in the thermodynamic limit L → ∞, except for the case of the boundary root α − BR = −i(ζ/2 + ξ − + − ) for which the coefficient v(α − BR ) diverges as the inverse of the boundary root deviation − : This divergence is compensated in (5.9) by the fact that the function 2L ξ itself diverges at α − BR , via the contribution g (α − BR ), as the inverse of the boundary root deviation − : In other words, the divergence in (5.10) is compensated by a divergence of the same order in the last row of the matrix M(α, α) (3.22) The presence of the factor (1 + δ ξ − ,ξ + ) in (5.11) or in (5.12), which is equal to 1 when the two boundary fields are different and to 2 when they are equal, is due to the fact that the term g , see eq. (3.24), is summed over the two boundary fields: hence, when the latter are equal, the boundary root approaches a pole for both factors. Finally, (5.13) in which the symbol δ α N ,α − BR indicates that the last term exists only when one of the Bethe roots (and by convention the last one) coincides with the boundary root α − BR . Hence, the thermodynamic limit of the boundary magnetization in the ground state is given by where σ z 1 0 denotes the contribution given by the dense distribution of real roots, which is whereas σ z 1 BR denotes the possible contribution from the boundary root α − BR given by Here we have introduced the function H h − ,h + which is 1 when the boundary root α − BR belongs to the set of Bethe roots parametrizing the ground state, and 0 otherwise. Note that the presence of the boundary root α + BR does not play a direct role here, since it does not correspond to a divergence in the form factor. However, we have seen in section 4.2 that the presence of the boundary root α − BR in the set of roots for the ground state depends in fact on the value of both boundary magnetic fields, so that the value of the boundary magnetization depends also indirectly on the boundary field h + at infinity in the thermodynamic limit through the function H h − ,h + (see Fig. 10 and Fig. 11 for few specific evaluations and for a comparison with numerical data).
For instance, if |h + | < h cr ], and H h − ,h + = 1 otherwise. Hence the thermodynamic limit of the boundary magnetization presents, at h − = h + , a discontinuity corresponding to the boundary root contribution (5.16): 1 − q 2(2n−1) 2 1 + e 2ξ − q 2n 2 1 + e −2ξ − q 2n 2 , (5.17) Figure 10: Plot of the magnetization at the left boundary σ z 1 in the thermodynamic limit L → ∞ at ∆ = 3 as function of h − and for different values of h + . The blue dashed line shows the contribution from the bulk rapidities σ z 1 0 , while the red dashed line is the contribution from the boundary root σ z 1 BR . The sum of the two is shown as a black line, giving the boundary magnetization. Above Left: h + = 0. The discontinuity at h − = 0 is due to the boundary root moving to the other side of the chain when h − < h + (see Fig. 8) and therefore not contributing to the magnetization on the left edge for h − < h + . Above Right: same as the left plot but with h + = 1, the discontinuity being now at h − = 1. Below : same as above but with h + = −3.5. There is no discontinuity in this case since the ground state does not contain a boundary root when h + = h − . which vanishes in the limit h + → ±h (1) cr . We recall that q = e −ζ , and that the boundary fields are parametrized in this regime |h ± | < h (1) cr as h ± = sinh ζ tanhξ ± . Note that the difference between taking the limit of equal field and evaluating at exactly the same field is given by the factor 1 + δ ξ − ,ξ + in the contribution (5.16) from the boundary root. In our convention we indeed have cr ], and H h − ,h + = 1 otherwise. In that case the thermodynamic limit of the boundary magnetization is continuous at h − = h + . By symmetry of the model under the reversal of all spins and change of sign of the boundary fields, this is also the case when h + > h (1) cr . The integral in (5.15) can be computed by closing the integration contour on the lower half-plane and evaluating the corresponding residues. It gives . (5.19) In particular, at h − = h + = 0 (i.e. for ξ + = ξ − = iπ/2), we have so that as it should be. Moreover, due to the factor δ ξ − ,ξ + in the contribution (5.16) of the boundary root, we have the relation which corresponds (up to the sign) to the square of the bulk magnetization [18], as already noticed in [16].

The form factor between the two quasi-degenerate ground states
We now consider the form factor of the σ z 1 operator between the two quasi-degenerate ground states in the regime h − = h + = h (namely ξ − = ξ + = ξ) with |h| < h (1) cr : which can be expressed by means of (3.13) and (3.21). In (5.23), {α + } and {α − } denote the two sets of Bethe roots associated to the two quasi-degenerate ground states identified in section 4.3. Let us first consider the first ratio. We recall that the Bethe roots of the two states only differ by exponentially small corrections in L, 24) so that most of the prefactors in (3.21) simplify up to exponentially small corrections in L: Here we have explicitly used that the two boundary complex roots α ± N ≡ α ± BR are of the form for each row such that α ± j are real roots, i.e. for 1 ≤ j ≤ N − 1. The N -th row has to be treated separately since in that case the complex root α ± N approaches, with an exponentially small deviation ± ∼ ± , the double pole of the function g (3.24) so that the corresponding diagonal coefficient is exponentially diverging with L, see (5.12), and we have whereas the off-diagonal coefficients M(α ± , α ± ) N k with k = N remain finite (and therefore are exponentially subleading with respect to (5.28)). Finally, we obtain from (5.25), (5.27) and (5.28) that Let us now consider the second ratio in (5.23). Using again (5.24) so as to simplify the prefactors, and the fact that P is a rank-one matrix so as to decompose the determinant in the numerator, we obtain that Note that, from the orthogonality property of two different Bethe states, the first term in (5.30) should in fact vanish. Moreover, it follows from (5.24) and from (3.18) that in which we have used the notations of (5.2), and from (5.24) and (3.17) that If α ± j are real roots, i.e. for j < N , we therefore obtain that so that we recover for the first N −1 rows the elements of the Gaudin matrix (3.22) up to exponentially small corrections in L. The row j = N has to be treated separately since the two complex roots α ± N are, in the leading order, symmetrically distributed around a zero of the function a (see (5.26)). In that case, and defined the regularized function: . (5.36) Now we use again the Bethe equations for α + N , which state that Hence we obtain in which we have defined Notice that, contrary to what happens for the Gaudin matrix M(α + , α + ) in the denominator of (5.30) (see (5.12)), there is no singularity in this last row associated to the complex root. Hence, in (5.30), all terms but the one with = N vanish as (i.e. exponentially fast with L) in the large L limit due to the fact that det N M(α + , α + ) diverges as 1/ . The only term in the sum (5.30) which does not vanish is the term with = N , since the corresponding matrix elements of P(α + , α − ) themselves diverge as 1/ . Therefore (5.40) in which Note that is equal (up to the sign) to the contribution σ z 1 BR to the boundary magnetization from the boundary root when ξ − = ξ + = ξ, see eq. (5.16). Finally, which is exactly half of the discontinuity of the boundary magnetization at h + = h − , see (5.17). In Fig. 12 and Fig. 13 we report this result at h = 0 and h = 1 and as function of ∆ compared to the values of the form factors in a finite size chain.

Conclusion
In this paper we have shown that the physics of the open XXZ chain in the massive regime ∆ > 1 is strongly influenced by the presence of its boundary modes. In the language of Bethe ansatz these modes correspond to isolated complex Bethe roots converging exponentially fast with L towards a zero of the boundary factor, and can be understood as excitations that are exponentially pinned at one of the two edges of the chain. We have shown that for chains of even size there exists a regime, |h ± | < h (1) cr , in which such a boundary root is present both in the Bethe solutions for the ground state and for the first excited state. The values of the two boundary magnetic fields determine at which edge the corresponding boundary excitation is localized in the ground state. As a consequence of this localization, the value of the spin magnetization at one of the boundaries of the chain also depends indirectly on the value of the magnetic field at the opposite edge, even in the thermodynamic limit: it presents in particular a discontinuity in this limit at h − = h + . Moreover we have shown that, when the two boundary fields are equal cr = ∆ − 1), the spectrum is gapped in the thermodynamic limit and the ground state is doubly degenerate up to exponentially small corrections in L. In this case the spin-spin autocorrelation function on one of the two edges relaxes at large time to a finite value, given by the contribution of the boundary root to the boundary magnetization, or equivalently by the discontinuity of this boundary magnetization at h − = h + .
It would be desirable to establish a more direct relation between the boundary root and the strong zero mode found in [9]. In particular few questions can be immediately formulated. Is the strong zero mode related to the creation operator of the boundary excitation corresponding to the boundary root at one edge of the chain? And what happens at finite temperatures? Does the degeneracy observed in [9] at all energies in the spectrum, which is due to the strong zero mode, translate into the presence of a boundary root also for finite temperature states? Are there two degenerate representative thermal states (via the thermodynamic Bethe ansatz [52]) distinguished by two opposite deviations of the boundary root, as it is the case for the ground state? We postpone the study of these interesting questions to future works.
Finally, it would also be interesting to investigate supersymmetric properties of the open chain (see for example [53,54]), and to understand how to formulate an ensemble of states in the presence of boundary roots. Namely how to determine for example the steady state after a quantum quench [55] (Generalised Gibbs Ensemble) close to the edge of the chain. It is indeed evident that the value of the conserved charges that are extensive in the system must be supplied with the information about the boundary [56] which at the moment is not clear how to include in the steady state of an interacting system.

A The Bethe equations in logarithmic form and allowed quantum numbers for the real roots
In this appendix we consider the Ising limit ∆ → +∞, i.e. ζ → +∞, of the logarithmic Bethe equations (4.11). More precisely, we suppose that ∆, and therefore ζ, are large but finite, and we write the logarithmic Bethe equations at leading order in ζ. At leading order in ζ the counting function (4.12) becomes linear, which enables us to determine the allowed quantum numbers n j for the real roots. It is easy to see from the expressions (4.13)-(4.16) that, for α ∈ R, and that where, for σ ∈ {+, −}, Let us now consider a solution {λ} ≡ {λ 1 , . . . , λ N } of the logarithmic Bethe equations (4.11). Let n w be the number of wide roots λ k such that ζ = o(| (λ k )|). Then, if α ∈ R, so that the logarithmic Bethe equations (4.11) for each real root λ j become, at leading order in ζ: Since the allowed real solutions are such that 0 < λ j < π 2 (we recall that we have to discard the obvious solutions 0 and π 2 , see footnote 3), the integers n j associated to real roots can then take only the possible values cr ]).
1. They are L − N + n w − 1 possible quantum numbers for the real roots.
2. The maximum number of real Bethe roots for a solution in the sector N = L 2 is N − 1. Such a solution therefore contains an additional isolated complex root, which may correspond either to one of the two possible boundary roots α σ BR (4.9) with σ ∈ {+, −}, or to a wide root. The corresponding quantum numbers n j (j = 1, . . . N − 1) for the real roots are such that   Other types of solutions in the sector N = L 2 contain more holes, except the solution with N − 2 real roots and a pair of bulk close roots (i.e. from [46] a 2-string), which has to be compared with the solution 2b since it also contains one hole.
3. In the sector N = L 2 − 1, there exists a solution with N real roots (and therefore no complex root) with quantum numbers to be distributed within the set {1, . . . , N + 1}. Hence this solution contains a hole at some position h ∈ {1, . . . , N + 1}. Other possible solutions in that sector or in sectors N < L 2 − 1 contain more holes.
Case B. One of the fields is in the interval delimited by the two critical fields h (1) cr and h (2) cr and the other is not.
1. They are L − N + n w possible quantum numbers for the real roots.
2. The maximum number of real Bethe roots for a solution in the sector N = L 2 is N . It corresponds to a full set of adjacent quantum numbers j = 1, . . . , N (no hole and no complex root). Other types of solutions with complex roots in that sector contain one or more hole(s).
3. Solutions in sectors N < L 2 contain at least two holes.
Case C. Both boundary fields are in the interval delimited by the two critical fields h (1) cr and h (2) cr (h (1) cr < h ± < h (2) cr ).
1. They are L − N + n w + 1 possible quantum numbers for the real roots.
2. In the sector N = L 2 , there exists a solution with N real roots (and therefore no complex root) with quantum numbers to be distributed within the set {1, . . . , N + 1}. Hence this solution contains a hole at some position h ∈ {1, . . . , N + 1}. Other possible solutions in that sector, i.e. with some complex roots, contain two or more holes.

B Controlling the finite-size corrections in the large L limit
In this appendix, we explain how to control the finite-size corrections to the integral over the density which come from sums over real Bethe roots in the large L limit. Let {λ} ≡ {λ 1 , . . . , λ N } be a solution of the Bethe equations (4.4). We suppose that this solution corresponds to an infinite number of real roots (i.e. of order L), with a finite number of complex roots and a finite number of holes in the thermodynamic limit. The logarithmic equation for the real roots can be written as in (4.18), in terms of the positions h 1 , . . . , h n of the holes in the adjacent set of quantum numbers for the real roots, with M given by (A.7) that we suppose to be of the same order as L. Note that the counting function ξ(α) ≡ ξ(α|{λ}) associated to this set of Bethe roots, which is defined as in (4.12), satisfies the following properties for α ∈ R: where O(L −∞ ) stand for exponentially small corrections in L.
Proof. The proof can be done with similar arguments as in [19] (see also [47]), adapted here to the case of the open chain and of general low-energy states.
Since f is π-periodic, and since f is even, We now make a change of variables using the function ξ defined from the counting function ξ as which is still odd and invertible and satisfies, instead of (B.2) and (B.3), the properties ξ(α + π) = ξ(α) + π, ξ(0) = 0, ξ π 2 = − ξ − π 2 = π 2 . (B.10) Hence, the function f • ξ −1 is also even and π-periodic, so that we have Multiplying by M/L and setting appart the contributions of the holes from the ones of the real roots we obtain (B.6). Let g be a C ∞ -function such that g is π-periodic. Then Proof. (B.12) is a direct consequence of (B.6). If g (x) is π-periodic then g(x) − c g x is also π-periodic, where c g = 1 π π 2 − π 2 g (x) dx = g( π 2 ) − g(− π 2 ) π = g(y + π) − g(y) π , ∀y. (B.14) Hence one can apply (B.12) to g(x) − c g x, 15) and the second integral vanishes due to the fact that ξ is an even function.

B.2 Finite size corrections to the counting function
We can in particular apply (B.13) to transform the sum over real roots in the definition (4.12) of the counting function: in which Z is the set of indices corresponding to the complex roots (i.e. (λ k ) = 0 if k ∈ Z). Deriving (B.16), we obtain the following integral equation for ξ : Hence, the expression (B.16) of the counting function can be decomposed in terms of the different contributions of the real roots, the complex roots and the holes as in (4.20). In (4.20), ξ 0 (α) is the common contribution of the "Fermi sea" of real roots. It is an odd function, and its derivative is defined as the solution of the integral equation . where ρ is the density (4.8) solution of (4.7), and where ξ open is the correction due to the 1/L terms in (B.18), which is defined as the solution to the integral equation . (B.20) The function ξ µ , which corresponds to the contribution to the counting function of an excitation (an additional complex root or a hole at position µ) with respect to the above Fermi sea of real roots, is also an odd function with derivative being the solution of the integral equation: This latter can easily be computed in Fourier modes by using the following lemma: where H denotes the Heaviside function and sgn the sign function.
In particular, the dressed energy of a hole with rapidityλ h ∈ (0, π 2 ) is given by (4.32), whereas the dressed energy of the boundary root (4.9) is given by (4.33).