Off-shell scalar products for the $XXZ$ spin chain with open boundaries

In this work we study scalar products of Bethe vectors associated with the $XXZ$ spin chain with open boundary conditions. The scalar products are obtained as solutions of a system of functional equations. The description of scalar products through functional relations follows from a particular map having the reflection algebra as its domain and a function space as the codomain. Within this approach we find a multiple contour integral representation for the scalar products in which the homogeneous limit can be obtained trivially.


Introduction
The inverse scattering method was formulated in the late sixties by Gardner, Greene, Kruskal and Miura [1] and it represented a large step towards the understanding of evolution equations. This method was originally conceived for solving the Korteweg-deVries (KdV) equation [2] but turned out to be a very general framework, capable of integrating a large number of non-linear differential equations. Although the application of this method is restricted to a special class of differential equations, its importance is twofold: on the one hand, it offers a fertile soil where novel and sophisticated mathematical structures have emerged and continue to do so; on the other hand, the inverse scattering provides a powerful tool for the exact analysis of physical quantities. For instance, although numerical studies of the KdV equation were already available in the sixties, the exact solution obtained in [1] elucidated several aspects concerning the solitonic behavior reported in [3]. The history repeated itself in [4] where the spin of a spin wave was finally understood through the exact solution of the XXX spin chain obtained by means of the Bethe ansatz [5]. The fundamental difference between the non-linear waves governed by the KdV equation and the spin waves of the XXX chain, is that the first is described by a partial differential equation while the latter arises from the spectrum of excitations of a quantum mechanical hamiltonian operator. The hamiltonian of the XXX model is not defined as a differential operator; hence its spectral problem is not a priori characterized by a differential equation. Nevertheless, the spectrum of the XXX spin chain was shown by Baxter to admit a description in terms of functional equations [6].
The existence of a direct correspondence between Baxter's functional equations and integrable non-linear differential equations is not apparent at first sight but several results suggest this relation should exist. This possibility has been concluded from different approaches and we refer the reader to [7][8][9][10][11][12][13] and references therein. In particular, in the work [13] we have presented a first principle mechanism able to convert the spectral problem of spin chains built from solutions of the Yang-Baxter equation [6,14,15] into the solution of linear partial differential equations. This mechanism is one of the outcomes of the Algebraic-Functional Method proposed in [16] and refined in the series of papers [17][18][19][20][21].
Among the results obtained through the algebraic-functional approach, we have demonstrated in [20] how scalar products of Bethe vectors can be described by functional equations. More precisely, in [20] we have shown how the Yang-Baxter algebra can be exploited in order to derive functional relations determining the scalar products associated with the XXZ chain with periodic boundary conditions. The algebraic structure underlying the XXZ chain with open boundary conditions is the so called reflection algebra [22] and here we intend to extend the results of [20] for that case. The reflection algebra is the analogue of the Yang-Baxter algebra for open spin chains and in [23] we have already demonstrated the feasibility of exploiting it along the lines of [16].
Scalar products of Bethe vectors are building blocks of correlation functions within the framework of the algebraic Bethe ansatz [24][25][26] and, as far as the XXX model with diagonal open boundaries is concerned, those scalar products have been previously studied in [27] using the method of [26]. The results of [27] are expressed as determinants and they have been generalized to the XXZ chain in [28]. Also, these results have been used in [29,28,30] to evaluate certain correlation functions for the XXX and XXZ models with open boundaries. In the works [28,30], in particular, the authors have obtained multiple integral representations for correlation functions in the case of a halfinfinite lattice. It is also important to remark here that for half-infinite lattices such correlation functions have been studied in [31,32] through the vertex-operator approach and in [33][34][35] using the q-Onsager algebra.
Here we use the terminology scalar product of Bethe vectors but we should keep in mind the conditions in which such quantities are computed. Bethe vectors are parameterized by a set of complex variables which are required to satisfy a set of algebraic equations, i.e. Bethe ansatz equations, in order to having a model eigenvector. The dual Bethe vector required to compute the aforementioned scalar products also carries another set of complex variables. In this way we have the so called off-shell scalar products when both sets of variables parameterizing the scalar product are arbitrary complex variables. On the other hand, the cases where only one or both sets of variables satisfy Bethe ansatz equations are usually refereed to as on-shell scalar products. That being said, we should also remark here that the scalar products obtained in [27,28] are on-shell. The case of off-shell scalar products for open spin chains has not been considered in the literature so far to the best of our knowledge, and we shall address this problem in the present paper. Although the computation of correlation functions usually require on-shell scalar prod-ucts, the off-shell case can still be regarded as a partition function of a vertex model with special boundary conditions [24,36] as discussed in [37]. More precisely, here we shall present a multiple integral representation for the off-shell scalar product of the XXZ spin chain with open boundary conditions. Interestingly, this integral representation allows one to obtain the homogeneous limit in a trivial way.
The XXZ spin-chain with open boundaries has been extensively discussed in the literature and this paper has been organized in such a manner to avoid large overlaps with the existing literature. Therefore, we shall simply collect the required definitions in Section 2 in order to clarify our notation. Section 3.1 is then devoted to the derivation of functional relations characterizing the desired scalar products. The solution of our functional equations is then obtained in Section 4. Concluding remarks are discussed in Section 5 and appendices A through E are devoted to technical details and proofs.

The open XXZ model and Bethe vectors
The anisotropic Heisenberg chain was proposed by Bloch in [38,39] as a model for remanent magnetization. The exact solution of the XXZ case was then obtained by Yang and Yang through Bethe ansatz in the series of papers [40][41][42][43]. In 1967 Sutherland [44] found a vertex model transfer matrix exhibiting the same eigenvectors firstly obtained by Yang and Yang for the XXZ spin chain. It is worth remarking here that, at that time, the relation between transfer matrices and spin chains solvable by Bethe ansatz was not clear yet. The aforementioned works considered the case with periodic boundary conditions while the case with parallel boundary fields was then solved in [45] through a generalization of the Bethe ansatz.
The model hamiltonian is a linear operator H : (C 2 ) ⊗L → (C 2 ) ⊗L with L ∈ Z >0 and it explicitly reads where J x = J y = 1 and J z = cosh (γ). Here h,h and γ are arbitrary complex parameters while σ x i , σ y i and σ z i are standard Pauli matrices acting on the i-th node of the tensor product space (C 2 ) ⊗L . The hamiltonian (2.1) can be embedded as the derivative of a commuting transfer matrix at a particular point. This construction will be discussed in what follows and here we shall also consider the definitions and notation employed in [23]. Moreover, we shall restrict ourselves to presenting only the required definitions which have not been described in [23].
Reflection matrices. Let K,K : C → End(C 2 ) be respectively refereed to as reflection matrix and dual reflection matrix. Then following [22], the construction of integrable spin chains with open boundary conditions through the Quantum Inverse Scattering Method (QISM) requires K andK to satisfy the so called reflection equations. In this way, K is governed by the equation whileK is required to satisfy the following dual relation Both relations (2.2) and (2.3) involve the operator-valued functions R ij : Here we have V i ∼ = C 2 as we are considering the XXZ spin chain hamiltonian (2.1). Also, the operator R 12 corresponds to the U q [ sl (2)] invariant R-matrix satisfying the Yang-Baxter equation as described in [23]. In addition to that, the symbol t i stands for the standard transposition on the space V i of a generic matrix in End(V i ⊗ V j ).
Here we shall consider the following solution of (2.2), where h ∈ C corresponds to an arbitrary parameter governing the field strength at one of the boundaries of the hamiltonian (2.1). The dual reflection matrixK is obtained with the help of the following lemma.
Proof. We apply the map d to (2.3) taking into account the following properties satisfied by the U q [ sl(2)] R-matrix: R 21 = R 12 , R t 1 12 = R t 2 12 , R 12 (λ)R 12 (−λ) ∝ id ⊗ id and R t 1 12 (λ)R t 1 12 (−λ − 2γ) ∝ id ⊗ id. Lemma 1 establishes a map d : K →K which can be exploited to determine a solution of (2.3) in a straightforward manner. From (2.4) we then find the dual reflection matrix Double-row transfer matrix. Let λ, µ j ∈ C be respectively refereed to as spectral parameter and inhomogeneity parameter. Then the double-row transfer matrix or simply transfer matrix T : C → End((C 2 ) ⊗L ) is defined as the following operator, The hamiltonian (2.1) then corresponds to (2.7) ABCD structure. The trace tr 0 in (2.6) is taken only over the space V 0 , while the term inside the brackets lives in End(V 0 ⊗ V Q ) where V Q ∼ = (C 2 )) ⊗L . Thus, it is convenient to employ the following representation, where A, B, C, D ∈ End(V Q ). In terms of these operators, the transfer matrix (2.6) reads Highest-weight property. The highest-weight representation theory of the sl(2) Lie algebra plays a prominent role for the system we are considering. For instance, the vector is a sl(2) highest-weight vector in V Q . Its dual is simply given by 0| := (1 0) ⊗L . The left and right action of the operators (2.8) on these vectors can be straightforwardly computed due to the structure of the U q [ sl (2)] invariant R-matrix entering in the definition (2.8).
In this way we have the following expressions: whereD(λ) := D(λ) − sinh (γ) sinh (2λ+γ) A(λ) has been defined for latter convenience. It is also useful to define the functions a(λ) := sinh (λ + γ), b(λ) := sinh (λ) and c(λ) := sinh (γ), in such a way that the coefficients Λ A and ΛD explicitly read Bethe vectors. One of the achievements of the QISM is the algebraic construction of eigenvectors of commuting transfer matrices. This is also the case for spin chains with open boundary conditions as shown in [22]. In particular, the eigenvectors associated with the spin chain hamiltonian (2.1) can be built with the help of the operators defined in (2.8). Those vectors |Ψ n ∈ span(V Q ) will be refereed to as off-shell Bethe vectors and they are defined as Strictly speaking, the vector |Ψ n is only an eigenvector of the transfer matrix (2.9) for particular choices of the set of variables {λ B j }. Dual eigenvectors can also be built in a similar way and we shall refer to them as Dual Bethe vectors. They are defined as (2.14) Here we shall assume that the set of variables {λ C j } is independent of the set {λ B j } in order to keep our results as general as possible.
Scalar product. The quantity we are interested in the present work is the scalar product of Bethe vectors S n : C 2n → C defined as, The quantity S n is usually refereed to as off-shell scalar product as the variables λ B,C j are arbitrary. The on-shell case corresponds to the situation where the variables λ B j , and/or λ C j , are constrained by Bethe ansatz equations. Here we are interested in the off-shell case, but for completeness we shall also present the on-shell constraint. It reads , (2.16) where θ(λ, ω) := sinh (λ+ω) sinh (λ−ω+γ) . The solutions of (2.16) are usually refereed to as Bethe roots. Also, dual Bethe vectors (2.14) will only be eigenvectors of (2.9) when the variables λ C j are Bethe roots. The scalar product S n is a key ingredient for the computation of correlation functions [27,28,30] and in what follows we shall investigate this quantity along the lines of [20].

Reflection algebra and functional equations
The operators A, B, C and D defined in (2.8) are subjected to the reflection algebra relations [22]. See also [23] for the conventions we are using here. This property is a direct consequence of the Yang-Baxter equation [46], in addition to the reflection equation (2.2). In the present paper we aim to use the reflection algebra relations as a source of functional equations characterizing the scalar product (2.15). For that it is convenient to introduce the following definitions.
Definition 1 (Higher-degree monodromy set). Let M(λ) := {A, B, C, D}(λ) for λ ∈ C be a set whose elements are the operators defined in (2.8). Then define the monodromy set of degree n as W n := M(λ 1 ) ⊗ M(λ 2 ) ⊗ . . . M(λ n ) and W n : Definition 2 (Higher-degree reflection relation). Let SK 2 ⊂ W 2 be the set of quadratic relations originated from the reflection algebra as described in [23]. We then define the reflection relations of degree n as the elements of SK n ⊂ W n defined recursively through the relation In other words, the relations in SK n consist of the expressions obtained from W n (n > 2) after the repeated use of the relations in SK 2 1 .
The set SK 2 is formed by the fundamental reflection algebra relations and they amount to sixteen relations in the case considered here. We shall only need a few relations in SK 2 for our purposes. For instance, we have B(λ 1 )B(λ 2 ) = B(λ 2 )B(λ 1 ) and C(λ 1 )C(λ 2 ) = C(λ 2 )C(λ 1 ) among them. These two relations motivates the use of the following simplified notation.  We can safely use the notation (3.2) since the operators B's commute for different values of their spectral parameters. The same argument also applies for the operators C's.
As far as the evaluation of the scalar product S n is concerned, we shall make use of the following relations in SK n+1 : In (3.3) and (3.4) we have considered the sets X i,j and Y i,j respectively defined as X i,j : In addition to that we have also employed the sets X i,j λ := X i,j \{λ} and Y i,j λ := Y i,j \{λ}. The relations (3.3) and (3.4) do not exhaust all possible relations in SK n+1 which can be used to determine the scalar product S n . For instance, here we shall also make use of the following ones: (3.6)

Algebraic-functional approach
In this section we aim to derive functional relations determining the scalar product (2.15). The method we shall employ consists of an extension of the work [20] for the case of open boundary conditions. The key idea is to use the reflection algebra as a source of functional relations as discussed in [23]. This method was first proposed in [16] for spectral problems and extended for vertex models partition functions in the works [17-19, 47, 48]. We refer to this approach as Algebraic-Functional Method and an important step within it is to find a suitable linear map π n : SK n → C[λ ±1 1 , λ ±1 2 , . . . , λ ±1 n ] able to convert a higher-degree reflection algebra relation into a complex multivariate function. Moreover, we would like to build a map π n yielding the simplest possible functional equation for the scalar product (2.15). In fact, the structure of the resulting functional equation will depend on two ingredients: the higher-degree relation in SK n we single out and the particular realization of π n we are considering.
Proposition 1 (The map π n ). The following scalar products are realizations of π n , Proof. From (3.7) we can readily see that π B,C n : . Therefore, as SK n ⊂ W n , our claim is automatically fulfilled.
In what follows we shall demonstrate how one can use (3.7) to convert the relations (3.3)-(3.6) into two functional equations characterizing the scalar product S n .

Equation type A
The functional equation obtained from the reflection algebra relations (3.3) and (3.4) shall be refereed to as equation type A. This equation is a consequence of the trivial identity which is the starting point to prove Theorem 1. For convenience we also introduce the following notation.
Definition 4. Let i, j ∈ Z such that i < j. Then define X i,j and Y i,j as the following (j − i + 1)-dimensional vectors, Also, consider λ ∈ X i,j andλ ∈ Y i,j , and additionally define the (j − i)-dimensional vectors (3.10) In (3.10) the underline stands for omission. Theorem 1. The scalar product S n , as defined in (2.15), satisfies the equation with coefficients Proof. Within the algebraic-functional approach, we firstly consider the action of the map π C n+1 , defined in (3.7), on the relation (3.3). By doing so we are left with terms of the form: . Those terms satisfy the reduction property π C n+1 → π C n due to the highest-weight property of the vector |0 entering in the definition (3.7). Next we consider the action of π B n+1 on the higher-degree relation (3.4). This procedure yields terms of the form Similarly, those terms also enjoy the property π B n+1 → π B n due to the highest-weight property of (2.10). The precise relations associated with π B,C n+1 → π B,C n can be found with the help of (2.11). Lastly, we use the identifications Interestingly, Eq. (3.11) exhibits the same structure as the equations derived in [20] for the scalar products of Bethe vectors associated with the XXZ spin chain with periodic boundary conditions. The only difference is in the particular form of the coefficients M 0 , N (B) λ and N (C) λ . In [20] we have presented two functional equations describing the desired scalar product and the situation is similar in the present paper. A second functional equation satisfied by the scalar product (2.15) will be described in the next section.

Equation type D
As anticipated in the previous section, there exists a second functional equation satisfied by the scalar product S n . This second equation exhibits the same structure as (3.11). In fact, a large number of functional equations could be derived by exploring different reflection algebra relations and different realizations of the maps π B,C n . This feature of the algebraic-functional method has been previously observed for partition functions of vertex models with domain-wall boundaries as discussed in [48]. However, we would like to deal with the simplest possible equations and a systematic method for solving equations with the structure of (3.11) has been already developed in [19,47,20].
This second equation for the scalar product S n will be obtained from the higherdegree relations (3.5) and (3.6). The derivation will also require an equivalent of (3.8) and we shall use the following one, Then the functional equation is satisfied by the scalar product S n .
Proof. We follow the same strategy used in the proof of Theorem 1. The first step is to apply the map π C n+1 on the relation (3.5). By doing so we are left with terms , which obey a reduction relation π C n+1 → π C n due to the sl(2) highest-weight property of (2.10). Next we apply π B n+1 to (3.6). The terms obtained through this procedure are of the following form: The reduction property π B,C n+1 → π B,C n then allows one to identify the scalar product S n through (3.13). The mechanism above described then yields the functional relation (3.16).
Remark 1. At this stage it is important to remark some structural properties of (3.11) and (3.16). Although we can readily see that those equations are invariant under the permutations of variables λ B i ↔ λ B j and λ C i ↔ λ C j for any i, j ∈ {1, 2, . . . , n}, the same does not hold for the permutations λ 0 ↔ λ B i and λ 0 ↔ λ C i . Therefore, the latter permutations indeed produces significantly more independent equations describing the scalar product S n .

The scalar product S n
In the previous section we have derived two functional equations satisfied by the scalar product S n . This set of equations is described in Theorems 1 and 2, and their derivation makes use of higher-order reflection algebra relations and the highest-weight property of the vector |0 entering in the definition (2.15). However, it is important to stress here that equations (3.11) and (3.16) carry no information about the particular representation of the operators B's and C's we are considering in the definition of the scalar product (2.15). The choice of representation characterizes the particular model for which the scalar products are being computed. In this way, it is expected that different choices of representations for the operators B's and C's will be manifested in different classes of solutions of the system (3.11, 3.16). Thus, we need to characterize the class of solutions of (3.11, 3.16) we are interested. The class of solution corresponding to the scalar product of Bethe vectors associated with the XXZ chain with open boundaries (2.15) will be discussed in what follows, but some remarks are in order before proceeding with that.
For instance, the system of equations (3.11, 3.16) exhibits the same structure as the one derived in [20] for the scalar product of Bethe vectors associated with the XXZ chain with periodic boundary conditions. The main difference between our present equations (3.11, 3.16) and the equations of [20] lies on the particular form of the equations coefficients. The steps required for solving (3.11, 3.16) need to be modified due to that difference but the general idea remains the same. In fact, the general strategy for solving this type of functional relations was firstly described in [19,47]. The resolution of (3.11, 3.16) will follow a number of systematic steps exploiting properties expected for the desired scalar product (2.15).
Lemma 2 (Doubly symmetric function). Let Z := (λ 1 , λ 2 , . . . , λ n ) be a generic n-tuple. Then the scalar product S n satisfies the symmetry property, In addition to that we also have the property Proof. See Appendix A.
Remark 2. Due to Lemma 2 we can also consider the simplified notation

4)
whereS n is a doubly symmetric polynomial of order 2L in each one of its variables.
Proof. See Appendix B.
Lemma 4 (Special zeroes). The function S n (X 1,n |Y 1,n ) vanishes for the specialization of variables Similarly, it also vanishes for Proof. See Appendix D.
Proof. See Appendix C.
The above lemmas pave the way for the resolution of the system of equations (3.11, 3.16). The desired solution is given by the following theorem. Theorem 3. The scalar product S n can be written as the following multiple contour integral, where . Proof. The proof follows from the resolution of the system of equations (3.11, 3.16) under the conditions established by Lemmas 3 and 5. In fact, these two lemmas are the only ones that do not follow from the functional equations (3.11) and (3.16), and they characterize the class of solution we are interested. Our procedure will be described by a sequence of systematic steps.
Step 6. Next we set λ C n = µ 1 in (4.13) and compare the result with (4.12). This allows us to conclude that V = W . The same conclusion can be obtained by setting λ B n = µ 1 in (4.15) and comparing the result with (4.11). It is worth remarking that a similar analysis using the results of Steps 4 and 5 also yields the same constraint.
Step 7. We set λ 0 = µ 1 in Eq. (3.11) and use the relations (4.16) and (4.18). From this point on, we are already assuming V = W as obtained in Step 6. This procedure allows us to write S n (X 1,n |Y 1,n ) = λ∈X 1,n λ ∈Y 1,n K λλ V (X 1,n λ |Y 1,n λ ) . The function K λλ is given by where is a constant, i.e. independent of λ andλ, and Γ λλ is a 2 × 2 matrix with entries . (4.23) Step 8. Next we look at Eq. (3.16) under the specialization λ 0 = µ 1 − γ. This particular specialization allows us to readily use formulae (4.13) and (4.15). By doing so we obtain the relation S n (X 1,n |Y 1,n ) = The termC 0 is a constant and it explicitly reads HereΓ λλ is also a 2 × 2 matrix with entries given by . Step 1 Step 2 Step 3 Step 4 Step 5 Step 6 Step 6 Step 7 Step 8 Step 9 Step 10 Step 9. Substitute formula (4.19) in Eq. (3.11) and set λ B n = µ 1 − γ and λ C n = µ 1 . This yields an equation for the function V exhibiting the same structure of (3.11), but with modified coefficients and n → n − 1. The explicit form of the coefficients will not be relevant for our purposes here as one can verify that the resulting equation consists of a linear combination of (3.11) and (3.16) under the mappings L → L − 1, n → n − 1 and µ i → µ i+1 .
Step 10. We repeat the procedure of Step 9 but substitute instead formula (4.19) in Eq. (3.16). We then set λ C n = µ 1 − γ and λ B n = µ 1 to find another equation for the function V with the same structure of (3.16). Similarly to the previous step, one can verify that the resulting equation consists of a linear combination of (3.11) and (3.16) taking into account the mappings L → L − 1, n → n − 1 and µ i → µ i+1 .
The sequence of steps 1 through 10 is interrelated and the dependence of a single step with the remaining ones is schematically depicted in Figure 1. Some comments are also required at this stage. For instance, in step 7 we have obtained formula (4.19) expressing the scalar product S n as a linear combination of auxiliary functions V . The coefficients of this linear combination is given by the function K λλ defined in (4.20). Alternatively, in step 8 we have also obtained a similar formula expressing S n in terms of the same functions V . This formula is given by (4.24) and it is expressed in terms of coefficients K λλ defined in (4.25). Moreover, steps 9 and 10 tell us that the function V satisfy the same system of equations (3.11, 3.16) subjected to the mappings L → L − 1, n → n − 1 and µ i → µ i+1 . This conclusion has been obtained by inserting the expression (4.19) back into the system (3.11, 3.16). Alternatively, this same conclusion can be obtained by replacing (4.24) into (3.11, 3.16). In its turn, the function V is essentially a polynomial of order 2(L−1) as discussed in step 1. Polynomial solutions of this type of linear functional equations are unique as demonstrated in [19] and this implies that V is essentially the scalar product S n−1 on a lattice of length L − 1. Therefore, we can use (4.19) or (4.24) to obtain S n recursively starting from the solution of (3.11, 3.16) for n = 1. In particular, the iteration procedure described by (4.19) and (4.24) seems to be naturally realized by multiple contour integrals.
Multiple contour integrals. In order to find an explicit expression realizing the iteration procedure described by (4.19), we assume the following ansatz for the scalar product S n , where H is a function to be determined. Also, the integration contour for each variable w i in the formula (4.28) must only enclose all variables in the set X 1,n . Analogously, each variablew i is integrated along a contour containing only all variables in the set Y 1,n . It is important to remark here that formula (4.28) concentrates the dependence with all variables λ B j and λ C j in the denominator of its integrand. In this way, we are left with an integral representation for S n if we are able to exhibit a function H implementing the iteration described by (4.19). This approach has been already discussed in [20] and here we shall restrict ourselves to presenting only the main points of this procedure. The substitution of (4.28) in (4.19) then yields the following recursion relation for H, The functionH in (4.29) consists of H, up to an overall multiplicative constant factor, under the mappings L → L − 1, n → n − 1 and µ i → µ i+1 . The function R 1 reads while the 2 × 2 matrix Φ (1) has entries given by for variables s m i := δ 1m w i + δ 2mwi . The iteration of (4.29) starting with the solution for the case n = 1 obtained in the Appendix E yields formula (4.8). In fact, equations (3.11) and (3.16) are homogeneous and they can only determine the scalar product S n up to an overall multiplicative constant. This constant is then fixed by the asymptotic behavior (4.7). This completes our proof.
Remark 3. The relation (4.24) could also have been used instead of (4.19) to produce a multiple contour integral representation for the scalar product S n . However, one can notice thatK λλ corresponds essentially to K λλ under the map µ i → −µ i . This property is only violated by the constant factors C 0 andC 0 . Since the overall constant factor needs to be fixed by the asymptotic behavior (4.7), this issue will not be relevant. In this way, the use of (4.24) produces the same representation (4.8) under the map µ i → −µ i . This also implies that the transformation µ i → −µ i is a discrete symmetry of the scalar product S n .

Concluding remarks
This work was devoted to the analysis of scalar products of Bethe vectors associated with the XXZ spin chain with open boundary conditions. In particular, our study is based on the description of scalar products by means of functional equations.
The desired scalar products are shown to satisfy a system of functional relations originated from the reflection algebra. This approach for scalar products of Bethe vectors was firstly proposed in [20] using the Yang-Baxter algebra as the source of functional equations. The feasibility of using the reflection algebra in a similar way has been put forward recently in [23]. Although the quantity computed in the present paper is different from the one considered in [23], the corresponding functional equations are derived from the same reflection algebra relation. More precisely, the higher-degree relation (3.3) is the same one employed in [23]. Here, however, we need to use a total of four reflection algebra relations, namely (3.3), (3.4), (3.5) and (3.6), and the main difference compared to [23] is the choice of map π n . In the present work, two different realizations π B,C n was required in order to describe the desired scalar product.
The methodology used to solve this type of functional equations has been introduced in [18] for an equation describing the partition function of a Solid-on-Solid model with domain wall boundaries. The key point of this method is that the location of special zeroes of the partition function or scalar products induces a separation of variables. This separation of variables can also be regarded as a recurrence relation allowing us to readily obtain the desired solution. Here the solution of our system of equations is obtained in a compact form. It is given by a multiple contour integral as stated in Theorem 3. Interestingly, this type of integral formula seems to be a general structure accommodating quantities such as scalar products and partition functions associated with integrable vertex models. For instance, similar integral formulae have appeared in [37,19,47,23] for partition functions with domain-wall boundaries and in [37,20] for scalar products.
Functional equations such as (3.11) and (3.16) also to encode a family of partial differential equations describing certain quantities of interest. This has been shown in [18,48,13,23] by recasting our functional relations in an operatorial form. We have not pursued that direction in the present paper but we hope to report on that in a future publication.

Acknowledgements
The author thanks G. Arutyunov for discussions and comments on this manuscript. The work of W.G. was supported by the German Science Foundation (DFG) under the Collaborative Research Center (SFB) 676 Particles, Strings and the Early Universe.

A S n as a doubly symmetric function
This appendix is concerned with the proof of Lemma 2. This lemma states that the function S n (λ C 1 , . . . , λ C n |λ B 1 , . . . , λ B n ) is invariant under the permutation of variables λ C i ↔ λ C j and λ B k ↔ λ B l independently. We thus say that S n is a doubly symmetric function. The definition of S n , as given in (2.15), and the commutation rules [B(λ 1 ), B(λ 2 )] = [C(λ 1 ), C(λ 2 )] = 0 encoded in the reflection algebra, already requires this property to be fulfilled. However, once we assume the scalar product to satisfy equations (3.11) and (3.16), we would like their solution to naturally inherit this symmetry property. That is precisely what we intend to show here and our proof will follow the same arguments employed in [20,23].
Firstly, we assume S n to be an analytic function and inspect (3.11) in the limit λ 0 → λ C k . From (3.12) we can see that M 0 and N are the only singular coefficients in that limit. Thus the integration of (3.11) with respect to the variable λ 0 around a contour enclosing solely the variable λ C k yields the following identity, Moreover, from formulae (3.12) we can readily verify the property which allows us to conclude that 2) Next we integrate Eq. (3.11) over the variable λ 0 around a contour containing solely the variable λ B k . Similarly to (A.2), this procedure allows us to conclude that 3) The relations (A.2) and (A.3) tell us that S n is invariant under cyclic permutations of λ C,B 1 , . . . , λ C,B k for any k in the interval 1 ≤ k ≤ n. Then, since the symmetric group of order n is generated by any cycle of length n, in addition to any single transposition, we can conclude that S n is invariant under the action of Sym (X 1,n ) ⊗ Sym (Y 1,n ). This concludes our proof. Although here we have considered only Eq. (3.11), it is worth remarking that the same results could have been obtained from a similar analysis of Eq. (3.16).
As we can see from (2.15), the dependence of S n with the set of variables X 1,n and Y 1,n is only due to the operators C and B respectively. More precisely, the whole dependence with the variable λ C i comes from the operator C(λ C i ) entering the definition (2.15), while the dependence with λ B i is characterized by the operator B(λ B i ). Therefore, Lemma 3 follows from the functional dependence of these operators with their variables. Using the notation x := e 2λ , it is then sufficient to show that The proof of (B.1) for the operator B can be found in [23] and we shall not repeat it here. On the other hand, the functional dependence (B.1) for the operator C is a direct consequence of the proof presented in [23] and the fact that C(λ) corresponds to B(λ) t under the map µ i → −µ i as demonstrated in [16]. In what follows we shall describe the main points leading to the desired property. Definition 6. Let A, B, C, D,Ā,B,C,D : C → End (C 2 ) ⊗L be operator-valued functions defined as Taking into account (2.8) and (B.2) we can then write Moreover, in [16] we have demonstrated that the operators defined in (B.2) satisfy the following properties: as a consequence of the crossing symmetry exhibited by the U q [ sl (2)] invariant R-matrix employed in the present work. Using (B.3) and (B.4) we can readily show that the relation C(λ) = B(λ) t | µ i →−µ i holds, which concludes our proof.

C Asymptotic behavior
The full determination of the scalar product S n as solution of the system of functional equations (3.11, 3.16) requires we are able to evaluate S n for a particular value of its variables. This is due to the fact that (3.11) and (3.16) are homogeneous equations and, therefore, they are only able to determine the solution up to an overall multiplicative factor independent of the variables λ B,C i . Here we find that the points x B,C i := e 2λ B,C i → ∞ are suitable for that purpose, and the first step to compute S n in that limit is the analysis of the behavior of B(x) and C(x) for x → ∞.
In [23] we have shown that the operator B(x) exhibits the asymptotic behavior Here we are using the conventions q := e γ , t := e h , y i := e 2µ i and, due to the property C(λ) = B(λ) t | µ i →−µ i shown in Appendix B, we readily find that in the limit x → ∞. The operatorsP ± j appearing in (C.3) are in their turn defined as The symbol id in (C.2) and (C.4) stands for the identity matrix in End(C 2 ), while X ± and K denote the generators of the U q [sl (2)] algebra in the fundamental representation as described in [23]. As a consequence of the U q [sl (2)] algebra satisfied by X ± and K, one can show that the operators P ± j are subjected to the following commutation rules: The operatorsP ± j obey similar relations, namelȳ For latter convenience we also define the operators in such a way that Q as a consequence of (C.5) and (C.6), Now the expressions (C.1), (C.3) and (C.7) allow us to writē in the limit x B,C i → ∞ for all i ∈ {1, 2, . . . , n}. Formula (C.9) is given in terms of the operators J andJ defined respectively as rs . (C.11) Next we use the methodology described in [23] to find that J andJ can be rewritten as where ∆ ± l := l m=0 q ±2m . In order to find an explicit expression for the leading coefficient (C.9) we still need to evaluate the quantity 0|JJ |0 . This task can be readily performed with the help of (C.2), (C.4), (C.12) and (C.13). By doing so we find which leads us directly to the statement of Lemma 5.

D Special zeroes
The resolution of the system of Eqs. (3.11, 3.16) follows a sequence of systematic steps relying heavily on Lemma 4. This lemma gives us the location of special zeroes of the scalar product S n inducing the separation of variables expressed by the relations (4.19) and (4.24). More precisely, we have considered three pairs of points, namely (λ C 1 , λ C 2 ) = (µ 1 − γ, µ 1 ), (λ C 1 , λ C 2 ) = (µ 1 − γ, −µ 1 − γ) and (λ C 1 , λ C 2 ) = (−µ 1 , µ 1 ), for which the specialization of S n vanishes. The same property with variables λ C i exchanged by λ B i also holds. The location of these special points is encoded in our system of equations and this is what we intend to unveil here. Also, we shall only present a detailed proof of this property for (λ C 1 , λ C 2 ) = (µ 1 − γ, µ 1 ) as the remaining cases can be obtained along the same lines. For illustrative purposes we shall firstly discuss the cases n = 2, 3, and only then describe the general case.
Case n = 2. We consider (3.11) under the specializations λ C 1 = µ 1 − γ and λ C 2 = µ 1 . Under this particular specialization one can verify that the coefficients N vanish. We are then left with the equation where | 1,2 denotes the prescribed specialization of λ C 1 and λ C 2 . Next we consider Remark 1 and apply the maps λ 0 ↔ λ B 1 and λ 0 ↔ λ B 2 taking into account Lemma 2. This produces two additional equations with the same structure of (D.1) but modified coefficients. We then define the coefficients in such a way that the resulting equations can be conveniently written as Now one can verify that det N i j = 0 for generic values of the variables and this allows us to conclude that S n (µ 1 − γ, µ 1 |λ B 1 , λ B 2 ) vanishes. It is worth remarking here that we have used only Eq. (3.11) for the case n = 2. For n > 2 we shall also need (3.16).
Case n = 3. In that case we look at (3.11) and (3.16) under the specializations vanish for this particular specialization and we are left with the following equations: Similarly to the case n = 2, in (D.4) and (D.5) we have used the symbol | 2,3 to denote the aforementioned specializations of λ C 2 and λ C 3 . We then eliminate the term S 3 (λ 0 , µ 1 − γ, µ 1 |λ B 1 , λ B 2 , λ B 3 ) from the system of equations formed by (D.4) and (D.5). In addition to that we also consider the maps λ 0 ↔ λ B i for 1 ≤ i ≤ 3 to produce three extra equations. We are thus left with a total of four linear equations which can be conveniently written as (D. 6) The coefficients N i j in (D.6) are then defined as and along the lines used for n = 2, one can also verify here that det N i j = 0 for generic values of the variables. In this way we can conclude that General case. We consider n levels of specializations of Eqs. (3.11) and (3.16), and at each k-level we set λ C n−k = µ 1 − kγ for 0 ≤ k ≤ n − 1, keeping the specializations at the previous levels. At the final level k = n − 1 we are left with the following equations, S n (X * * |Y 1,n ) + Next we eliminate the term S n (X * * |Y 1,n ) from the system of equations (D.8) and consider the n additional equations obtained from the map λ 0 ↔ λ B i for 1 ≤ i ≤ n. The system of equations obtained through this procedure can then be written as Now one can verify in (D.10) that det N i j = 0 for arbitrary values of the variables, and we can conclude that S n (X * |Y 1,n ) = 0. This is still not the relation we want to prove. For that we revisit each level of specialization backwards, taking into account the result obtained at the final level, and perform a similar analysis. By doing so we finally obtain the desired property S n (µ 1 − γ, µ 1 , λ C 3 , . . . , λ C n |λ B 1 , . . . , λ B n ) = 0.

E Solution for n = 1
The relations (4.19) and (4.24) consist of a separation of variables induced by the location of certain zeroes of the scalar product S n . Those zeroes are given in Lemma 4 and, alternatively, formulae (4.19) and (4.24) can also be regarded as recurrence relations. This is due to the fact that the function V corresponds to the function S n−1 under the map L → L − 1 and µ i → µ i+1 as shown in Section 4. Here we shall consider that n ≤ L in such a way that the last step of the iteration procedure described by (4.19) and (4.24) will require the solution of (3.11, 3.16) for n = 1 and arbitrary L. In that case the solution can be obtained through simple algebraic manipulations. The sequence of steps required to obtain the desired solution will be described in what follows. The equations (3.11) and (3.16) for n = 1 simply read with coefficients explicitly given by and We then proceed by eliminating the term S 1 (λ 0 |λ B 1 ) from (E.1). By doing so we find the separated relation, with function ψ(λ,λ) reading (E.5) The variable λ C 1 plays the role of a parameter in (E.4) and we can readily conclude that , (E.6) where F : C → C is a function yet to be determined. Next we solve (E.1) for S 1 (λ C 1 |λ B 1 ) instead and use formula (E.6). This leaves us with the following equation for the function F , Eq. (E.7) is readily solved by F (λ) = κ b(2λ) a(2λ) where κ is a constant. Gathering our results we then find the solution . (E.8) In (E.8) we have already considered κ = c in accordance with the asymptotic behavior (4.7).