<figure><img src="image_1.png" /></figure>

HAL
open science

# Boundary overlap in the open XXZ spin chain

Charbel Abetian, Nikolai Kitanine, Véronique Terras

► To cite this version:

Charbel Abetian, Nikolai Kitanine, Véronique Terras. Boundary overlap in the open XXZ spin chain. SciPost Physics, 2025, 18, pp.026. 10.21468/SciPostPhys.18.1.026 . hal-04974743

HAL Id: hal-04974743

https://hal.science/hal-04974743v1

Submitted on 21 Nov 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Boundary overlap in the open XXZ spin chain

Charbel Abetian¹, Nikolai Kitanine²* and Veronique Terras¹

1 Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France

2 Institut de Mathématiques de Bourgogne, UMR 5584, CNRS and Université Bourgogne Europe, F-21000 Dijon, France

★ Nikolai.Kitanine@u-bourgogne.fr

## Abstract

In this paper we compute the overlaps of the ground states for the open spin chains after a change of one of the boundary magnetic fields. It can be considered as the first step toward the study of the boundary quench problem: behaviour of an open spin chain after an abrupt change of one boundary magnetic field.

<figure><img src="image_3.png" /></figure>

Copyright C. Abetian *et al.*

This work is licensed under the Creative Commons

Attribution 4.0 International License.

Published by the SciPost Foundation.

Received 2024-10-02

Accepted 2024-12-04

Published 2025-01-21

doi:10.21468/SciPostPhys.18.1.026

<figure><img src="image_4.png" /></figure>

Check for
updates

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>The open XXZ spin chain with longitudinal boundary fields</td><td>3</td></tr><tr><td>2.1</td><td>Solution of the model by algebraic Bethe Ansatz</td><td>3</td></tr><tr><td>2.2</td><td>Scalar products of Bethe states</td><td>6</td></tr><tr><td>2.3</td><td>Description of the ground state</td><td>6</td></tr><tr><td>2.4</td><td>Changing the boundary field h</td><td>9</td></tr><tr><td>3</td><td>Cauchy determinant representation for the overlap</td><td>9</td></tr><tr><td>3.1</td><td>Determinant representation</td><td>9</td></tr><tr><td>3.2</td><td>Gaudin extraction</td><td>10</td></tr><tr><td>3.3</td><td>Product formula for the overlap</td><td>13</td></tr><tr><td>4</td><td>Taking the thermodynamic limit</td><td>14</td></tr><tr><td>4.1</td><td>The case of real Bethe roots</td><td>14</td></tr><tr><td>4.1.1</td><td>General strategy</td><td>14</td></tr><tr><td>4.1.2</td><td>Computation of the overlap</td><td>16</td></tr><tr><td>4.2</td><td>Cases with presence of a boundary root</td><td>18</td></tr><tr><td>4.2.1</td><td>Presence of two boundary roots λ<sup>σ1</sup><sub>BR</sub> and μ<sup>σ2</sup><sub>BR</sub></td><td>18</td></tr><tr><td>4.2.2</td><td>Presence of only one boundary root λ<sup>σ1</sup><sub>BR</sub></td><td>20</td></tr><tr><td>5</td><td>The overlap in the thermodynamic limit: Summary of the results</td><td>22</td></tr><tr><td>5.1</td><td>The overlap for a chain with an odd number of sites</td><td>22</td></tr><tr><td>5.2</td><td>The overlap for a chain with an even number of sites</td><td>24</td></tr></table>
<table><tr><td>6 Conclusion</td><td>26</td></tr></table>

<table><tr><td>A Convergence of numerical results</td><td>26</td></tr></table>

<table><tr><td>References</td><td>29</td></tr></table>

# 1 Introduction

Spin chains with open boundaries [1] are interesting example of integrable systems interacting with their environment: the interaction is in that case implemented through the boundary magnetic fields acting on the first and the last site of the chain. It is therefore important to understand what happens if we change one of these boundary magnetic fields. In particular, it can be interesting to study the dynamics after an abrupt change of the magnetic field (boundary quench), or to consider a time-dependent magnetic field (boundary driven system). Although it is a priori only a local perturbation, it was observed that a change of the boundary magnetic field can induce a macroscopic change of the ground state [2] or a phase transition [3].

Integrability gives a unique possibility to study such systems in details. However, a necessary condition to deal with such kind of quench problems in the framework of quantum integrability is the possibility to compute overlaps [4] between the eigenstates before and after the quench. In this paper, we propose the first step in this direction: we compute the overlaps between the ground states of the chain for different values of one the boundary magnetic field.

Typically, being able to exactly calculate the relevant overlaps is crucial to making the study of quench dynamics possible [5–7]. This is often a complicated problem, and the main difficulty comes in general from the fact that the eigenstates before and after the quench are constructed using different algebras: this is for instance the case if the quench comes from a change of the coupling constant. For a boundary quench in which only one boundary magnetic field is modified, however, we can use the same reflection algebra [8] to describe the eigenstates before and after the quench. This gives a possibility to express overlaps by means of the open chain version of Slavnov formula for the scalar products [9–12].

We consider in this paper the XXZ spin chain in the massive antiferromagnetic regime with diagonal boundary terms (i.e. boundary magnetic fields parallel to the z-anisotropy direction). It was shown [2, 13] that low energy properties of the model in this case are determined by two boundary magnetic fields and changing one of them can change the ground state in a macroscopic way. In the diagonal boundary terms case, the eigenstates can be obtained from the boundary version of the Algebraic Bethe Ansatz (ABA) [1], and the Slavnov formula for the overlaps leads to a determinant expression in which the size of the matrix is proportional to the length L of the spin chain [12].

To compute the thermodynamic limit of this determinant representation for the overlaps, we use the method first introduced in [14] which we call *Gaudin extraction*: it permits to express the ratio of Slavnov and Gaudin determinants as a Cauchy determinant (up to an exponentially small correction in L). The resulting Cauchy determinant can be easily computed as a product over the ground state Bethe roots.

Using the ground state configuration of the Bethe roots [2, 15], we obtain in the thermodynamic limit a very simple expression for the overlap in terms of the anisotropy parameter and of the two values of the boundary magnetic field. We consider here two possible configurations of Bethe roots: (1) the ground state Bethe roots are all real or (2) there is one complex Bethe root (boundary string). In both configurations there are no real holes, the solutions with holes
(corresponding to excited states) will be treated elsewhere. Our final result is valid in the half-infinite chain limit, however we are able to control the finite size corrections to this result, and to show that all corrections are exponentially small in $L$. We also observe numerically a very rapid convergence towards our analytic result (see appendix A).

We would like to mention that similar results [16, 17] in the massive case can be obtained from the $q$-vertex operator approach [18–20]. It is noteworthy that, in general, the Gaudin extraction method introduced in [14] permits to reproduce in the ABA framework the results for the form factors obtained from the $q$-vertex operator approach. In particular, it permits to relate two possible descriptions for the excited states.

The paper is organised as follows. In section 2, we give a brief technical introduction to the computation of boundary overlaps in the ABA framework: we remind the construction of eigenstates, the structure of the ground state Bethe roots according to the different values of the boundary magnetic fields, and the Slavnov determinant formula for the scalar products. In section 3, we express the overlap as a ratio of determinants and show how to apply to this ratio the Gaudin extraction technique leading to a Cauchy determinant representation. Section 4 is devoted to the computation of the overlaps in the thermodynamic limit: we study separately the case where all Bethe roots are real and the case of a presence of a boundary complex root. Our final result, concerning the thermodynamic limit of the overlap between ground states before and after a boundary quench, is presented in section 5. In the Appendix A we illustrate the rapid convergence of numerical results to our analytic formula for the overlap.

## 2 The open XXZ spin chain with longitudinal boundary fields

We consider here the open XXZ spin-1/2 chain with longitudinal boundary fields. The Hamiltonian of this model is

$$ \mathbf{H} = \mathbf{H}_{h^{-},h^{+}} = \sum_{i=1}^{L-1} [\sigma_{i}^{x} \sigma_{i+1}^{x} + \sigma_{i}^{y} \sigma_{i+1}^{y} + \Delta(\sigma_{i}^{z} \sigma_{i+1}^{z} - 1)] + h^{-} \sigma_{1}^{z} + h^{+} \sigma_{N}^{z}, \quad (1) $$

in which $\sigma_i^\alpha$, $\alpha = x, y, z$, denote the Pauli matrices at site $i$, $\Delta$ is the anisotropy of the coupling along the $z$-direction, and $h_-, h_+$ are the values of the boundary fields. We restrict our study to the case in which these boundary fields are oriented along the direction $z$ of the anisotropy, and to the antiferromagnetic regime for which $\Delta > 1$. In the following, we shall parametrize the anisotropy parameter $\Delta$ as

$$ \Delta = \frac{q+q^{-1}}{2} = \cosh \zeta, \quad q = e^{-\zeta}, \quad \zeta > 0, \quad (2) $$

and the boundary fields $h^\pm$ as

$$ h^{\pm} = -\sinh \zeta \coth \xi^{\pm}. \quad (3) $$

### 2.1 Solution of the model by algebraic Bethe Ansatz

This model is integrable. It can be solved by a boundary version of the algebraic Bethe Ansatz [1], from the construction of a boundary monodromy matrix $U(\lambda)$. The latter is a 2 × 2 matrix of quantum operators depending on a spectral parameter $\lambda$:

$$ U(\lambda) = \begin{pmatrix} A(\lambda) & B(\lambda) \\ C(\lambda) & D(\lambda) \end{pmatrix}, \quad A(\lambda), B(\lambda), C(\lambda), D(\lambda) \in \text{End} \mathbb{H}, \quad (4) $$
where $\mathbb{H} = \otimes_{n=1}^{L} \mathbb{H}_n$, $\mathbb{H}_n = \mathbb{C}^2$, is the $2^L$-dimensional quantum space of states of the model. It satisfies, with the 6-vertex R-matrix,

$$R(\lambda) = \begin{pmatrix} \sin(\lambda - i\zeta) & 0 & 0 & 0 \\ 0 & \sin(\lambda) & \sin(-i\zeta) & 0 \\ 0 & \sin(-i\zeta) & \sin(\lambda) & 0 \\ 0 & 0 & 0 & \sin(\lambda - i\zeta) \end{pmatrix}, \qquad (5)$$

the reflection equation introduced in [21]:

$$R(\lambda - \mu) (\mathcal{U}^t(-\lambda) \otimes \mathrm{Id}_2) R(\lambda + \mu + i\zeta) (\mathrm{Id}_2 \otimes \mathcal{U}^t(-\mu)) \\ = (\mathrm{Id}_2 \otimes \mathcal{U}^t(-\mu)) R(\lambda + \mu + i\zeta) (\mathcal{U}^t(-\lambda) \otimes \mathrm{Id}_2) R(\lambda - \mu). \quad (6)$$

Here, $\mathrm{Id}_2$ stands for the 2 × 2 identity matrix, and $\mathcal{U}^t$ should be understood as the transposed of $\mathcal{U}$ as a 2 × 2 matrix. The reflection equation (6) gives the commutation relations for the operators $\mathcal{A}(\lambda)$, $\mathcal{B}(\lambda)$, $\mathcal{C}(\lambda)$, $\mathcal{D}(\lambda)$ elements of the matrix (4).

The boundary monodromy matrix $\mathcal{U}(\lambda)$ can be constructed from simpler scalar solutions of (6), the boundary K-matrices. Let us introduce the following diagonal scalar solution of (6),

$$K(\lambda; \xi) = \begin{pmatrix} \sin(\lambda + i\zeta/2 + i\xi) & 0 \\ 0 & \sin(i\xi - \lambda - i\zeta/2) \end{pmatrix}, \qquad (7)$$

and define

$$K^{-}(\lambda) = K(\lambda; \xi^{-}), \quad K^{+}(\lambda) = K(\lambda - i\zeta; \xi^{+}), \qquad (8)$$

where $\xi^{\pm}$ parametrize the boundary fields at both ends of the chain, see (3). Let us also introduce the bulk monodromy matrix as the following product of R-matrices along the chain:

$$T(\lambda) \equiv T_0(\lambda) = R_{0L}(\lambda - \zeta_L) \dots R_{01}(\lambda - \zeta_1). \qquad (9)$$

Here, the notation $R_{0n}$, for $1 \le n \le L$, means that the corresponding R-matrix acts on $\mathbb{V}_0 \otimes \mathbb{H}_n$, where $\mathbb{H}_n$ is the quantum space at site $n$, whereas $\mathbb{V}_0 = \mathbb{C}^2$ is a 2-dimensional auxiliary space, so that $T(\lambda)$ can be considered as 2 × 2 matrix (on the auxiliary space $\mathbb{V}_0$) of quantum operators. The parameters $\zeta_1, \dots, \zeta_L$ are inhomogeneity parameters which may be introduced for convenience. Then the boundary monodromy matrix can be constructed from $K^+$ and $T$ as

$$\mathcal{U}^t(\lambda) \equiv \mathcal{U}_0^{t_0}(\lambda) = (-1)^L T_0^{t_0}(\lambda) [K_0^+]^{t_0}(\lambda) [\sigma_0^\gamma T_0^{t_0}(-\lambda) \sigma_0^\gamma]^{t_0}. \qquad (10)$$

Taking the trace, on the auxiliary space $\mathbb{V}_0$, of the product of this monodromy matrix with the other boundary matrix $K^-$, one obtains a one-parameter family of commuting transfer matrices, which are generating functions of the commuting conserved charges of the model:

$$\mathcal{T}(\lambda) = \mathrm{tr}_0 [K_0^-(\lambda) \mathcal{U}_0(\lambda)]. \qquad (11)$$

The Hamiltonian (1) of the open spin-1/2 chain can then be expressed in terms of this transfer matrix in the homogeneous limit $\zeta_\ell = -i\zeta/2$, $\ell = 1, \dots, L$, as

$$\mathbf{H} = \frac{-i \sinh \zeta}{\mathcal{T}(\lambda)} \frac{d}{d\lambda} \mathcal{T}(\lambda)|_{\lambda = -i\zeta/2} + \frac{1}{\cosh \zeta} - 2L \cosh \zeta. \qquad (12)$$

In this algebraic framework, the common eigenstates of the transfer matrices (and of the Hamiltonian in the homogeneous limit) can be constructed as Bethe states, i.e. by using the operators $\mathcal{B}(\lambda)$ and $\mathcal{C}(\lambda)$ of (4) as generalised creation and annihilation operators when acting
on a reference state. Let us define $|0\rangle$ (respectively $\langle 0|$) to be the reference state (respectively the dual reference state) corresponding to the ferromagnetic state where all spins are up, and let us consider, for a given set of spectral parameters $\{\lambda\} \equiv \{\lambda_1, \dots, \lambda_N\}$, Bethe states of the form

$$|\{\lambda\}\rangle = \prod_{j=1}^{N} \mathcal{B}(\lambda_j) |0\rangle, \quad \langle\{\lambda\}| = \langle 0| \prod_{j=1}^{N} \mathcal{C}(\lambda_j). \qquad (13)$$

Then, the Bethe states (13) are eigenstates of the transfer matrix $\mathcal{T}(\nu)$ (11), with eigenvalue

$$\tau(\nu, \{\lambda\}) = \mathbf{a}(\nu) \frac{\sin(2\nu - i\zeta)}{\sin(2\nu)} \prod_{i=1}^{N} \frac{\mathfrak{s}(\nu + i\zeta, \lambda_i)}{\mathfrak{s}(\nu, \lambda_i)} + \mathbf{a}(-\nu) \frac{\sin(2\nu + i\zeta)}{\sin(2\nu)} \prod_{i=1}^{N} \frac{\mathfrak{s}(\nu - i\zeta, \lambda_i)}{\mathfrak{s}(\nu, \lambda_i)}, \quad (14)$$

if the corresponding set of parameters $\{\lambda\}$ satisfies the system of Bethe equations

$$\mathbf{a}(\lambda_j) \frac{\sin(2\lambda_j - i\zeta)}{\sin(2\lambda_j)} \prod_{k=1}^{N} \mathfrak{s}(\lambda_j + i\zeta, \lambda_k) + \mathbf{a}(-\lambda_j) \frac{\sin(2\lambda_j + i\zeta)}{\sin(2\lambda_j)} \prod_{k=1}^{N} \mathfrak{s}(\lambda_j - i\zeta, \lambda_k) = 0, \quad (15)$$

for $j = 1, \dots, N$. Here we have used the shortcut notation

$$\mathfrak{s}(\lambda, \mu) = \sin(\lambda + \mu) \sin(\lambda - \mu) = \sin^2 \lambda - \sin^2 \mu, \qquad (16)$$

and defined

$$\mathbf{a}(\lambda) = (-1)^L a(\lambda) d(-\lambda) \sin(\lambda + i\xi^+ + i\zeta/2) \sin(\lambda + i\xi^- + i\zeta/2), \quad (17)$$

with

$$a(\lambda) = \prod_{\ell=1}^{L} \sin(\lambda - \zeta_{\ell} - i\zeta), \quad d(\lambda) = \prod_{\ell=1}^{L} \sin(\lambda - \zeta_{\ell}). \qquad (18)$$

Bethe states of the form (13) are called on-shell Bethe states if the corresponding set of parameters $\{\lambda\}$ satisfies the Bethe equations (15), and off-shell Bethe states otherwise. On-shell Bethe states are therefore eigenstates of the transfer matrix (11). They are also eigenstates of the Hamiltonian (1) in the homogeneous limit $\zeta_{\ell} = -i\zeta/2$, $\ell = 1, \dots, L$. The corresponding energy is

$$E(\{\lambda\}) = h^+ + h^- - \sum_{j=1}^{N} \frac{4 \sinh^2 \zeta}{\cosh \zeta - \cos(2\lambda_j)}. \qquad (19)$$

In the following, it will be convenient to rewrite the Bethe equations (15) as<fnref n="1" />

$$a(\lambda_j | \{\lambda\}) = 1, \quad 1 \le j \le N, \qquad (20)$$

by introducing the exponential counting function

$$a(\nu|\{\lambda\}) = \frac{\mathbf{a}(\nu) \sin(i\zeta - 2\nu)}{\mathbf{a}(-\nu) \sin(i\zeta + 2\nu)} \prod_{\ell=1}^{N} \frac{\mathfrak{s}(\nu + i\zeta, \lambda_\ell)}{\mathfrak{s}(\nu - i\zeta, \lambda_\ell)}. \qquad (21)$$

In terms of this exponential counting function, the eigenvalue (14) of the transfer matrix can be expressed as

$$\tau(\nu, \{\lambda\}) = -\mathbf{a}(-\nu) \frac{\sin(2\nu + i\zeta)}{\sin(2\nu)} \prod_{i=1}^{N} \frac{\mathfrak{s}(\nu - i\zeta, \lambda_i)}{\mathfrak{s}(\nu, \lambda_i)} [a(\nu|\{\lambda\}) - 1]. \qquad (22)$$

<footnote n="1">¹In this rewriting, one should exclude the root 0, $\pi/2$ of the function $a(u|\{\lambda\}) - 1$.</footnote>
## 2.2 Scalar products of Bethe states

In this algebraic framework, one can express scalar products between an on-shell and an off-shell Bethe vectors of the form (13) in terms of a determinant of usual functions, see [11, 12]. This representation is the analog, in the open case, of the determinant representation obtained in [9] in the periodic case. We recall it here.

For {$\lambda$} ≡ {$\lambda_1, ..., \lambda_N$} a solution of the Bethe equations and {$\mu$} ≡ {$\mu_1, ..., \mu_N$} an arbitrary set of parameters, the scalar product $\langle\{\lambda\}|\{\mu\}\rangle$ is equal to the scalar product $\langle\{\mu\}|\{\lambda\}\rangle$, and admits the following determinant representation:

$$ \begin{aligned} \langle\{\lambda\}|\{\mu\}\rangle &= \langle\{\mu\}|\{\lambda\}\rangle \\ &= \prod_{j=1}^{N} \left[ (-1)^L a(\lambda_j) d(-\lambda_j) \frac{\sin(2\lambda_j - i\zeta) \sin(2\mu_j - i\zeta)}{\sin(2\mu_j)} \frac{\sin(\lambda_j + i\xi_+ + i\frac{\zeta}{2})}{\sin(\lambda_j - i\xi_- - i\frac{\zeta}{2})} \right] \\ & \quad \times \prod_{j<k} \left[ \frac{\sin(\lambda_j + \lambda_k - i\zeta)}{\sin(\lambda_j + \lambda_k + i\zeta)} \frac{1}{s(\lambda_j, \lambda_k)s(\mu_k, \mu_j)} \right] \det_N[H(\lambda, \mu)], \end{aligned} \quad (23) $$

where the elements of the $N \times N$ matrix $H(\lambda, \mu)$ are

$$ [H(\lambda, \mu)]_{jk} = \frac{\sin(-i\zeta)}{s(\mu_k, \lambda_j)} \left[ a(\mu_k) \prod_{\ell \neq j} s(\mu_k + i\zeta, \lambda_\ell) - a(-\mu_k) \prod_{\ell \neq j} s(\mu_k - i\zeta, \lambda_\ell) \right], \quad (24) $$

for $\lambda \equiv (\lambda_1, ..., \lambda_N)$ and $\mu \equiv (\mu_1, ..., \mu_N)$. When {$\lambda$} and {$\mu$} coincide, (23) becomes

$$ \begin{aligned} \langle\{\lambda\}|\{\mu\}\rangle &= \prod_{j=1}^{N} \left[ (-1)^L a(\lambda_j) d(-\lambda_j) \frac{\sin(\lambda_j + i\xi_+ + i\frac{\zeta}{2})}{\sin(\lambda_j - i\xi_- - i\frac{\zeta}{2})} \right] \\ & \quad \times \prod_{j<k} \frac{\sin(\lambda_j + \lambda_k - i\zeta)}{\sin(\lambda_j + \lambda_k + i\zeta)} \prod_{k=1}^{N} \frac{a(-\lambda_k) \prod_{\ell=1}^{N} s(\lambda_k - i\zeta, \lambda_\ell)}{i \sin^2(2\lambda_k) \prod_{\ell \neq k} s(\lambda_k, \lambda_\ell)} \det_N[\mathcal{M}(\lambda, \lambda)], \end{aligned} \quad (25) $$

where the elements of the $N \times N$ matrix $\mathcal{M}(\lambda, \lambda)$ are

$$ [\mathcal{M}(\lambda, \lambda)]_{jk} = i \delta_{jk} a'(\lambda_j | \{\lambda\}) - 2\pi [K(\lambda_j - \lambda_k) - K(\lambda_j + \lambda_k)]. \quad (26) $$

Here we have defined

$$ K(\lambda) = \frac{\sinh(2\zeta)}{2\pi \sin(\lambda + i\zeta) \sin(\lambda - i\zeta)} = \frac{1}{2\pi} [t(\lambda) + t(-\lambda)], \quad (27) $$

with

$$ t(\nu) = \frac{\sinh \zeta}{\sin \nu \sin(\nu - i\zeta)}, \quad (28) $$

and $a'(\nu|\{\lambda\})$ denotes the derivative, with respect to the variable $\nu$, of the exponential counting function (21).

## 2.3 Description of the ground state

The study of the solutions of Bethe equations in the homogeneous limit $\zeta_\ell = -i\frac{\xi_\ell}{2}$, $\ell = 1, ..., L$, and in particular of the ground state of the Hamiltonian (1) in the thermodynamic limit $L \to \infty$, has been performed in [2, 15, 22]. We briefly recall here the results which will be useful for our study. As mentioned above, we restrict ourselves to the antiferromagnetic regime
$\Delta > 1$, with $\Delta$ parametrized as in (2). The boundary fields $h^\sigma$, $\sigma = \pm$, are parameterized as in (3), with

$$ \xi^{\sigma} = -\tilde{\xi}^{\sigma} + i\delta^{\sigma} \frac{\pi}{2}, \quad \tilde{\xi}^{\sigma} \in \mathbb{R}, \quad \delta^{\sigma} = \begin{cases} 1, & \text{if } |h^{\sigma}| < \sinh \zeta, \\ 0, & \text{if } |h^{\sigma}| > \sinh \zeta. \end{cases} \qquad (29) $$

Due to the parity and periodicity of the Bethe equations, it is enough to consider solutions of (20) such that $0 < \Re(\lambda_j) < \frac{\pi}{2}$ or $(\Re(\lambda_j) = 0, \frac{\pi}{2}$ and $\Im(\lambda_j) < 0)$. Low-energy states are given by solutions of (20) close to half-filling for which nearly all Bethe roots are real. Complex roots appear by complex conjugated pairs $\lambda_j, \bar{\lambda}_j$, except if $\Re(\lambda_j) = 0, \frac{\pi}{2}$.

The Bethe equations for real roots can be rewritten in logarithmic form as

$$ \hat{\xi}(\lambda_j|\{\lambda\}) = \frac{\pi n_j}{L}, \quad n_j \in \mathbb{Z}, \qquad (30) $$

in which $\hat{\xi}(\mu|\{\lambda\}) = -\frac{i}{2L} \log a(\mu|\{\lambda\})$ is the counting function:

$$ \hat{\xi}(\mu|\{\lambda\}) = p(\mu) + \frac{g(\mu)}{2L} - \frac{\theta(2\mu)}{2L} + \frac{1}{2L} \sum_{k=1}^{N} [\theta(\mu - \lambda_k) + \theta(\mu + \lambda_k)]. \qquad (31) $$

The functions $p, \theta$ and $g$ appearing in (31) are defined on the real axis by

$$ p(\lambda) = i \log \frac{\sin(i\zeta/2 + \lambda)}{\sin(i\zeta/2 - \lambda)}, \quad \theta(\lambda) = i \log \frac{\sin(i\zeta - \lambda)}{\sin(i\zeta + \lambda)}, \qquad (32) $$

$$ g(\lambda) = i \log \prod_{\sigma=\pm} \frac{\sin(\lambda - i\xi^{\sigma} - i\zeta/2)}{\sin(\lambda + i\xi^{\sigma} + i\zeta/2)}, \qquad (33) $$

and by appropriate analytical continuation around this axis.

In the thermodynamic limit, the real Bethe roots for the ground state (and more generally for low-energy states) form a dense distribution on the interval $(0, \frac{\pi}{2})$, which can be extended by parity on the interval $(-\frac{\pi}{2}, \frac{\pi}{2})$. The corresponding Bethe equations (30) turn at the leading order into an integral equation for their density $\rho(\lambda)$,

$$ \rho(\lambda) + \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} K(\lambda - \beta) \rho(\beta) d\beta = \frac{p'(\lambda)}{\pi}, \qquad (34) $$

with solution

$$ \rho(\lambda) = \frac{\vartheta'_{1}(0) \vartheta_{3}(\lambda)}{\pi \vartheta_{2}(0) \vartheta_{4}(\lambda)}. \qquad (35) $$

Here and in the following, $\vartheta_i(\lambda) \equiv \vartheta_i(\lambda, q)$, $i \in \{1, 2, 3, 4\}$, denote the Theta functions of nome $q$ defined as in [23]. The function $p'$ in (34) is the derivative of the function $p$ (32) and corresponds, up to a shift in the argument, to the function $t$ (28), whereas the function $K$ (27) corresponds, up to a scalar factor, to the derivative of the function $\theta$:

$$ p'(\lambda) = t(\lambda + i\frac{\zeta}{2}), \quad K(\lambda) = -\frac{1}{2\pi}\theta'(\lambda). \qquad (36) $$

A detailed analytical description of the ground state was performed in [2], and different regimes can be distinguished according to the values of $h^+$ and $h^-$ with respect to the two boundary critical fields $h_{cr}^{(1)}$ and $h_{cr}^{(2)}$ defined as [15, 19]

$$ h_{cr}^{(1)} = \Delta - 1, \quad h_{cr}^{(2)} = \Delta + 1. \qquad (37) $$

In particular, for

$$ h^{-}, h^{+} > h_{cr}^{(1)}, \quad \text{or} \quad -h^{-}, -h^{+} > h_{cr}^{(2)}, \quad \text{if } L \text{ is even,} \qquad (38) $$
or for

$$h^{-}, -h^{+} > h_{cr}^{(1)}, \text{ or } -h^{-}, h^{+} > h_{cr}^{(1)}, \text{ if } L \text{ is odd,} \quad (39)$$

the spectrum becomes gapless in the thermodynamic limit. In the following, we shall restrict our study to values of the boundary fields for which the spectrum remains gapped in the thermodynamic limit, which means that we exclude the consideration of the configurations (38) and (39). We therefore have to distinguish the following different cases [2]:

### For L even:

$$(A) \quad |h^{\sigma_1}| < h_{cr}^{(1)} \text{ with } h^{\sigma_1} > h^{\sigma_2} (\{\sigma_1, \sigma_2\} = \{+,-\}).$$

The ground state is in the sector $N = \frac{L}{2}$ (with magnetization 0). It is given by $\frac{L}{2} - 1$ real roots with adjacent quantum numbers $n_j = 1, \dots, \frac{L}{2} - 1$ and an isolated complex root of the form

$$\lambda_{BR}^{\sigma_1} = -i(\zeta/2 + \xi^{\sigma_1} + \epsilon^{\sigma_1}) = -i(\zeta/2 - \xi^{\sigma_1} + \epsilon^{\sigma_1}) + \delta^{\sigma_1} \frac{\pi}{2}, \quad (40)$$

where $\epsilon^{\sigma_1}$ is a finite length correction which becomes exponentially small for large L, i.e. $\epsilon^{\sigma_1} = O(L^{-\infty})$.

$$(B) \quad h^{\sigma_2} < h_{cr}^{(1)} < h^{\sigma_1} < h_{cr}^{(2)} (\{\sigma_1, \sigma_2\} = \{+,-\}).$$

The ground state is in the sector $N = \frac{L}{2}$ (with magnetization 0). It is given by $\frac{L}{2}$ real roots with adjacent quantum numbers $n_j = 1, \dots, \frac{L}{2}$.

$$(C) \quad h^{\sigma_2} < h_{cr}^{(1)} < h_{cr}^{(2)} < h^{\sigma_1} (\{\sigma_1, \sigma_2\} = \{+,-\}).$$

The ground state is in the sector $N = \frac{L}{2}$ (with magnetization 0). It is given by $\frac{L}{2} - 1$ real roots with adjacent quantum numbers $n_j = 1, \dots, \frac{L}{2} - 1$ and an isolated complex root of the form (40).

### For L odd:

$$(A') \quad h^{\pm} < h_{cr}^{(1)} \text{ with } h^{+} + h^{-} < 0.$$

The ground state is in the sector $N = \frac{L-1}{2}$ (with magnetization +1/2). It is given by $N = \frac{L-1}{2}$ real Bethe roots with adjacent quantum numbers.

$$(B') \quad h^{\pm} > -h_{cr}^{(1)} \text{ with } h^{+} + h^{-} > 0.$$

The ground state is in the sector $N = \frac{L+1}{2}$ (with magnetization -1/2). Computations of physical quantities in this case can simply be obtained from the previous case (A') by using the invariance of the model under the reversal of all spins together with a change of sign of the boundary fields $h^{\pm}$.

Note that there are no holes in the above configurations of Bethe roots (A) to (B').

Computations of physical quantities in the case (C) can also be obtained from cases (A) and (B) by using the invariance of the model under the reversal of all spins together with a change of sign of the boundary fields $h^{\pm}$.

As in [2], the isolated complex root of the form (40) will be called *boundary root* (it was called *boundary 1-string* in [15]). It was shown in [2] that such a boundary root plays an important role in the computation of physical quantities, such as the boundary magnetization at the thermodynamic limit.
## 2.4 Changing the boundary field $h_-$

In the following, we shall consider the effects of a change of the boundary field $h_-$ in the Hamiltonian (1). Hence, we consider a local quench in which we pass from a Hamiltonian $H_1$ with left boundary field $h_1^-$ to a Hamiltonian $H_2$ with left boundary field $h_2^-$, the other parameters $\Delta$ and $h^+$ of the Hamiltonian remaining unchanged.

It is important to remark that such a quench leaves the boundary monodromy matrix (4) unchanged, since the latter does not depend on $h^-$, see (10). Hence, the Bethe states take the same algebraic form (13) for both Hamiltonians. It is therefore possible to compute their overlap by using the determinant representation (23). This is the purpose of the next sections. In particular, our aim is to compute the normalized overlap between the ground state $|\{\lambda\}\rangle$, with Bethe roots that we shall denote by $\{\lambda\}$, of the Hamiltonian $H_1$, and the ground state $|\{\mu\}\rangle$, with Bethe roots that we shall denote by $\{\mu\}$, of the Hamiltonian $H_2$.

In the following, we shall adopt a subscript 1, respectively a subscript 2, for all quantities which depend on the boundary field $h_1^-$ (parametrized by $\xi_1^-$) of the Hamiltonian $H_1$, respectively on the boundary field $h_2^-$ (parametrized by $\xi_2^-$) of the Hamiltonian $H_2$. In particular, the Bethe equations for the ground state of $H_1$ are

$$a_1(\lambda_j|\{\lambda\}) = 1, \quad j = 1, \dots, N_1, \qquad (41)$$

whereas the Bethe equations for the ground state of $H_2$ are

$$a_2(\mu_j|\{\mu\}) = 1, \quad j = 1, \dots, N_2. \qquad (42)$$

The respective transfer matrix eigenvalues will be denoted by $\tau_1(v|\{\lambda\})$ and $\tau_2(v|\{\mu\})$.

# 3 Cauchy determinant representation for the overlap

The main goal of this section is to provide an exact representation for the overlap between the ground states of the two Hamiltonians $H_1$ and $H_2$, in the form of a generalised Cauchy determinant. More precisely, we consider the following normalized version of the overlap:

$$S(\{\lambda\}, \{\mu\}) = \frac{\langle\{\lambda\}|\{\mu\}\rangle \langle\{\mu\}|\{\lambda\}\rangle}{\langle\{\lambda\}|\{\lambda\}\rangle \langle\{\mu\}|\{\mu\}\rangle}. \qquad (43)$$

We recall that $\{\lambda\} = \{\lambda_1, \dots, \lambda_{N_1}\}$, respectively $\{\mu\} = \{\mu_1, \dots, \mu_{N_2}\}$, denotes the set of Bethe roots corresponding to the ground state of the Hamiltonian $H_1$ with first site boundary field $h_1^-$, respectively of the Hamiltonian $H_2$ with first site boundary field $h_2^-$. Note that the overlap (43) vanishes identically if $N_1 \neq N_2$, so that we restrict ourselves to cases for which the two ground states have the same number of Bethe roots $N_1 = N_2 = N$.

## 3.1 Determinant representation

We recall that, since the Hamiltonians $H_1$ and $H_2$ differ only by the value of the boundary field $h^-$ at first site (the boundary field $h^+$ remaining the same for both Hamiltonians), they share the same boundary monodromy matrix (4). Hence, the scalar products appearing in (43) can be represented by means of the formulas of section 2.2. More precisely, representing the first ratio of (43) by means of (23) and (25), and the second ratio by means of similar formulas in which we interchange the role of $\{\lambda\}$ and $\{\mu\}$, we obtain, after some slight simplifications:
$$ S(\{\lambda\}, \{\mu\}) = \prod_{k=1}^{N} \frac{a_1(-\mu_k)}{a_2(-\mu_k)} \frac{a_2(-\lambda_k)}{a_1(-\lambda_k)} \prod_{k=1}^{N} \prod_{l=1}^{N} \frac{s(\mu_k - i\zeta, \lambda_l)s(\lambda_k - i\zeta, \mu_l)}{s(\lambda_k - i\zeta, \lambda_l)s(\mu_k - i\zeta, \mu_l)} \\ \times \frac{\det_N \hat{H}_1(\lambda, \mu)}{\det_N \mathcal{M}_1(\lambda, \lambda)} \frac{\det_N \hat{H}_2(\mu, \lambda)}{\det_N \mathcal{M}_2(\mu, \mu)}. \quad (44) $$

Here we have defined, for two N-tuples $v = (v_1, ..., v_N)$ and $\omega = (\omega_1, ..., \omega_N)$ of complex numbers, the matrix $\hat{H}(v, \omega)$ with elements

$$ [\hat{H}(v, \omega)]_{jk} = a(\omega_k | \{v\}) [t(-\omega_k + v_j) - t(-\omega_k - v_j)] + t(\omega_k - v_j) - t(\omega_k + v_j), \quad (45) $$

in terms of the function (28). As mentioned in section (2.4), the index 1, respectively 2, means that the corresponding quantity depends explicitly on the boundary field $h_1^-$, respectively $h_2^-$. In other words,

$$ [\hat{H}_{\ell}(v, \omega)]_{jk} = a_{\ell}(\omega_k | \{v\}) [t(-\omega_k + v_j) - t(-\omega_k - v_j)] + t(\omega_k - v_j) - t(\omega_k + v_j), \quad (46) $$

for $\ell = 1, 2$, with $a_\ell$ being the exponential counting function (21) expressed in terms of the function $a_\ell$ (17) with boundary field $h_\ell^-$ parametrized by $\xi_\ell^-$. Similar notations are used for the Gaudin matrix (26), i.e.

$$ [M_{\ell}(v, v)]_{jk} = i \delta_{jk} a'_{\ell}(v_j | \{v\}) - 2\pi [K(v_j - v_k) - K(v_j + v_k)], \quad \ell = 1, 2. \quad (47) $$

Note that, to obtain the expression (44) from (23)-(24), we have used the following simple identity:

$$ \frac{\sinh \zeta}{s(\lambda, \mu)s(\lambda - i\zeta, \mu)} = \frac{t(\lambda - \mu) - t(\lambda + \mu)}{\sin(2\lambda - i\zeta) \sin(2\mu)}. \quad (48) $$

## 3.2 Gaudin extraction

In this subsection, we explain how to transform the ratios of determinants appearing in (44), adapting to the open case a procedure introduced in [14]. This will enable us to rewrite the overlap in terms of some generalised Cauchy determinant.

Let us consider the following ratio of determinants:

$$ \frac{\det_N[\hat{H}(v, \omega)]}{\det_N[\mathcal{M}(v, v)]} = \det_N[\mathcal{M}(v, v)^{-1} \hat{H}(v, \omega)], \quad (49) $$

involving the matrices (45) and (26). Here $v$ stands for a N-tuple of on-shell Bethe roots, whereas $\omega$ stands, for the moment, for any arbitrary N-tuple of complex numbers. Our aim is to compute the elements of the matrix $\mathcal{F}(v, \omega) = \mathcal{M}(v, v)^{-1} \hat{H}(v, \omega)$. The latter are given as the unique solution to the following system of linear equations:

$$ i a'(v_j | \{v\}) [\mathcal{F}(v, \omega)]_{jk} - 2\pi \sum_l [K(v_j - v_l) - K(v_j + v_l)] [\mathcal{F}(v, \omega)]_{lk} \\ = a(\omega_k | \{v\}) [t(-\omega_k + v_j) - t(-\omega_k - v_j)] + t(\omega_k - v_j) - t(\omega_k + v_j). \quad (50) $$

Following the strategy used in [14], let us look for a function $G(u, w) \equiv G_{\{v\}}(u, w)$ of two complex variables $u, w \in D_\zeta = \{z \in \mathbb{C} \mid |\Im(z)| < \zeta\}$, solution of the equation

$$ G(u, w) + 2\pi i \sum_{\ell} [K(u - v_{\ell}) - K(u + v_{\ell})] \frac{G(v_{\ell}, w)}{a'(v_{\ell} | \{v\})} \\ = a(w | \{v\}) [t(u - w) - t(-u - w)] + t(-u + w) - t(u + w). \quad (51) $$
<figure><img src="image_5.png" /><figcaption>Figure 1: Contour $\Gamma_{\epsilon}$.</figcaption></figure>

Note that, if we manage to construct such a solution to (51), it will give the elements of $\mathcal{F}(\nu, \omega)$ by setting

$$i\alpha'(\nu_j|\{\nu\})[\mathcal{F}(\nu, \omega)]_{jk} = G(\nu_j, \omega_k), \quad (52)$$

We shall moreover look for a solution of (51) such that

(i) $u \mapsto G(u, w)$ is a meromorphic function in the strip $D_{\zeta}$, with no poles at the points $\nu_{\ell}, 1 \le \ell \le N$,

(ii) $u \mapsto G(u, w)$ is odd: $G(-u, w) = -G(u, w)$.

Under the assumptions (i) and (ii), (51) can be rewritten as

$$G(u, w) + \oint_{\Gamma_{\nu, \mathbb{R}}^{\pm}} K(u-z) \frac{G(z, w)}{a(z|\{\nu\})-1} dz + 2\pi i \sum_{\ell \in \mathbb{Z}_{\nu}} [K(u-\nu_{\ell}) - K(u+\nu_{\ell})] \frac{G(\nu_{\ell}, w)}{a'(\nu_{\ell}|\{\nu\})} \\ = a(w|\{\nu\})[t(u-w) - t(-u-w)] + t(-u+w) - t(u+w), \quad (53)$$

in which the contour $\Gamma_{\nu, \mathbb{R}}^{\pm}$ surrounds, with index 1, the elements of the set $(\{\nu\} \cup \{-\nu\}) \cap \mathbb{R}$ (and no other poles of the integrand), whereas the remaining sum runs over the set $\mathbb{Z}_{\nu}$ of indices corresponding to the subset of complex Bethe roots: $\{\nu_{\ell}\}_{\ell \in \mathbb{Z}_{\nu}} = \{\nu\} \cap (\mathbb{C} \setminus \mathbb{R})$.

It is clear that any function $G(u, w)$ solution of (53) should be a $\pi$-periodic function in $u$: $G(u + \pi, w) = G(u, w)$. Moreover, taking into account the requirement that $u \mapsto G(u, w)$ is a meromorphic function, (53) implies that the only poles of this function that may be on the real axis are those at $u = \pm w \mod \pi$ (if $w \in \mathbb{R}$). The corresponding residues can be easily derived from the R.H.S. of (53):

$$\mathrm{Res}_{u=\pm w} [G(u, w)] = i(a(w|\{\nu\}) - 1). \quad (54)$$

Hence, for $\epsilon > 0$ and small enough, we can deform the contour $\Gamma_{\nu, \mathbb{R}}^{\pm}$ into a rectangle $\Gamma_{\epsilon}$ with vertices at $(-\pi/2, i\epsilon)$, $(-\pi/2, -i\epsilon)$, $(\pi/2, i\epsilon)$, $(\pi/2, -i\epsilon)$. This can be done provided we subtract the contribution of the potential other poles that become encircled by the contour. If there are no holes, i.e. if the function $a(u|\{\nu\}) - 1$ has no other roots on the interval $(0, \pi/2)$ than the real Bethe roots $\nu_j$, the contour integral can be rewritten as

$$\oint_{\Gamma_{\nu, \mathbb{R}}^{\pm}} K(u-z) \frac{G(z, w)}{a(z|\{\nu\})-1} dz = \oint_{\Gamma_{\epsilon}} K(u-z) \frac{G(z, w)}{a(z|\{\nu\})-1} dz \\ - 2\pi i \delta_{w \in \mathbb{R}} [iK(u-w) - iK(u+w)a(w|\{\nu\})], \quad (55)$$
in which the notation $\delta_{w \in \mathbb{R}}$ means that the last term is non-zero only if $w \in \mathbb{R}$. Hence (53) can be rewritten as

$$ G(u, w) + \oint_{\Gamma_e} K(u-z) \frac{G(z, w)}{a(z|\{v\})-1} dz + 2\pi i \sum_{l \in \mathbb{Z}_v} [K(u-v_l) - K(u+v_l)] \frac{G(v_l, w)}{a'(v_l|\{v\})} \\ = \begin{cases} [t(u-w) + t(u+w)][a(w|\{v\})-1], & \text{if } w \in \mathbb{R}, \\ a(w|\{v\})[t(u-w) - t(-u-w)] + t(-u+w) - t(u+w), & \text{if } w \notin \mathbb{R}. \end{cases} \quad (56) $$

Let us first suppose that all the Bethe roots $v_1, \dots, v_N$ are real. It is not difficult to see, following standards arguments [2,24], that the quantity $a(w|\{v\})$ has the following behaviour in $L$ for $w$ above and below the real axis:

$$ a(w|\{v\}) = O(L^{+\infty}), \quad \text{for} \quad -\zeta < \Im(w) < 0, \quad (57) $$

$$ a(w|\{v\}) = O(L^{-\infty}), \quad \text{for} \quad 0 < \Im(w) < \zeta. \quad (58) $$

Here and in the following, the notation $O(L^{+\infty})$ (respectively $O(L^{-\infty})$) means that the corresponding quantity is diverging exponentially with $L$ (respectively is vanishing exponentially with $L$). Hence, the equation (56) simplifies drastically for large $L$:

$$ G(u, w) + \int_{-\frac{\pi}{2}+i0}^{\frac{\pi}{2}+i0} K(u-z) G(z, w) dz \\ = \begin{cases} [t(u-w) + t(u+w)][a(w|\{v\})-1] + O(L^{-\infty}), & \text{if } w \in \mathbb{R}, \\ t(-u+w) - t(u+w) + O(L^{-\infty}), & \text{if } 0 < \Im(w) < \zeta, \\ a(w|\{v\})[t(u-w) - t(-u-w) + O(L^{-\infty})], & \text{if } -\zeta < \Im(w) < 0. \end{cases} \quad (59) $$

This integral equation has to be compared with the Lieb equation (34). Note that all the functions involved in (34) are holomorphic functions in the strip $D_{\zeta/2} = \{z \in \mathbb{C} \mid |\Im(z)| < \zeta/2\}$, so that (34) remains valid in the whole strip $\lambda \in D_{\zeta/2}$. Moreover, the integration contour in (34) can be shifted by a fixed imaginary value as soon as we do not cross poles. Hence, up to exponentially small corrections in $L$, the integral equation (59) admits a solution which is given by a linear combination of solutions of (34):

$$ G(u, w) = i[a(w|\{v\}) - 1][\bar{\rho}(u, w) + O(L^{-\infty})], \quad (60) $$

in which $\bar{\rho}$ is given in term of the function $\rho$ (35) as

$$ \bar{\rho}(u, w) = -i\pi[\rho(u-w-i\zeta/2) + \rho(u+w-i\zeta/2)] \\ = \frac{\vartheta'_1(0)}{\vartheta_2(0)} \left[ \frac{\vartheta_2(u-w)}{\vartheta_1(u-w)} + \frac{\vartheta_2(u+w)}{\vartheta_1(u+w)} \right]. \quad (61) $$

Note that (60)-(61) provides a solution for all cases considered in (59) thanks to the parity and quasi-periodicity properties of the function $\rho$ (35). Note also that the function (61) satisfies the properties (i) and (ii), so that it also provides a solution of (51). Hence the elements of the matrix $\mathcal{F}(v, \omega) = \mathcal{M}(v, v)^{-1} \hat{H}(v, \omega)$ are computed by (52) up to exponentially small corrections in $L$, and

$$ \frac{\det_N[\hat{H}(v, \omega)]}{\det_N[\mathcal{M}(v, v)]} = \det_N[\mathcal{F}(v, \omega)] = \prod_{j=1}^{N} \frac{a(\omega_j|\{v\})-1}{a'(v_j|\{v\})} \det_N[\bar{\rho}(v_j, \omega_k)] (1+O(L^{-\infty})). \quad (62) $$
Let us now suppose that the set of Bethe roots {$v$} contains only real roots except for one single boundary complex root $v_{BR}^{\sigma} = -i(\frac{\zeta}{2} + \xi^{\sigma} + \epsilon^{\sigma})$, $\sigma = \pm$, with $\epsilon^{\sigma} = O(L^{-\infty})$. Then, the large $L$ behaviour (57)-(58) is still valid provided that $w$ remains at finite distance from the zeroes and the poles of the function $\alpha$, which we suppose from now on: this is notably the case when $w$ is in the vicinity of the real axis, or when it coincides with a boundary root built on a different boundary parameter as those appearing in $\alpha$. The difference with the previous case is that we have a priori to add to the left hand side of (59) an additional term of the form

$$[K(u - v_{BR}^{\sigma}) - K(u + v_{BR}^{\sigma})] \frac{G(v_{BR}^{\sigma}, w)}{a'(v_{BR}^{\sigma}|\{v\})}. \qquad (63)$$

Let us however remark that

$$a'(v_{BR}^{\sigma}|\{v\}) = 2Li \xi'(v_{BR}^{\sigma}|\{v\}) \underset{L \to \infty}{\sim} i g'(v_{BR}^{\sigma}) = O(\frac{1}{\epsilon^{\sigma}}), \qquad (64)$$

which diverges exponentially with $L$. This means that the correction induced by the presence of the additional term (63) to the solution of (59) is in fact exponentially small in $L$ with respect to the leading order, i.e. that the solution (60)-(61)-(62) still holds even in the presence of a boundary root.

This result can be applied directly to the computation of the two ratios in (44):

$$\frac{\det_N \hat{H}_1(\lambda, \mu)}{\det_N M_1(\lambda, \lambda)} = \prod_{i=1}^{N} \frac{a_1(\mu_i|\{\lambda\}) - 1}{a'_1(\lambda_i|\{\lambda\})} \det_N [\bar{\rho}(\lambda_j, \mu_k)] (1 + O(L^{-\infty})), \qquad (65)$$

$$\frac{\det_N \hat{H}_2(\mu, \lambda)}{\det_N M_2(\mu, \mu)} = \prod_{i=1}^{N} \frac{a_2(\lambda_i|\{\mu\}) - 1}{a'_2(\mu_i|\{\mu\})} \det_N [\bar{\rho}(\mu_j, \lambda_k)] (1 + O(L^{-\infty})). \qquad (66)$$

### 3.3 Product formula for the overlap

The determinant in (62) (or equivalently in (65)-(66)) is a generalised Cauchy determinant. It can easily be computed by means of usual arguments for elliptic functions:

$$\det_N [\bar{\rho}(v_j, \omega_k)] = \left(\frac{\vartheta'_1(0)}{\vartheta_2(0)}\right)^N \det_N \left[\frac{\vartheta_2(v_j - \omega_k)}{\vartheta_1(v_j - \omega_k)} + \frac{\vartheta_2(v_j + \omega_k)}{\vartheta_1(v_j + \omega_k)}\right] \\
= \left(\frac{\vartheta'_1(0)}{\vartheta_2(0)}\right)^N \prod_{i=1}^N [\vartheta_1(2v_i, q^2)\vartheta_4(2\omega_i, q^2)] \\
\times \frac{\prod_{k<j} \vartheta_1(v_j + v_k)\vartheta_1(v_j - v_k)\vartheta_1(\omega_k + \omega_j)\vartheta_1(\omega_k - \omega_j)}{\prod_{k=1}^N \prod_{j=1}^N \vartheta_1(v_j + \omega_k)\vartheta_1(v_j - \omega_k)}. \qquad (67)$$

This enables us to express the overlap (44) as a simple product. Introducing

$$\chi(u) = \frac{\tau_2(u|\{\mu\})}{\tau_1(u|\{\lambda\})} \\
= \frac{[a_2(u|\{\mu\}) - 1] \sin(u - i\xi_2 - i\zeta/2)}{[a_1(u|\{\lambda\}) - 1] \sin(u - i\xi_1 - i\zeta/2)} \prod_{l=1}^{N} \frac{s(u - i\zeta, \mu_l) s(u, \lambda_l)}{s(u - i\zeta, \lambda_l) s(u, \mu_l)}, \qquad (68)$$

to be the ratio of the two transfer matrix eigenvalues, and defining the following regularised theta function

$$\varphi(\lambda) = \frac{\vartheta_1(\lambda)}{\sin \lambda} = 2q^{1/4} \prod_{n=1}^{\infty} (1 - q^{2n} e^{2i\lambda})(1 - q^{2n} e^{-2i\lambda})(1 - q^{2n}), \qquad (69)$$
we finally obtain

$$S(\{\lambda\}, \{\mu\}) = \prod_{i=1}^{N} \frac{\chi(\lambda_i)}{\chi(\mu_i)} \prod_{i=1}^{N} \prod_{j=1}^{N} \frac{\varphi(\lambda_i + \lambda_j)\varphi(\lambda_i - \lambda_j)}{\varphi(\lambda_i + \mu_j)\varphi(\lambda_i - \mu_j)} \frac{\varphi(\mu_i + \mu_j)\varphi(\mu_i - \mu_j)}{\varphi(\mu_i + \lambda_j)\varphi(\mu_i - \lambda_j)} + O(L^{-\infty}). \quad (70)$$

This simple product formula is valid for large but finite L and N (we recall that for the ground state N is of order $\frac{L}{2}$), up to exponentially small corrections in L. To derive it, we have essentially used here the fact that the Bethe states that we consider are characterised by real Bethe roots, with possibly some boundary root. We have also used the absence of holes in the distribution of real roots.

In the next section, we will compute the thermodynamic limit of this product, taking more specifically into account the properties of the ground states that we consider.

# 4 Taking the thermodynamic limit

In this section, we explain how to take the thermodynamic limit of the product formula (70) obtained for the overlap in the previous section. We have to distinguish several cases according to the configurations of the two sets of ground state Bethe roots {$\lambda$} and {$\mu$}, see section 2.3.

## 4.1 The case of real Bethe roots

Let us first present the computation of the overlap (70) in the simplest case in which all Bethe roots are real. According to the description of the ground state recalled in section 2.3, this situation occurs for the following respective configurations of the boundary magnetic fields:

(a) case (A') for both states for L odd: $h^+, h_1^-, h_2^- < h_{cr}^{(1)}$ and $h_1^-, h_2^- < -h_+$;

(b) case (B) for both states for L even, which means that we are in one of the two following situations:

$$ (b1) \quad h_1^-, h_2^- < h_{cr}^{(1)} < h^+ < h_{cr}^{(2)}, $$

$$ (b2) \quad h^+ < h_{cr}^{(1)} < h_1^-, h_2^- < h_{cr}^{(2)}. $$

### 4.1.1 General strategy

When considering the thermodynamic limit of (70), we have to deal with factors of the form

$$F_f(v|\{\lambda\}, \{\mu\}) = \prod_{j=1}^{N} \frac{f(v-\lambda_j)f(v+\lambda_j)}{f(v-\mu_j)f(v+\mu_j)}, \quad (71)$$

in which f is a $\pi$-periodic analytic function which does not vanish on the real axis. The thermodynamic limit of products of the form (71), in this simplest case without complex roots nor holes, can be computed as

$$ \begin{aligned} F_f(v|\{\lambda\}, \{\mu\}) &= \exp \left\{ -\sum_{j=1}^{N} \left[ \log f(v-\mu_j) + \log f(v+\mu_j) - \log f(v-\lambda_j) - \log f(v+\lambda_j) \right] \right\} \\ &= \exp \left\{ -\frac{L}{\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \log f(v-w) \left[ \hat{\xi}'_2(w|\{\lambda\}) - \hat{\xi}'_1(w|\{\mu\}) \right] dw + O(L^{-\infty}) \right\}, \quad (72) \end{aligned} $$
in which we have used the sum-to-integral transform, see Corollary B.1 of [2]. The difference of the two counting functions is itself given as

$$ \begin{aligned} \hat{\xi}_2(w|\{\mu\}) - \hat{\xi}_1(w|\{\lambda\}) &= \frac{g_2(w) - g_1(w)}{2L} \\ &\quad + \frac{1}{2L} \sum_{k=1}^{N} [\theta(w - \mu_k) + \theta(w + \mu_k) - \theta(w - \lambda_k) - \theta(w + \lambda_k)] \\ &= \frac{1}{L} \hat{\xi}_d(w) + O(L^{-\infty}), \end{aligned} \quad (73) $$

in which $\hat{\xi}_d$ satisfies the equation

$$ \hat{\xi}_d(\lambda) = \frac{1}{2\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \theta(\lambda - \mu) \hat{\xi}'_d(\mu) d\mu + \frac{g_2(\lambda) - g_1(\lambda)}{2}. \quad (74) $$

In the cases we consider here, $\hat{\xi}_1(\frac{\pi}{2}) = -\hat{\xi}_1(-\frac{\pi}{2}) = \hat{\xi}_2(\frac{\pi}{2}) = -\hat{\xi}_2(-\frac{\pi}{2}) = \frac{(N+1)\pi}{L}$, so that $\hat{\xi}_d(\frac{\pi}{2}) = -\hat{\xi}_d(-\frac{\pi}{2}) = 0$. Therefore (74) can alternatively be rewritten as

$$ \hat{\xi}_d(\lambda) + \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} K(\lambda - \mu) \hat{\xi}_d(\mu) d\mu = \frac{g_2(\lambda) - g_1(\lambda)}{2}, \quad (75) $$

and (72) as

$$ F_f(v|\{\lambda\}, \{\mu\}) = \exp \left[ -\frac{1}{\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{f'(\lambda - \mu)}{f(\lambda - \mu)} \hat{\xi}_d(\mu) d\mu \right] [1 + O(L^{-\infty})]. \quad (76) $$

The r.h.s. of (75) is a $2\pi$-periodic odd function:

$$ g_2(\lambda) - g_1(\lambda) = g_+(\lambda) - g_+(\lambda), \quad (77) $$

with

$$ g_+(\lambda) = -i\epsilon \log \frac{1 - e^{2i\lambda}(p_2q)^\epsilon}{1 - e^{2i\lambda}(p_1q)^\epsilon}, \quad \text{with} \quad p_j = e^{-2\xi_j^{-}}, \quad (78) $$

and $\epsilon = 1$ if $|p_jq| < 1$, i.e. if we are in case (a) or (b1), whereas $\epsilon = -1$ if $|p_jq| > 1$, i.e. if we are in case (b2). Hence, the integral equation (75) can be solved in terms of Fourier series, and the solution $\hat{\xi}_d(\lambda)$ of (75) can also be written in the form

$$ \hat{\xi}_d(\lambda) = \hat{\xi}_+(\lambda) - \hat{\xi}_+(\lambda), \quad (79) $$

where $\hat{\xi}_+(\lambda)$ is a power series of $e^{2i\lambda}$ (with only positive powers) which satisfies the integral equation

$$ \hat{\xi}_+(\lambda) + \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} K(\lambda - \mu) \hat{\xi}_+(\mu) d\mu = \frac{g_+(\lambda)}{2}. \quad (80) $$

Taking into account the explicit form of the kernel $K(\lambda)$ and its Fourier coefficients, one can rewrite this equation as

$$ \hat{\xi}_+(\lambda) + \hat{\xi}_+(\lambda + i\zeta) = \frac{g_+(\lambda)}{2}. \quad (81) $$

Introducing

$$ a_+(u) = \exp[2i\hat{\xi}_+(\lambda)], \quad u = e^{2i\lambda}, \quad (82) $$
we get the following functional equation equivalent to the integral equation (80):

$$a_+(u) a_+(uq^2) = \left( \frac{1 - u(p_2q)^\epsilon}{1 - u(p_1q)^\epsilon} \right)^\epsilon. \qquad (83)$$

The unique solution of this equation can be easily expressed

$$a_+(u) = \left( \frac{(u(p_2q)^\epsilon; q^4)_\infty}{(uq^2(p_2q)^\epsilon; q^4)_\infty} \frac{(uq^2(p_1q)^\epsilon; q^4)_\infty}{(u(p_1q)^\epsilon; q^4)_\infty} \right)^\epsilon, \qquad (84)$$

in terms of the the q-Pochhammer symbol

$$(x; \alpha)_\infty = \prod_{n=0}^{\infty} (1 - x\alpha^n). \qquad (85)$$

It remains to report the expression (84) into (76) and to compute the integral, so as to obtain an evaluation of the product (71) in the thermodynamic limit, up to exponentially small corrections in L.

In particular, if f is of the form

$$f_{p,\pm}(\lambda) = 1 - p e^{\pm 2i\lambda}, \quad \text{with } |p| < 1, \qquad (86)$$

the integral in (76) can easily be computed. Using the representation of $f'_{p,\pm}(\lambda)/f_{p,\pm}(\lambda)$ in Fourier series, we straightforwardly obtain that

$$\exp \left[ -\frac{1}{\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{f'_{p,\pm}(\lambda - \mu)}{f_{p,\pm}(\lambda - \mu)} \hat{\xi}_d(\mu) d\mu \right] = a_+(p u^{\pm 1}), \quad u = e^{2i\lambda}. \qquad (87)$$

### 4.1.2 Computation of the overlap

We can now apply the previous strategy so as to compute the overlap (70) in the case where all Bethe roots are real.

In particular, from the infinite product representation of the function $\varphi$ (69) in terms of functions of the form (86), we readily obtain

$$F_{\varphi}(v|\{\lambda\}, \{\mu\}) = \prod_{j=1}^{N} \frac{\varphi(v - \lambda_j) \varphi(v + \lambda_j)}{\varphi(v - \mu_j) \varphi(v + \mu_j)} \\ = \left[ \prod_{n=1}^{\infty} a_+(q^{2n}u) a_+(q^{2n}u^{-1}) \right] [1 + O(L^{-\infty})]. \qquad (88)$$

Let us now consider the ratio of transfer matrix eigenvalues $\chi(\lambda_j)$ evaluated at a Bethe root $\lambda_j$. Using the Bethe equations $a_1(\lambda_i|\{\lambda\}) = 1$ we can express the quantity $a_2(\lambda_j|\{\mu\})$ as

$$a_2(\lambda_j|\{\mu\}) = \frac{\sin(\lambda_j + i\xi_2^- + i\zeta/2) \sin(\lambda_j - i\xi_1^- - i\zeta/2) \phi(\lambda_j - i\zeta)}{\sin(\lambda_j + i\xi_1^- + i\zeta/2) \sin(\lambda_j - i\xi_2^- - i\zeta/2) \phi(\lambda_j + i\zeta)}, \qquad (89)$$

in which we have introduced the shortcut notation

$$\phi(v) = \prod_{k=1}^{N} \frac{\sin(v - \lambda_k) \sin(v + \lambda_k)}{\sin(v - \mu_k) \sin(v + \mu_k)}. \qquad (90)$$
Therefore (68) takes the following form

$$
\begin{align}
\chi(\lambda_j) &= \frac{\phi'(\lambda_j)}{a'_1(\lambda_j|\{\lambda\})} \left[ \frac{\sin(\lambda_j + i\xi_2^- + i\zeta/2)}{\sin(\lambda_j + i\xi_1^- + i\zeta/2)} \frac{1}{\phi(\lambda_j + i\zeta)} - \frac{\sin(\lambda_j - i\xi_2^- - i\zeta/2)}{\sin(\lambda_j - i\xi_1^- - i\zeta/2)} \frac{1}{\phi(\lambda_j - i\zeta)} \right] \nonumber \\
&= \frac{\phi'(\lambda_j)}{a'_1(\lambda_j|\{\lambda\})} \left[ \frac{(1-u_j(p_2q)^{-1})(1-u_j^{-1}(p_2q)^{-1})}{(1-u_j(p_1q)^{-1})(1-u_j^{-1}(p_1q)^{-1})} \right]^{\frac{1-\epsilon}{2}} \left(\frac{p_1}{p_2}\right)^{\frac{\epsilon}{2}} \nonumber \\
&\quad \times \left[ \left(\frac{1-u_j(p_2q)^{\epsilon}}{1-u_j(p_1q)^{\epsilon}}\right)^{\epsilon} \frac{1}{a_+(q^2 u_j)} - \left(\frac{1-u_j^{-1}(p_2q)^{\epsilon}}{1-u_j^{-1}(p_1q)^{\epsilon}}\right)^{\epsilon} \frac{1}{a_+(q^2 u_j^{-1})} \right] [1+O(L^{-\infty})], \tag{91}
\end{align}
$$

in which we have set $u_j = e^{2i\lambda_j}$ and used that

$$
\phi(\lambda \pm i\zeta) = F_{f_{q2,+}}(\lambda|\{\lambda\}, \{\mu\}) = a_+(q^2 e^{\pm 2i\lambda}) [1 + O(L^{-\infty})]. \quad (92)
$$

This expression can be simplified by means of the functional equation (83), and we obtain

$$
\chi(\lambda_j) = \frac{\phi'(\lambda_j)}{2iL \xi_1'(\lambda_j|\{\lambda\})} \left[ \frac{(1-u_j(p_2q)^{-1})(1-u_j^{-1}(p_2q)^{-1})}{(1-u_j(p_1q)^{-1})(1-u_j^{-1}(p_1q)^{-1})} \right]^{\frac{1-\epsilon}{2}} \left(\frac{p_1}{p_2}\right)^{\frac{\epsilon}{2}} \\
\times \left[ a_+(u_j) - a_+(u_j^{-1}) \right] [1 + O(L^{-\infty})]. \quad (93)
$$

Let us now consider the function

$$
\psi(v) = 2iL \frac{\xi_1'(v|\{\lambda\})}{\phi'(v)}. \qquad (94)
$$

We can evaluate it by means of similar arguments to those used in [25]. Let us first remark
that

$$
2\pi \sum_{a=1}^{N} [K(\mu - \lambda_a) - K(\mu + \lambda_a)] \frac{1}{\phi'(\lambda_a)} = i \left[ \frac{1}{\phi(\mu + i\zeta)} - \frac{1}{\phi(\mu - i\zeta)} \right], \quad (95)
$$

which can easily be proved by comparing poles and residues of the r.h.s and l.h.s as well as their
behaviour at infinity. Replacing as usual the sum over Bethe roots by an integral by means of
Proposition B.1 of [2] and using (92), we obtain the following equation for the function ψ(λ):

$$
\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} K(\mu - \nu) \psi(\nu) d\nu = - \left[ \frac{1}{a_+(e^{2i\mu q^2})} - \frac{1}{a_+(e^{-2i\mu q^2})} \right] [1 + O(L^{-\infty})]. \quad (96)
$$

Using the explicit form of the kernel $K(\lambda)$ and analyticity of the function $1/a_+(u)$ inside the
unit circle we can easily solve this equation:

$$
\psi(\mu) = - \left[ \frac{1}{a_+(e^{2i\mu})} - \frac{1}{a_+(e^{-2i\mu})} \right] [1 + O(L^{-\infty})]. \quad (97)
$$

Inserting this result into (93) we obtain the following very simple formula for the ratio of
transfer matrix eigenvalues evaluated in $\lambda_j$:

$$
\chi(\lambda_j) = \left(\frac{p_1}{p_2}\right)^{\frac{\epsilon}{2}} \left[a_+(q^{1-\epsilon}u_j) a_+(q^{1-\epsilon}u_j^{-1})\right]^{\epsilon} [1 + O(L^{-\infty})]. \quad (98)
$$
Combining this result with (88) we obtain

$$ \begin{aligned} \chi(\lambda_j) F_{\varphi}(\lambda_j | \{\lambda\}, \{\mu\}) &= \left(\frac{p_1}{p_2}\right)^{\frac{\epsilon}{2}} \left[ \prod_{n=1}^{\infty} a_+(q^{2(n-\epsilon)} u_j) a_+(q^{2(n-\epsilon)} u_j^{-1}) \right] [1 + O(L^{-\infty})] \\ &= \left(\frac{p_1}{p_2}\right)^{\frac{\epsilon}{2}} \left[ \frac{(u_j q^{2-\epsilon} p_2^{\epsilon}; q^4)_{\infty}}{(u_j q^{2-\epsilon} p_1^{\epsilon}; q^4)_{\infty}} \frac{(u_j^{-1} q^{2-\epsilon} p_2^{\epsilon}; q^4)_{\infty}}{(u_j^{-1} q^{2-\epsilon} p_1^{\epsilon}; q^4)_{\infty}} \right]^{\epsilon} [1 + O(L^{-\infty})], \quad (99) \end{aligned} $$

in which we have used (84).

The ratio $\chi(\mu_j)$ of transfer matrix eigenvalues evaluated in $\mu_j$ can be expressed similarly as in (98), by simply replacing $u_j = e^{2i\lambda_j}$ by $v_j = e^{2i\mu_j}$ in the expression. Hence, the overlap is given as a product over $j$ of ratios of the form (99) in terms of $\lambda_j$ and $\mu_j$, that we can again evaluate following what has been done in section 4.1.1, using infinite product representation of the q-Pochhammer symbol in terms of functions of the form (86):

$$ \begin{aligned} S(\{\lambda\}, \{\mu\}) &= \left[ \prod_{n=0}^{\infty} \frac{a_+(p_2^{\epsilon} q^{2-\epsilon+4n})}{a_+(p_1^{\epsilon} q^{2-\epsilon+4n})} \right]^{\epsilon} [1 + O(L^{-\infty})] \\ &= \frac{(p_1^{2\epsilon} q^2; q^4, q^4)_{\infty} (p_2^{2\epsilon} q^2; q^4, q^4)_{\infty} ((p_1 p_2)^{\epsilon} q^4; q^4, q^4)_{\infty}^2}{(p_1^{2\epsilon} q^4; q^4, q^4)_{\infty} (p_2^{2\epsilon} q^4; q^4, q^4)_{\infty} ((p_1 p_2)^{\epsilon} q^2; q^4, q^4)_{\infty}^2} + O(L^{-\infty}). \end{aligned} \quad (100) $$

We recall that $p_i = e^{-2\xi_i^-}$, $i = 1, 2$, and that $\epsilon = 1$ if we are in case (a) or (b1), whereas $\epsilon = -1$ if we are in case (b2). Here $(x; q_1, q_2)_{\infty}$ denotes the double q-Pochhammer symbol defined as

$$ (x; q_1, q_2)_{\infty} = \prod_{n_1, n_2=0}^{\infty} (1 - x q_1^{n_1} q_2^{n_2}). \quad (101) $$

## 4.2 Cases with presence of a boundary root

We now explain how the previous analytical computations are modified when one of the two states, or both, are described by a set of Bethe roots with a boundary root of the form (40).

### 4.2.1 Presence of two boundary roots $\lambda_{BR}^{\sigma_1}$ and $\mu_{BR}^{\sigma_2}$

When considering a product of the form (71), we have now to take into account the presence of these boundary roots $\lambda_{BR}^{\sigma_1} = -i(\zeta_2 + \xi_1^{\sigma_1} + \epsilon_1^{\sigma_1})$ and $\mu_{BR}^{\sigma_2} = -i(\zeta_2 + \xi_2^{\sigma_2} + \epsilon_2^{\sigma_2})$, with $\sigma_1, \sigma_2 \in \{+, -\}$ and $\epsilon_1^{\sigma_1}, \epsilon_2^{\sigma_2} = O(L^{-\infty})$:

$$ F_f(v|\{\lambda\}, \{\mu\}) = F_{f,\text{BR}}(v) \exp \left[ -\frac{1}{\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{f'(v-\mu)}{f(v-\mu)} \hat{\xi}_d(\mu) d\mu \right] [1 + O(L^{-\infty})], \quad (102) $$

where $\hat{\xi}_d$ now satisfies the equation<fnref n="2" />

$$ \begin{aligned} \hat{\xi}_d(\lambda) + \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} K(\lambda - \mu) \hat{\xi}_d(\mu) d\mu &= \frac{g_2(\lambda) - g_1(\lambda)}{2} \\ &\quad + \frac{1}{2} \sum_{\sigma=\pm} \left[ \theta(\lambda - i\sigma(\frac{\zeta}{2} + \xi_2^{\sigma_2})) - \theta(\lambda - i\sigma(\frac{\zeta}{2} + \xi_1^{\sigma_1})) \right], \end{aligned} \quad (103) $$

and

$$ F_{f,\text{BR}}(v) = \frac{f(v + i\frac{\zeta}{2} + i\xi_1^{\sigma_1}) f(v - i\frac{\zeta}{2} - i\xi_1^{\sigma_1})}{f(v + i\frac{\zeta}{2} + i\xi_2^{\sigma_2}) f(v - i\frac{\zeta}{2} - i\xi_2^{\sigma_2})}. \quad (104) $$

<footnote n="2">²We have here $\hat{\xi}_1(\frac{\pi}{2}) = -\hat{\xi}_1(-\frac{\pi}{2}) = \hat{\xi}_2(\frac{\pi}{2}) = -\hat{\xi}_2(-\frac{\pi}{2}) = \frac{N\pi}{L}$, so that $\hat{\xi}_d(\frac{\pi}{2}) = -\hat{\xi}_d(-\frac{\pi}{2}) = 0$.</footnote>
**Case with $\lambda_{BR}^+$ and $\mu_{BR}^+$:** Let us first consider the case in which the two boundary roots are $\lambda_{BR}^+$ and $\mu_{BR}^+$. This occurs when

$$h_1^-, h_2^- < h^+, \text{ with } |h^+| < h_{cr}^{(1)}, \text{ or } h_1^-, h_2^- < h_{cr}^{(1)} < h_{cr}^{(2)} < h^+. \quad (105)$$

This case is particularly simple since these two boundary roots coincide up to exponentially small corrections in $L$. In other words, the extra factor $F_{f,\text{BR}}(v)$ (104) in (102) is identically equal to 1, whereas the term in the second line of (103) vanishes identically. Hence, the presence of the boundary roots does not induce any changes in the computations with respect to the case with only real roots, and the overlap is given by the expression (100) with $\epsilon = 1$.

**Case with $\lambda_{BR}^-$ and $\mu_{BR}^-$:** Let us now consider the case in which the two boundary roots are $\lambda_{BR}^-$ and $\mu_{BR}^-$. This occurs when

$$h^+ < h_1^-, h_2^- \text{, with } |h_1^-|, |h_2^-| < h_{cr}^{(1)} \text{, or } h^+ < h_{cr}^{(1)} < h_{cr}^{(2)} < h_1^- \text{, } h_2^- \quad (106)$$

We can apply in this case similar arguments as in section 4.1.1 and decompose $\hat{\xi}_d$ as in (79), which leads to the following integral equation for $\hat{\xi}_+$:

$$\hat{\xi}_+(\lambda) + \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} K(\lambda - \mu) \hat{\xi}_+(\mu) d\mu = -\frac{i}{2} \log \left( \frac{1 - e^{2i\lambda} p_2 q}{1 - e^{2i\lambda} p_1 q} \frac{1 - e^{2i\lambda} p_2^{-1} q}{1 - e^{2i\lambda} p_1^{-1} q} \frac{1 - e^{2i\lambda} p_2 q^3}{1 - e^{2i\lambda} p_1 q^3} \right). \quad (107)$$

Note that, in the regimes (106) considered here, we have both $|p_j q| < 1$ and $|p_j^{-1} q| < 1$, $j = 1, 2$. This equation is equivalent to

$$a_+(u) a_+(u q^2) = \frac{1 - up_2 q}{1 - up_1 q} \frac{1 - up_2^{-1} q}{1 - up_1^{-1} q} \frac{1 - up_2 q^3}{1 - up_1 q^3}, \quad (108)$$

in which we have defined $a_+$ in terms of $\hat{\xi}_+$ as in (84). The unique solution of this equation is given by

$$a_+(u) = \frac{(up_2^{-1}q; q^4)_\infty (up_1^{-1}q^3; q^4)_\infty}{(up_1^{-1}q; q^4)_\infty (up_2^{-1}q^3; q^4)_\infty} \frac{1 - up_2q}{1 - up_1q}. \quad (109)$$

For $f$ of the form $f_{p,\pm}$ (86), the identity (87) still holds in terms of the function $a_+$ (109). Combining this expression with the extra factor (104) due to the presence of the boundary roots in (102), which in this case takes the form

$$F_{f_{p,\pm},\text{BR}}(\lambda) = \frac{(1 - pp_1 qu^{\pm 1})(1 - p(p_1 q)^{-1} u^{\pm 1})}{(1 - pp_2 qu^{\pm 1})(1 - p(p_2 q)^{-1} u^{\pm 1})}, \quad (110)$$

we obtain

$$F_{f_{p,\pm},\text{BR}}(\lambda) \exp \left[ -\frac{1}{\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{f'_{p,\pm}(\lambda - \mu)}{f_{p,\pm}(\lambda - \mu)} \hat{\xi}_d(\mu) d\mu \right] = \tilde{a}_+(pu^{\pm 1}), \quad (111)$$

with

$$\tilde{a}_+(u) = \frac{(up_2^{-1}q; q^4)_\infty (u(p_1q)^{-1}; q^4)_\infty}{(up_1^{-1}q; q^4)_\infty (u(p_2q)^{-1}; q^4)_\infty}. \quad (112)$$

Note that $\tilde{a}_+(u)$ coincides with the function obtained in (84) for $\epsilon = -1$, which therefore satisfies (83) with $\epsilon = -1$.

The rest of the computation is therefore nearly identical to what has been done in section 4.1.2 in the case $\epsilon = -1$. The functions $F_\varphi(v|\{\lambda\}, \{\mu\})$, $\phi(\lambda \pm i\zeta)$ and the ratio $\chi(\lambda_j)$ are then expressed respectively as in (88), (92) and (93) (for $\epsilon = -1$), in terms of the function $\tilde{a}_+$ (112) (instead of in terms of $a_+$ (84) for $\epsilon = -1$).<fnref n="3" />

Finally, the overlap is given by the expression (100) with $\epsilon = -1$.

<footnote n="3">³To compute the ratio of transfer matrix eigenvalues $\chi(\lambda_j)$, it may be convenient to isolate the contribution of the boundary roots, and to compute separately $\chi(\lambda_{BR}^-)$.</footnote>
**Case with $\lambda_{BR}^+$ and $\mu_{BR}^-$:** We finally consider the case in which the set of Bethe roots $\{\lambda\}$ contains the boundary root $\lambda_{BR}^+$, and the set $\{\mu\}$ the boundary root $\mu_{BR}^-$. This occurs when

$$h_1^- < h^+ < h_2^-, \quad \text{with} \quad |h^+|, |h_2^-| < h_{cr}^{(1)}, \quad \text{or} \quad |h^+| < h_{cr}^{(1)} < h_{cr}^{(2)} < h_1^- \quad (113)$$

Let us consider in that case the ratio of transfer matrices eigenvalues (68) evaluated at $\lambda_{BR}^+$:

$$\chi(\lambda_{BR}^{+}) = \frac{a_2(\lambda_{BR}^{+}|\{\mu\}) - 1 \sin(\lambda_{BR}^{+} - i\xi_2^{-} - i\zeta/2)}{a_1'(\lambda_{BR}^{+}|\{\lambda\}) \sin(\lambda_{BR}^{+} - i\xi_1^{-} - i\zeta/2)} \frac{\phi'(\lambda_{BR}^{+})}{\phi(\lambda_{BR}^{+} - i\zeta)} \quad (114)$$

As mentioned in (57), $a_2(\lambda|\{\mu\})$ diverges exponentially with the system size if $-\zeta < \mathfrak{J}(\lambda) < 0$. However, at $\lambda = \lambda_{BR}^+ = -i\xi_+ - i\zeta/2$, this divergence is compensated by the presence of the factor $\sin(\lambda + i\xi_+ + i\zeta/2)$, which appears in the Bethe equations for both Hamiltonians $H_1$ and $H_2$. It also follows from (64) that $a_1'(\lambda_{BR}^+|\{\lambda\}) = O(L^{+\infty})$. Hence, $\chi(\lambda_{BR}^+)$ is of order $O(L^{-\infty})$.

Proceeding with the computations similarly to what has been done in previous sections, we can show that all the other factors in (70) remain finite. Hence, the overlap is of order $O(L^{-\infty})$ in this case.

### 4.2.2 Presence of only one boundary root $\lambda_{BR}^{\sigma_1}$

We finally consider the case in which only one set of Bethe roots, say $\{\lambda\}$, contains a boundary root $\lambda_{BR}^{\sigma_1}$, the other set of Bethe roots $\{\mu\}$ involving only real roots. This occurs either when

$$h^+ < h_1^-, \quad \text{with} \quad h^+ < h_{cr}^{(1)} < h_2^- < h_{cr}^{(2)}, \quad \text{and} \quad h_1^- \in ]-h_{cr}^{(1)}, h_{cr}^{(1)}[ \cup ]h_{cr}^{(2)}, +\infty[, \quad (115)$$

in which case the boundary root is $\lambda_{BR}^-$, or when

$$h_1^- < h^+, \quad \text{with} \quad |h^+| < h_{cr}^{(1)} < h_2^- < h_{cr}^{(2)}, \quad (116)$$

in which case the boundary root is $\lambda_{BR}^+$. Note that in both cases we have

$$|p_1 q| < 1, \quad |(p_2 q)^{-1}| < 1, \quad \text{and} \quad q^2(p_1^{\sigma_1} q)^{\pm 1} < 1, \quad (117)$$

in which we have defined $p_1^{\sigma_1} = e^{-2i\xi_1^{\sigma_1}}$.

When considering a product of the form (71), we can still write

$$F_f(v|\{\lambda\}, \{\mu\}) = F_{f,BR}(v) \exp \left[ -\frac{1}{\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \log f(v-w) \hat{\xi}_d'(w) dw \right] [1 + O(L^{-\infty})], \quad (118)$$

with an extra factor $F_{f,BR}$ given by

$$F_{f,BR}(v) = f(v + i\frac{\zeta}{2} + i\xi_1^{\sigma_1}) f(v - i\frac{\zeta}{2} - i\xi_1^{\sigma_1}), \quad (119)$$

and the difference of the two counting function $\hat{\xi}_d$ satisfies the equation

$$\hat{\xi}_d'(\lambda) - \frac{1}{2\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \theta(\lambda - \mu) \hat{\xi}_d'(\mu) d\mu = \frac{g_2(\lambda) - g_1(\lambda)}{2} - \frac{1}{2} \sum_{\sigma=\pm} \theta(\lambda - i\sigma(\frac{\zeta}{2} + \xi_1^{\sigma_1})). \quad (120)$$

It is convenient, in that case, to consider the derivative of (120),

$$\hat{\xi}_d'(\lambda) + \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} K(\lambda - \mu) \hat{\xi}_d'(\mu) d\mu = \frac{g_2'(\lambda) - g_1'(\lambda)}{2} - \frac{1}{2} \sum_{\sigma=\pm} \theta'(\lambda - i\sigma(\frac{\zeta}{2} + \xi_1^{\sigma_1})), \quad (121)$$
$\hat{\xi}_d$ being the primitive function of the solution $\xi'_d$ of (121) which vanishes at 0. Note that the right hand side of (121) is a $2\pi$-periodic even function of $\lambda$. Hence, using similar arguments as in section 4.1.1, we can look for the solution of (121) as a Fourier series in the form

$$ \hat{\xi}'_d(\lambda) = \hat{\xi}'_0 + \hat{\xi}'_+(\lambda) + \hat{\xi}'_-(\lambda), \quad (122) $$

in which $\hat{\xi}'_0$ is a constant and $\hat{\xi}'_+(\lambda)$ is a power series of $e^{2i\lambda}$ with only positive powers, and (121) can be rewritten as

$$ \hat{\xi}'_0 + \hat{\xi}'_+(\lambda) + \hat{\xi}'_+(\lambda + i\zeta) = \frac{1}{2} \left[ \frac{1 + e^{2i\lambda}p_1q}{1 - e^{2i\lambda}p_1q} + \frac{1 + e^{2i\lambda}(p_2q)^{-1}}{1 - e^{2i\lambda}(p_2q)^{-1}} + \sum_{\sigma=\pm 1} \frac{1 + e^{2i\lambda}q^2(p_1^{\sigma_1}q)^\sigma}{1 - e^{2i\lambda}q^2(p_1^{\sigma_1}q)^\sigma} \right]. \quad (123) $$

It follows from (123) that $\hat{\xi}'_0 = 2$. Hence

$$ \hat{\xi}_d(\lambda) = 2\lambda + \hat{\xi}_+(\lambda) - \hat{\xi}_+(\lambda), \quad (124) $$

in which $\hat{\xi}_+$ is any primitive function of $\hat{\xi}'_+$. Up to a constant which can conveniently be chosen to be zero, the equation (123) can be integrated as

$$ \hat{\xi}_+(\lambda) + \hat{\xi}_+(\lambda + i\zeta) = -\frac{1}{2i} \log \left[ (1 - e^{2i\lambda}p_1q)(1 - e^{2i\lambda}(p_2q)^{-1}) \prod_{\sigma=\pm 1} (1 - e^{2i\lambda}q^2(p_1^{\sigma_1}q)^\sigma) \right], \quad (125) $$

or equivalently

$$ a_+(u) a_+(uq^2) = \frac{1}{1 - up_1q} \frac{1}{1 - u(p_2q)^{-1}} \frac{1}{1 - u(p_1^{\sigma_1})^{-1}q} \frac{1}{1 - up_1^{\sigma_1}q^3}, \quad (126) $$

in which we have defined $a_+(u) = \exp[2i\hat{\xi}_+(\lambda)]$, $u = e^{2i\lambda}$.

In particular, if $f$ is of the form (86), we have

$$ \exp\left[-\frac{1}{\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \log f_{p,\pm}(v-w) \hat{\xi}'_d(w) dw\right] = \exp\left[-\frac{1}{\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{f'_{p,\pm}(v-\mu)}{f_{p,\pm}(v-\mu)} \left[\hat{\xi}_+(\mu) - \hat{\xi}_+(\mu)\right] d\mu\right] \\ = a_+(pu^{\pm 1}), \quad (127) $$

in which we have used the decomposition (124)-(126) of $\hat{\xi}_d$.

**Case with $\lambda_{BR}^{-}$:** Let us consider more particularly the case in which the boundary root is $\lambda_{BR}^{-}$, i.e. the case (115). Then $p_1^{\sigma_1} = p_1$, and the unique solution of the functional equation (126) is

$$ a_+(u) = \frac{(up_2^{-1}q; q^4)_\infty}{(up_2^{-1}q^{-1}; q^4)_\infty} \frac{(up_1^{-1}q^3; q^4)_\infty}{(up_1^{-1}q; q^4)_\infty} \frac{1}{1 - up_1q}. \quad (128) $$

It follows that

$$ \tilde{F}_{f_{p,\pm},\text{BR}}(\lambda) \exp \left[ -\frac{1}{\pi} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \frac{f'_{p,\pm}(\lambda-\mu)}{f_{p,\pm}(\lambda-\mu)} \hat{\xi}_d(\mu) d\mu \right] = \tilde{a}_+(pu^{\pm 1}), \quad (129) $$

with $\tilde{a}_+(u)$ given by (112) and coinciding with the function obtained in (84) for $\epsilon = -1$.

The remaining part of the computation is then similar to what has been done in the other cases, and the overlap can be expressed as (100) for $\epsilon = -1$.
**Case with $\lambda_{BR}^+$:** In the case in which the boundary root is $\lambda_{BR}^+$, i.e. the case (116), one can apply similar arguments as for (113), and one finds that the overlap vanishes up to exponentially small corrections in $L$.

# 5 The overlap in the thermodynamic limit: Summary of the results

We now summarize our final results concerning the value of the overlap in the thermodynamic limit, according to the different configurations of the boundary magnetic fields.

## 5.1 The overlap for a chain with an odd number of sites

The study of the thermodynamic limit of the overlap for $L$ odd is quite simple, since we have only very few cases to distinguish, see section 2.3:

1. If $h^+, h_1^-, h_2^- < h_{cr}^{(1)}$ and $h_1^-, h_2^- < -h_+$, all the Bethe roots are real and the overlap is given by the expression (100) with $\epsilon = 1$:

$$S(\{\lambda\}, \{\mu\}) = \frac{(p_1^2 q^2; q^4, q^4)_\infty (p_2^2 q^2; q^4, q^4)_\infty (p_1 p_2 q^4; q^4, q^4)_\infty^2}{(p_1^2 q^4; q^4, q^4)_\infty (p_2^2 q^4; q^4, q^4)_\infty (p_1 p_2 q^2; q^4, q^4)_\infty^2} + O(L^{-\infty}). \quad (130)$$

2. If $|h^+|, h_1^-, -h_2^- < h_{cr}^{(1)}$ and $h_1^- < -h_+ < h_2^-$, the two ground states are in sectors of different magnetisation, so that the overlap vanishes identically:

$$S(\{\lambda\}, \{\mu\}) = 0. \quad (131)$$

3. If $h^+, h_1^-, h_2^- > -h_{cr}^{(1)}$ and $h_1^-, h_2^- > -h_+$, the overlap can be obtained by spin-reversal symmetry from the first case, using that

$$(\sigma^x)^{\otimes N} \mathbf{H}_{h^-, h^+} (\sigma^x)^{\otimes N} = \mathbf{H}_{-h^-, -h^+}. \quad (132)$$

Hence the overlap in this case is simply given by the expression (130) in which we have replaced $p_i$ by $p_i^{-1}$:

$$S(\{\lambda\}, \{\mu\}) = \frac{(p_1^{-2}q^2; q^4, q^4)_\infty (p_2^{-2}q^2; q^4, q^4)_\infty ((p_1p_2)^{-1}q^4; q^4, q^4)_\infty^2}{(p_1^{-2}q^4; q^4, q^4)_\infty (p_2^{-2}q^4; q^4, q^4)_\infty ((p_1p_2)^{-1}q^2; q^4, q^4)_\infty^2} + O(L^{-\infty}). \quad (133)$$

We recall that in all these expressions we have defined $p_i = e^{-2\xi_i^-}$, $i = 1, 2$.

In Fig. 2, we have plotted our analytical result and compared it with numerical results obtained by exact diagonalisation of the Hamiltonian using the Quspin package [26].
<figure><img src="image_6.png" /><figcaption>Figure 2: The overlap via exact diagonalisation for a chain of length $L = 17$ compared to the ABA exact result at the thermodynamic limit obtained in section 5.1. Here $\zeta = 1.8$, $h_+ = -1$, $h_1^- = 0$, and the value of the overlap $S(\{\lambda\}, \{\mu\})$ is plotted for different values of $h_2^-$: when $h_2^- < -h_+ = 1$, we are in the case 1 of section 5.1, and the overlap is given by (130); when $h_2^- > -h_+ = 1$, we are in the case 2 of section 5.1 and the overlap vanishes, see (131).</figcaption></figure>
## 5.2 The overlap for a chain with an even number of sites

The analytical study in this case is slightly more complicated than for $L$ odd. Indeed, we have a priori to consider many different configurations for the sets of Bethe roots describing the ground states, according to whether they contain or not a boundary root, see section 2.3. We can nevertheless distinguish the following three main different cases:

1. Case $h_1^-, h_2^- < h^+$. Under the hypothesis that the two ground states are in a configuration (A), (B) or (C), see section 2.3, this may occur in the following situations:

(i) $h_1^-, h_2^- < h^+$ with $|h^+| < h_{cr}^{(1)}$ or $h_1^-, h_2^- < h_{cr}^{(1)} < h_{cr}^{(2)} < h^+$, i.e. we are in the situation considered in (105).

(ii) $h_1^-, h_2^- < h_{cr}^{(1)} < h^+ < h_{cr}^{(2)}$, i.e. we are in the situation (b1) considered in section 4.1.

For all these situations, the overlap is given by the expression (100) with $\epsilon = 1$, i.e.

$$S(\{\lambda\}, \{\mu\}) = \frac{(p_1^2 q^2; q^4, q^4)_{\infty} (p_2^2 q^2; q^4, q^4)_{\infty} (p_1 p_2 q^4; q^4, q^4)_{\infty}^2}{(p_1^2 q^4; q^4, q^4)_{\infty} (p_2^2 q^4; q^4, q^4)_{\infty} (p_1 p_2 q^2; q^4, q^4)_{\infty}^2} + O(L^{-\infty}). \quad (134)$$

2. Case $h_1^- < h^+ < h_2^-$. Under the hypothesis that the two ground states are in a configuration (A), (B) or (C), we may be in one of the two following situations:

(i) $|h^+|, |h_2^-| < h_{cr}^{(1)}$, or $|h^+| < h_{cr}^{(1)} < h_{cr}^{(2)} < h_2^-$, i.e. we are in the situation considered in (113).

(ii) $|h^+| < h_{cr}^{(1)} < h_2^- < h_{cr}^{(2)}$, i.e. we are in the situation considered in (116).

Then the overlap vanishes up to exponentially small corrections in $L$:

$$S(\{\lambda\}, \{\mu\}) = O(L^{-\infty}). \quad (135)$$

3. Case $h^+ < h_1^-, h_2^-$. The overlap can be obtained by spin-reversal symmetry from the first case, i.e.

$$S(\{\lambda\}, \{\mu\}) = \frac{(p_1^{-2}q^2; q^4, q^4)_{\infty} (p_2^{-2}q^2; q^4, q^4)_{\infty} ((p_1p_2)^{-1}q^4; q^4, q^4)_{\infty}^2}{(p_1^{-2}q^4; q^4, q^4)_{\infty} (p_2^{-2}q^4; q^4, q^4)_{\infty} ((p_1p_2)^{-1}q^2; q^4, q^4)_{\infty}^2} + O(L^{-\infty}). \quad (136)$$

Note that this case corresponds also to several situations in which the overlap has been explicitly computed, namely the situation considered in (106), or in (115), or in case (b2) of section 4.1. For all these situations, the overlap is given by the expression (100) with $\epsilon = -1$, which indeed coincides with (136).

In figures 3 and 4, we compare the analytic expression of the overlap we have obtained in several of these different cases with numerical results obtained by exact diagonalisation using the Quspin package [26].

Note that, as expected, the thermodynamic limit of the overlap does not depend on the parity of the length of the chain, up to a change of $h^+$ into $-h^+$. We would also like to mention that our final expression for the overlap in the thermodynamic limit coincides with a similar result obtained by R. Weston within the framework of the q-operator approach [17].
<figure><img src="image_7.png" /><figcaption>Figure 3: The overlap via exact diagonalisation for a chain of length $L = 18$ compared to the ABA exact result at the thermodynamic limit obtained in section 5.2. Here $\zeta = 1.5$, $h_+ = 2$, $h_1^- = -1$, and the value of the overlap $S(\{\lambda\}, \{\mu\})$ is plotted for different values of $h_2^-$ for $h_2^- < h_{cr}^{(1)}$: we are here in the configuration (ii) of the case 1 from section 5.2, and the overlap is given by (134).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4: The overlap via exact diagonalisation for a chain of length $L = 18$ compared to the exact result at the thermodynamic limit obtained in section 5.2. Here $\zeta = 1.8$, $h_+ = 0$, $h_2^- = 1$, and the value of the overlap $S(\{\lambda\}, \{\mu\})$ is plotted for different values of $h_1^-$: when $h_1^- < h_+$, we are in configuration (i) of case 2 of the section 5.2 (up to an exchange of $h_1^-$ and $h_2^-$), and the overlap vanishes up to exponentially small corrections in $L$; when $h_1^- > h_+$, we are in case 3 of section 5.2, and the overlap is given by (136).</figcaption></figure>
# 6 Conclusion

We have considered a boundary quench in the open XXZ spin chain with boundary magnetic fields parallel to the anisotropy axis, i.e. a change of the value of one of the boundary magnetic fields, and computed the overlaps between the ground states before and after the quench. Our approach is based on the Slavnov determinant representation of the overlaps, and on the Gaudin extraction technique proposed in [14]. In the massive antiferromagnetic regime of the chain $\Delta > 1$, and for all configurations of the magnetic fields for which the spectrum remains gapped (i.e. for which the ground state solutions of the Bethe equations do not include real holes), we have computed the thermodynamic limit of these overlaps up to exponentially small corrections in the length $L$ of the chain.

The fact that we limited our consideration to the massive antiferromagnetic regime enabled us to avoid issues with the convergence of infinite products. It seems however to be possible to apply similar technics in the massless case. In particular, the XXX case can directly be obtained from the massive one as a limit.

In our next publication, we intend to show how to generalise this computation to excited states, and namely how to deal with the presence of real holes. This should effectively open the way to a study of the quench dynamics, or of the behaviour of boundary driven spin chains. Further development of this approach should include all the excited states close to the ground state which implies treatment of the complex roots of the Bethe equations as it was done for the periodic XXX chain [27].

# Acknowledgments

We would like to thank R. Weston for interesting discussions, and in particular for informing us that these kind of results can also be obtained in the framework of the *q*-vertex operator approach [17]. N.K. would like to thank the LPTMS laboratory (Université Paris-Saclay) for hospitality which made this collaboration possible. N.K. is also grateful to the LPTHE laboratory for hospitality.

**Funding information** The institution of N.K. (IMB) receives support from the EIPHI Graduate School (contract ANR-17-EURE-0002). V.T. is supported by CNRS.

# A Convergence of numerical results

In this Appendix we illustrate by two plots and corresponding value tables the rapid convergence of numerical results for the overlap toward our analytic formulas (with and without boundary roots).
<figure><img src="image_9.png" /></figure>

<table><caption>Figure 5: The overlap via exact diagonalization for different chain sizes compared to the exact result from ABA at the thermodynamic limit. Here, the values of the overlap are plotted for $\zeta = 1.5$, $h_+ = 2$ and $h_1^- = -1$. They are also presented in a table to show the rapid (exponential) convergence of the numerical results to the analytic values. In this plot, we are in the exact same configuration as in figure 3.</caption><thead><tr><th>h2-</th><th>L=4</th><th>L=8</th><th>L=12</th><th>L=16</th><th>L=20</th><th>Analytic</th></tr></thead><tbody><tr><td>-5.000</td><td>0.9735854</td><td>0.9727717</td><td>0.9726791</td><td>0.9726645</td><td>0.9726618</td><td>0.9726611</td></tr><tr><td>-4.711</td><td>0.9753404</td><td>0.9745309</td><td>0.9744377</td><td>0.9744230</td><td>0.9744202</td><td>0.9744195</td></tr><tr><td>-4.421</td><td>0.9771818</td><td>0.9763802</td><td>0.9762869</td><td>0.9762720</td><td>0.9762692</td><td>0.9762685</td></tr><tr><td>-4.132</td><td>0.9791114</td><td>0.9783226</td><td>0.9782294</td><td>0.9782145</td><td>0.9782117</td><td>0.9782110</td></tr><tr><td>-3.842</td><td>0.9811292</td><td>0.9803593</td><td>0.9802669</td><td>0.9802520</td><td>0.9802492</td><td>0.9802484</td></tr><tr><td>-3.553</td><td>0.9832324</td><td>0.9824892</td><td>0.9823985</td><td>0.9823836</td><td>0.9823808</td><td>0.9823801</td></tr><tr><td>-3.263</td><td>0.9854142</td><td>0.9847076</td><td>0.9846194</td><td>0.9846048</td><td>0.9846021</td><td>0.9846013</td></tr><tr><td>-2.974</td><td>0.9876613</td><td>0.9870035</td><td>0.9869193</td><td>0.9869052</td><td>0.9869024</td><td>0.9869017</td></tr><tr><td>-2.684</td><td>0.9899506</td><td>0.9893565</td><td>0.9892781</td><td>0.9892648</td><td>0.9892622</td><td>0.9892615</td></tr><tr><td>-2.395</td><td>0.9922443</td><td>0.9917316</td><td>0.9916615</td><td>0.9916493</td><td>0.9916469</td><td>0.9916463</td></tr><tr><td>-2.105</td><td>0.9944822</td><td>0.9940703</td><td>0.9940116</td><td>0.9940011</td><td>0.9939990</td><td>0.9939985</td></tr><tr><td>-1.816</td><td>0.9965692</td><td>0.9962769</td><td>0.9962331</td><td>0.9962251</td><td>0.9962234</td><td>0.9962230</td></tr><tr><td>-1.526</td><td>0.9983570</td><td>0.9981953</td><td>0.9981695</td><td>0.9981646</td><td>0.9981636</td><td>0.9981633</td></tr><tr><td>-1.237</td><td>0.9996138</td><td>0.9995693</td><td>0.9995617</td><td>0.9995601</td><td>0.9995598</td><td>0.9995597</td></tr><tr><td>-0.947</td><td>0.9999777</td><td>0.9999746</td><td>0.9999740</td><td>0.9999739</td><td>0.9999739</td><td>0.9999739</td></tr><tr><td>-0.658</td><td>0.9988822</td><td>0.9986971</td><td>0.9986584</td><td>0.9986495</td><td>0.9986474</td><td>0.9986467</td></tr><tr><td>-0.368</td><td>0.9954430</td><td>0.9945147</td><td>0.9942927</td><td>0.9942366</td><td>0.9942220</td><td>0.9942167</td></tr><tr><td>-0.079</td><td>0.9882915</td><td>0.9853062</td><td>0.9844650</td><td>0.9842224</td><td>0.9841519</td><td>0.9841227</td></tr><tr><td>0.211</td><td>0.9753640</td><td>0.9673904</td><td>0.9646510</td><td>0.9637123</td><td>0.9633942</td><td>0.9632346</td></tr><tr><td>0.500</td><td>0.9537220</td><td>0.9346186</td><td>0.9263587</td><td>0.9228312</td><td>0.9213589</td><td>0.9203628</td></tr></tbody></table>
<figure><img src="image_10.png" /></figure>

<table><caption>Figure 6: The overlap via exact diagonalization for different chain sizes compared to the exact result from ABA at the thermodynamic limit. Here, the values of the overlap are plotted for $\zeta = 1.8$, $h_+ = 0$ and $h_1^- = 1$. They are also presented in a table to show the rapid (exponential) convergence of the numerical results to the analytic values. In this plot, we are in the same exact configuration as in figure 4.</caption><thead><tr><th>h2-</th><th>L=4</th><th>L=8</th><th>L=12</th><th>L=16</th><th>L=20</th><th>Analytic (ABA)</th></tr></thead><tbody><tr><td>-2.000</td><td>0.352847</td><td>0.061289</td><td>0.008120</td><td>0.001027</td><td>0.000129</td><td>0.000000</td></tr><tr><td>-1.724</td><td>0.383196</td><td>0.067860</td><td>0.009015</td><td>0.001141</td><td>0.000143</td><td>0.000000</td></tr><tr><td>-1.448</td><td>0.422264</td><td>0.077293</td><td>0.010321</td><td>0.001307</td><td>0.000164</td><td>0.000000</td></tr><tr><td>-1.172</td><td>0.473265</td><td>0.091814</td><td>0.012388</td><td>0.001571</td><td>0.000198</td><td>0.000000</td></tr><tr><td>-0.897</td><td>0.540083</td><td>0.116480</td><td>0.016088</td><td>0.002047</td><td>0.000257</td><td>0.000000</td></tr><tr><td>-0.621</td><td>0.625784</td><td>0.164779</td><td>0.024234</td><td>0.003112</td><td>0.000392</td><td>0.000000</td></tr><tr><td>-0.345</td><td>0.728218</td><td>0.279204</td><td>0.051331</td><td>0.006885</td><td>0.000872</td><td>0.000000</td></tr><tr><td>-0.069</td><td>0.833767</td><td>0.564982</td><td>0.304849</td><td>0.086419</td><td>0.013293</td><td>0.000000</td></tr><tr><td>0.207</td><td>0.919872</td><td>0.881636</td><td>0.955694</td><td>0.989780</td><td>0.995323</td><td>0.996146</td></tr><tr><td>0.483</td><td>0.972496</td><td>0.978721</td><td>0.994855</td><td>0.998123</td><td>0.998558</td><td>0.998620</td></tr><tr><td>0.759</td><td>0.995370</td><td>0.997620</td><td>0.999409</td><td>0.999701</td><td>0.999738</td><td>0.999744</td></tr><tr><td>1.034</td><td>0.999928</td><td>0.999971</td><td>0.999992</td><td>0.999995</td><td>0.999995</td><td>0.999995</td></tr><tr><td>1.310</td><td>0.995482</td><td>0.998394</td><td>0.999500</td><td>0.999658</td><td>0.999678</td><td>0.999681</td></tr><tr><td>1.586</td><td>0.987357</td><td>0.995837</td><td>0.998566</td><td>0.998945</td><td>0.998991</td><td>0.998998</td></tr><tr><td>1.862</td><td>0.978157</td><td>0.993113</td><td>0.997412</td><td>0.997997</td><td>0.998069</td><td>0.998079</td></tr><tr><td>2.138</td><td>0.969047</td><td>0.990464</td><td>0.996144</td><td>0.996906</td><td>0.996998</td><td>0.997011</td></tr><tr><td>2.414</td><td>0.960501</td><td>0.987962</td><td>0.994821</td><td>0.995731</td><td>0.995841</td><td>0.995856</td></tr><tr><td>2.690</td><td>0.952674</td><td>0.985622</td><td>0.993479</td><td>0.994513</td><td>0.994638</td><td>0.994655</td></tr><tr><td>2.966</td><td>0.945578</td><td>0.983437</td><td>0.992142</td><td>0.993280</td><td>0.993417</td><td>0.993435</td></tr><tr><td>3.241</td><td>0.939170</td><td>0.981398</td><td>0.990825</td><td>0.992050</td><td>0.992197</td><td>0.992217</td></tr><tr><td>3.517</td><td>0.933383</td><td>0.979491</td><td>0.989537</td><td>0.990837</td><td>0.990993</td><td>0.991014</td></tr><tr><td>3.793</td><td>0.928147</td><td>0.977704</td><td>0.988286</td><td>0.989649</td><td>0.989812</td><td>0.989834</td></tr><tr><td>4.069</td><td>0.923398</td><td>0.976026</td><td>0.987073</td><td>0.988492</td><td>0.988661</td><td>0.988684</td></tr><tr><td>4.345</td><td>0.919076</td><td>0.974449</td><td>0.985902</td><td>0.987369</td><td>0.987544</td><td>0.987567</td></tr><tr><td>4.621</td><td>0.915131</td><td>0.972962</td><td>0.984774</td><td>0.986282</td><td>0.986462</td><td>0.986486</td></tr><tr><td>4.897</td><td>0.911516</td><td>0.971559</td><td>0.983687</td><td>0.985232</td><td>0.985416</td><td>0.985441</td></tr><tr><td>5.172</td><td>0.908194</td><td>0.970232</td><td>0.982642</td><td>0.984220</td><td>0.984407</td><td>0.984433</td></tr><tr><td>5.448</td><td>0.905132</td><td>0.968975</td><td>0.981637</td><td>0.983244</td><td>0.983435</td><td>0.983461</td></tr><tr><td>5.724</td><td>0.902299</td><td>0.967784</td><td>0.980671</td><td>0.982305</td><td>0.982499</td><td>0.982525</td></tr><tr><td>6.000</td><td>0.899672</td><td>0.966652</td><td>0.979744</td><td>0.981401</td><td>0.981597</td><td>0.981624</td></tr></tbody></table>
References

[1] E. K. Sklyanin, Boundary conditions for integrable quantum systems, J. Phys. A: Math. Gen. 21, 2375 (1988), doi:10.1088/0305-4470/21/10/015.

[2] S. Grijalva, J. De Nardis and V. Terras, Open XXZ chain and boundary modes at zero temperature, SciPost Phys. 7, 023 (2019), doi:10.21468/SciPostPhys.7.2.023.

[3] P. R. Pasnoori, J. Lee, J. H. Pixley, N. Andrei and P. Azaria, Boundary quantum phase transitions in the spin-1/2 Heisenberg chain with boundary magnetic fields, Phys. Rev. B 107, 224412 (2023), doi:10.1103/PhysRevB.107.224412.

[4] J.-S. Caux and F. H. L. Essler, Time evolution of local observables after quenching to an integrable model, Phys. Rev. Lett. 110, 257203 (2013), doi:10.1103/PhysRevLett.110.257203.

[5] M. Brockmann, J. De Nardis, B. Wouters and J.-S. Caux, A Gaudin-like determinant for overlaps of Néel and XXZ Bethe states, J. Phys. A: Math. Theor. 47, 145003 (2014), doi:10.1088/1751-8113/47/14/145003.

[6] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H. L. Essler and T. Prosen, Complete generalized Gibbs ensembles in an interacting theory, Phys. Rev. Lett. 115, 157201 (2015), doi:10.1103/PhysRevLett.115.157201.

[7] L. Piroli, B. Pozsgay and E. Vernier, What is an integrable quench?, Nucl. Phys. B 925, 362 (2017), doi:10.1016/j.nuclphysb.2017.10.012.

[8] I. V. Cherednik, Some finite-dimensional representations of generalized Sklyanin algebras, Funct. Anal. Appl. 19, 77 (1985), doi:10.1007/BF01086039.

[9] N. A. Slavnov, Calculation of scalar products of wave functions and form factors in the framework of the algebraic Bethe Ansatz, Theor. Math. Phys. 79, 502 (1989), doi:10.1007/BF01016531.

[10] N. Kitanine, J. M. Maillet and V. Terras, Form factors of the XXZ Heisenberg finite chain, Nucl. Phys. B 554, 647 (1999), doi:10.1016/S0550-3213(99)00295-3.

[11] Y.-S. Wang, The scalar products and the norm of Bethe eigenstates for the boundary XXX Heisenberg spin-1/2 finite chain, Nucl. Phys. B 622, 633 (2002), doi:10.1016/S0550-3213(01)00610-1.

[12] N. Kitanine, K. K. Kozlowski, J. M. Maillet, G. Niccoli, N. A. Slavnov and V. Terras, Correlation functions of the open XXZ chain: I, J. Stat. Mech.: Theory Exp. P10009 (2007), doi:10.1088/1742-5468/2007/10/P10009.

[13] P. R. Pasnoori, Y. Tang, J. Lee, J. H. Pixley, N. Andrei and P. Azaria, Spin fractionalization and zero modes in the spin-1/2 XXZ chain with boundary fields, (arXiv preprint) doi:10.48550/arXiv.2312.05970.

[14] N. Kitanine and G. V. Kulkarni, Thermodynamic limit of the two-spinon form factors for the zero field XXX chain, SciPost Phys. 6, 076 (2019), doi:10.21468/SciPostPhys.6.6.076.

[15] A. Kapustin and S. Skorik, Surface excitations and surface energy of the antiferromagnetic XXZ chain by the Bethe Ansatz approach, J. Phys. A: Math. Gen. 29, 1629 (1996), doi:10.1088/0305-4470/29/8/011.
[16] R. Weston, *Correlation functions and the boundary qKZ equation in a fractured XXZ chain*, J. Stat. Mech.: Theory Exp. P12002 (2011), doi:10.1088/1742-5468/2011/12/P12002.

[17] R. Weston, *Private communication* (2023).

[18] M. Jimbo and T. Miwa, *Algebraic analysis of solvable lattice models*, American Mathematical Society, Providence, USA, ISBN 9780821803202 (1994), doi:10.1090/cbms/085.

[19] M. Jimbo, R. Kedem, T. Kojima, H. Konno and T. Miwa, *XXZ chain with a boundary*, Nucl. Phys. B **441**, 437 (1995), doi:10.1016/0550-3213(95)00062-W.

[20] M. Jimbo, R. Kedem, H. Konno, T. Miwa and R. Weston, *Difference equations in spin chains with a boundary*, Nucl. Phys. B **448**, 429 (1995), doi:10.1016/0550-3213(95)00218-H.

[21] I. V. Cherednik, *Factorizing particles on a half-line and root systems*, Theor. Math. Phys. **61**, 977 (1984), doi:10.1007/BF01038545.

[22] F. C. Alcaraz, M. N. Barber, M. T. Batchelor, R. J. Baxter and G. R. W. Quispel, *Surface exponents of the quantum XXZ, Ashkin-Teller and Potts models*, J. Phys. A: Math. Gen. **20**, 6397 (1987), doi:10.1088/0305-4470/20/18/038.

[23] I. S. Gradshteyn and I. M. Ryzhik, *Table of integrals, series, and products*, Elsevier, Amsterdam, Netherlands, ISBN 9780123736376 (2007).

[24] O. Babelon, H. J. de Vega and C. M. Viallet, *Analysis of the Bethe Ansatz equations of the XXZ model*, Nucl. Phys. B **220**, 13 (1983), doi:10.1016/0550-3213(83)90131-1.

[25] N. Kitanine, J. M. Maillet and V. Terras, *Correlation functions of the XXZ Heisenberg spin-chain in a magnetic field*, Nucl. Phys. B **567**, 554 (2000), doi:10.1016/S0550-3213(99)00619-7.

[26] P. Weinberg and M. Bukov, *QuSpin: A Python package for dynamics and exact diagonalisation of quantum many body systems part I: Spin chains*, SciPost Phys. **2**, 003 (2017), doi:10.21468/SciPostPhys.2.1.003.

[27] N. Kitanine and G. Kulkarni, *Form factors of the Heisenberg spin chain in the thermodynamic limit: Dealing with complex Bethe roots*, Symmetry Integr. Geom.: Methods Appl. **17**, 112 (2021), doi:10.3842/SIGMA.2021.112.