Reflection algebra and functional equations

In this work we investigate the possibility of using the reflection algebra as a source of functional equations. More precisely, we obtain functional relations determining the partition function of the six-vertex model with domain-wall boundary conditions and one reflecting end. The model's partition function is expressed as a multiple-contour integral that allows the homogeneous limit to be obtained straightforwardly. Our functional equations are also shown to give rise to a consistent set of partial differential equations for the partition function.


D Special zeroes 22 E Asymptotic behavior 24 1 Introduction
The study of correlation functions of quantum-integrable systems is intrinsically related to partition functions of vertex models with special boundary conditions. In particular, the case of domain-wall boundaries is of fundamental importance and this fact was first noticed for the one-dimensional Heisenberg chain with periodic boundary conditions and the associated six-vertex model [1]. Among the results of [1] there is a recurrence relation determining the partition function of the six-vertex model with domain-wall boundary conditions, which was later on solved by Izergin in terms of a determinant [2]. Subsequently, a determinant representation for scalar products of Bethe vectors under certain specializations of the parameters (so called on-shell scalar products) was obtained by Slavnov [3]. Those results represented the first steps towards the computation of exact correlation functions of quantum-integrable systems, but it is worth remarking that the problem of computing norms of Bethe wave functions was first considered by Gaudin in the case of the non-linear Schrödinger model [4]. Gaudin also formulated the hypothesis that such norms would be given by certain Jacobian determinants. This hypothesis was subsequently proved by Korepin [5] even for more general models. Scalar products of Bethe vectors are the building blocks of correlation functions, and having them expressed as determinants for the Heisenberg chain paved the way for the use of well-established techniques. For instance, the aforementioned determinant formulae in the thermodynamic limit become determinants of integral operators, i.e. Fredholm determinants, allowing the derivation of differential equations describing correlation functions [6]. The asymptotic behavior of correlation functions was then derived from the analysis of the respective differential equations. Interestingly, the language of differential equations seems to be quite appropriate for the description of correlation functions. A remarkable example of this is provided by the Knizhnik-Zamolodchikov equation describing multi-point correlation functions of primary fields in conformal field theories [7].
As far as the one-dimensional Bose gas is concerned, the results of [6] are largely due to the existence of particular determinant representations for correlation functions. However, the recent works [8][9][10][11] have raised doubts on the existence of such determinant representations for models based on higher-rank symmetry algebras. Nevertheless, it is still important to remark here that several advances in the computation of correlation functions have been obtained through the algebraic Bethe ansatz in [12][13][14][15][16][17] and through Sklyanin's separation of variables in [18,19].
Integrable spin chains can also be addressed directly in the thermodynamic limit through the vertex-operator approach due to the quantum affine symmetry exhibited in that limit [20,21]. Alternatively, one can also consider the q-Onsager approach described in [22]. Within the vertex-operator approach the description of correlation functions is done by means of the quantized Knizhnik-Zamolodchikov equation [23][24][25]. The latter is not a differential equation but rather a functional equation describing the matrix coefficients of a product of intertwining operators for a given affine Lie algebra [26]. Interestingly, the partition function of the six-vertex model with domain-wall boundaries, which initiated the study of correlation functions for one-dimensional quantum spin chains, can also be described by a functional equation resembling certain aspects of the Knizhnik-Zamolodchikov equation [27]. Moreover, the functional equation for that partition function can further be translated into a partial differential equation, as was shown in [28,27].
The derivation of such functional equations is based on an algebraic-functional approach initially proposed for spectral problems in [29] and extended for correlations in the series of works [30,28,[31][32][33]. In particular, those equations allow for the derivation of integral representations for partition functions with domain-wall boundaries [31,32] and scalar products of Bethe vectors [33]. A fundamental ingredient within this algebraicfunctional approach is the Yang-Baxter algebra, which is a common algebraic structure underlying quantum-integrable systems. On the other hand, integrable systems with open boundary conditions are governed by the reflection algebra [34] and in the present paper we show that this algebra can also be exploited along the lines of [27].
Correlation functions for the Heisenberg chain with open boundary conditions have been studied in the literature through a variety of approaches [35][36][37]. In particular, a partition function with domain-wall boundaries has also been defined in that case [38]. As a matter of fact, the partition function introduced in [38] considers both domain-wall boundaries and one reflecting end, and it can also be expressed as a determinant along the lines of [2]. Moreover, this partition function has also found interesting applications in the computation of physical properties of the XXZ spin chain at finite temperature. For instance, the surface free energy of the XXZ spin chain has been expressed in [39] as the expectation value of a product of projection operators. This expectation value was then demonstrated in [40] to be precisely the partition function of the six-vertex model with one reflecting end and domain-wall boundaries described in [38]. Here we shall study this partition function through the algebraic-functional method developed in [31,32]. This approach will not only allow us to find a new representation for the model's partition function, but it will also unveil a set of partial differential equations satisfied by the latter.
This paper is organized as follows. In Section 2 we introduce definitions and conventions which shall be employed throughout this work. In Section 3 we describe the algebraic-functional approach in terms of the reflection algebra, and use this method to derive functional equations characterizing the partition function of the six-vertex model with one reflecting end and domain-wall boundaries. The solution of our equation is then given in Section 3. We proceed with the analysis of our functional equation in Section 3.3 where we extract a set of partial differential equations satisfied by the model's partition function. Concluding remarks are discussed in Section 4 while proofs and technical details are left for the Appendices.

Definitions and conventions
The most well studied cases of integrable lattice systems are those defined on a finite interval with periodic boundary conditions, although more general types of boundaries can also be considered. Among those systems, a prominent class is formed by onedimensional spin chains with open ends and particular terms characterizing the reflection at the boundaries. Those models can also be solved by means of the Bethe ansatz and the first results on that direction have been obtained in [41,42]. Boundary terms are normally not expected to modify the infinite-volume properties of physical systems, although counterexamples for that common belief have been reported in the literature [43]. Nevertheless, even for the cases where infinite-volume properties remain the same, boundary terms can still change the finite-size corrections of massless systems defined on a strip of width L [42]. The latter is able to provide fundamental information concerning the underlying conformal field theory as shown in [44].
The study of lattice systems with open boundary conditions gained a large impulse with the formulation of the Quantum Inverse Scattering Method (QISM) for that class of models [34]. The approach developed in [34] is based on Cherednik's condition for factorised scattering with reflection [45], although it can also be regarded in the context of vertex models of Statistical Mechanics [46]. The latter is the perspective to be adopted here, and in what follows we shall briefly describe the ingredients required for the construction of the six-vertex model with domain walls and one reflecting end as defined in [38].
The U q [ sl(2)] R-matrix. Integrable vertex models are characterized by an R-matrix R : C → End(V ⊗ V) satisfying the Yang-Baxter equation, namely In Eq. (2.1) we are using the standard notation R ij ∈ End(V i ⊗ V j ) and λ i denote arbitrary complex parameters. For the six-vertex model we have V i := V ∼ = C 2 and where a(λ) = sinh (λ + γ), b(λ) = sinh (λ) and c(λ) = sinh (γ) with anisotropy parameter γ ∈ C.
Monodromy matrices. Let V 0 ∼ = C 2 and V Q := (C 2 ) ⊗L for lattice length L ∈ Z >0 . Also let λ and µ j (1 ≤ j ≤ L) be arbitrary complex parameters. Then we consider operators τ,τ : C → End(V 0 ⊗ V Q ) defined as the following ordered products: The operators τ andτ are usually referred to as monodromy matrices while V 0 and V Q are called auxiliary and quantum spaces respectively. Since V 0 ∼ = C 2 , the monodromy matrices (2.3) can be recasted as with entries in End(V Q ).
Yang-Baxter algebra. Let R be given by (2.2) and consider U ∈ {τ,τ }. Then the following quadratic algebra is fulfilled by U, due to the Yang-Baxter equation (2.1). Here we are employing the notation Y 1 := Y⊗id 2 and Y 2 := id 1 ⊗ Y for any operator Y ∈ End(V), where the symbol id j stands for the identity operator in End(V j ). The relation (2.5) is commonly referred to as Yang-Baxter algebra and it describes commutation relations for the entries of the monodromy matrices (2.4).

Reflection equation.
Within the framework of the QISM developed in [34], integrable boundary conditions are characterized by a matrix K satisfying the so called reflection equation. This equation, which was first proposed in [45] in the context of factorised scattering, reads Here we shall restrict ourselves to a particular solution of (2.6) associated to the R-matrix (2.2). This solution is explicitly given by where κ ± (λ) = sinh (h ± λ) and h ∈ C is an arbitrary parameter describing the interaction at one of the boundaries.

Remark 1.
Within the context of one-dimensional integrable spin chains, the K-matrix (2.7) describes the reflection only at one of the ends of an open system. If we would want to describe reflection at the opposite end we would also need to introduce a matrixK satisfying an independent equation isomorphic to (2.6). More details on this construction can be found in [34].
Reflection algebra. Let the operator T : Then, following [34], one can show that T satisfies the following quadratic algebra, The relation (2.9) will be referred to as reflection algebra and it follows from the reflection equation obeyed by the matrix K together with properties satisfied by τ andτ . The operator T is commonly referred to as double-row monodromy matrix and, similarly to (2.4), it can be recasted as In this way (2.9) encodes commutation relations for the operators A, B, C, D ∈ End(V Q ).
Highest/lowest weight vectors. The vectors |0 , |0 ∈ V Q defined as Partition function. Following the work [38], the partition function of the six-vertex model with one reflecting end and domain-wall boundaries is given by This partition function is a multivariate function depending on L spectral parameters λ j , L inhomogeneity parameters µ j , the anisotropy parameter γ and the boundary parameter h. In Section 3 we shall describe how the reflection algebra (2.9) can be exploited in order to derive functional equations determining the partition function (2.13).
Diagrammatic representation. The lattice system described by the partition function (2.13) can be more intuitively depicted in terms of diagrams representing the action of the Yang-Baxter and reflection algebras elements. For that it is convenient to write R = R α β αβ e αα ⊗ e ββ , K = K α α e αα and T = T α α e αα , where summation over repeated indices is assumed. Here α, α , β, β ∈ {1, 2} label the basis vectors of V ∼ = C 2 , while e αβ is the matrix with entries (e αβ ) ij = δ αi δ βj . The diagrammatic representation of R and K is given in Figure 1 while T is depicted in Figure 2. Using these conventions the partition function (2.13) is illustrated in Figure 3 with external indices assuming the domain-wall configurations α j , β j = 1 and α j , β j = 2 for all j.

Algebraic-functional approach
In the works [29,30,28,[31][32][33]27] we have described a mechanism yielding functional equations satisfied by quantities of physical interest as a direct consequence of the Yang-Baxter algebra. This approach has been employed for the determination of spectra [29,32] and partition functions [31,32] of integrable vertex models. One issue arising within this method is that the algebraic relations we are considering might not suffice to determine Figure 1: Diagrammatic representation of the R-and K-matrices. the desired quantities. Furthermore, it would be desirable to have the simplest possible equations such that finding its solutions can be achieved without much effort. Up to the present moment, we have only considered the Yang-Baxter algebra and its dynamical counterpart as a source of functional relations [27] and here we aim to show that the reflection algebra (2.9) can also be exploited along the same lines. For this we need to introduce the following definitions.
be the ring of meromorphic functions in the variables λ 1 , . . . , λ n and defineW n : To obtain functional relations from the reflection algebra we also need to introduce an appropriate linear map A suitable realization of (3.1) will be given shortly.
Reflection relation of degree n. The reflection algebra (2.9) encodes a set of sixteen commutation relations governing the elements of (2.10). It is clear from (2.9) that those commutation rules are quadratic and here they are referred to as reflection relations of degree two. The repeated use of (2.9) then yields relations inW n which shall be referred to as reflection relations of degree n.

Functional equations
In Definition 1 we have introduced a map π n assigning multivariate complex functions to the elements of the set W n . Here our goal is to evaluate the partition function (2.13), and this can be achieved from the study of suitable functional equations derived through the application of the map (3.1) on reflection relations of higher degree. This procedure will require the following ingredients: a suitable realization of the map π n and a convenient reflection-algebra relation. As a matter of fact, different functional relations can be derived for the partition function Z by changing these ingredients.
Realization of π n . The operatorial formulation of the partition function Z as given by (2.13) suggests that a suitable realization of π n is given by the following scalar product: for F ∈W n and vectors |0 , |0 ∈ V Q defined in (2.11).
Reflection-algebra relation. Next we look for appropriate reflection relations of higher degree from which we can find functional relations satisfied by the partition function Z. In order to build such higher-degree relations we start from the following fundamental commutation rules contained in (2.9): Note that the above relation is given in terms of the operatorD(λ) = D(λ) − c(2λ) a(2λ) A(λ). Next we describe a suitable functional relation satisfied by (2.13). Although the partition function Z is a multivariate function depending on L spectral parameters λ i , in addition to parameters µ i , h and γ, here we shall obtain a functional equation determining Z where only λ i play the role of variables. Theorem 1. The partition function of the six-vertex model with one reflecting end and domain-wall boundaries obeys the functional equation with coefficients M 0 and M i given by The functions Λ A , ΛD andΛ A were defined in (2.12).
Proof. Consider the following element of W n+1 , under the light of the reflection algebra (2.9). The repeated use of (3.3) yields the following reflection relation of order n + 1, (3.7) Next we set n = L and apply the map π L+1 given by (3.2) to (3.7). The left-hand side of (3.7) then yields the term ). Note that π L+1 reduces to π L due to the sl(2) highest/lowest weight properties exhibited by the realization (3.2). More precisely we have: Now we can identify the partition function Z(λ 1 , . . . , λ L ) = π L ( − → 1≤j≤L B(λ j )) on the right-hand side of (3.8). Thus the relations (3.7) and (3.8) under the above mentioned conditions result in the functional equation (3.4). This proves Theorem 1.

Remark 2. The functional equation (3.4) is invariant under the permutation of variables
. . , L}. This conclusion follows directly from Lemma 2 which will be stated below. However, the permutation λ 0 ↔ λ j yields a different functional equation for Z. The resulting equation exhibits the same structure as (3.4), with modified coefficients though. In this way (3.4) actually encodes a set of L + 1 equations.

The partition function Z
This section is devoted to the determination of the partition function (2.13) as a particular solution of the functional equation (3.4). A priori we do not have any guarantee that (3.4) is enough for that but direct inspection reveals that this is indeed the case for small values of the lattice length L. The general strategy for solving (3.4) will follow the same steps described in [32]. This is anticipated since the structure of (3.4) resembles that of the functional equation derived in [32] for the partition function of the elliptic SOS model with domain-wall boundaries. However, here we shall need to exploit some further properties of (3.4) which were not required in [32]. In order to clarify our methodology let us first stress some characteristics of our functional equation. Firstly, Eq. (3.4) is an equation for a complex multivariate function Z formed by a linear combination of terms containing Z(λ 1 , λ 2 , . . . , λ L ) with one of the variables λ i replaced by the variable λ 0 . Thus (3.4) runs over the set of variables {λ 0 , λ 1 , . . . , λ L }. In addition to that, our equation is homogeneous in the sense that if Z is a solution then ωZ also solves (3.4) for any ω ∈ C independent of the variables λ i . This property anticipates that we shall need to evaluate the partition function (2.13) for a particular value of its variables in order to having the desired solution completely fixed. Moreover, due to the linearity of Eq. (3.4), we need to address the question of uniqueness of the solution. The partition function (2.13) consists of a particular polynomial solution and the uniqueness within such class of solutions was proved in [31] under very general conditions.
Considering the above discussion the following lemmas will assist us through the determination of the partition function Z.
Proof. The proof is obtained by induction and can be found in Appendix B. Proof. This property follows from the structure of poles appearing in (3.5). See Appendix C for details.
Lemma 3 (Special zeroes). For L ≥ 2 the partition function Z vanishes for the specialization of variables λ 1 = µ 1 − γ and λ 2 = µ 1 . The same holds for the specialization Proof. The proof follows from the inspection of (3.4) under these specializations of variables, taking into account Remark 2. See Appendix D for details.
Lemma 4 (Asymptotic behavior). In the limit where all variables x i → ∞, the function Z behaves asZ
Proof. As x i → ∞ the generators (2.10) tend to the generators of the U q [sl (2)] algebra. The properties of the latter can be employed to demonstrate (3.9) as is shown in Appendix E.

Multiple integral representation
The resolution of (3.4) will follow a sequence of systematic steps based on Lemmas 1 to 4. The general procedure consists in finding suitable specializations of the variables λ 0 and λ L allowing us to invoke the above lemmas. The desired solution of (3.4) is given by the following theorem.
Theorem 2. The partition function of the six-vertex model with one reflecting end and domain-wall boundaries (2.13) can be written as (3.10) Proof. The proof of Theorem 2 follows from the resolution of (3.4) taking into account certain properties of (2.13). The procedure consists of three steps.
Step 1. We first set λ 0 = µ 1 − γ in Eq. (3.4). Under this specialization the coefficient M 0 is reduced to a single product. This specialization also produces terms of the form Z(X 2,L ) whereX i,j := {µ 1 − γ} ∪ X i,j . Now due to Lemmas 1 to 3 we can write where V is a polynomial of degree 2(L−1) in each variable x i up to an overall exponential factor. Thus this particular specialization yields the expression with coefficients κ and m i given by (3.13) Step 2. We substitute formula (3.12) back into the original equation (3.4). By doing so we are left with an equation involving only functions V. Next we set λ L = µ 1 in the resulting equation which then further simplifies tõ (3.14) The explicit form of the coefficientsM 0 andM i is not enlightening but it is worth remarking that for L = 2 we find that (3.14) corresponds to (3.4) with L = 1 and µ 1 replaced by µ 2 . This fact suggests that (3.14) should coincide with (3.4) after replacing L by L − 1 and µ i by µ i+1 . Unfortunately this is not the case for general values of L and we actually find that (3.14) consists of a linear combination of (3.4) along the lines of Remark 2. Nevertheless, this still ensures that V is essentially our partition function under the maps L → L − 1 and µ i → µ i+1 since polynomial solutions are unique.
Step 3. The results of Step 2 allows us to obtain an explicit representation for our partition function from the relation (3.12) in a recursive manner. In fact, formula (3.12) suggests the following ansatz for Z Z(X 1,L ) = . . .
where H is a function yet to be determined. In particular, here we also assume that the integration contours in (3.15) enclose all the poles at w i = λ j and that H contains no poles inside those integration contours. Then we consider the mechanism described in [33] to find the following relation determining the function H, (3.16) The functionH in (3.16) corresponds to H under the maps L → L − 1, µ i → µ i+1 up to an overall constant factor. In this way the relation (3.16) can be iterated once we know the function H(w 1 ). This function can be directly read from the solution of (3.4) for L = 1 which can be found in Appendix F. Thus the iteration of (3.16) yields the following expression for the function H, where Θ i is given by (3.10). Formula (3.17) already takes into account the asymptotic behavior stated in Lemma 4 and this completes the proof of Theorem 2.

Partial differential equations
In Section 3.1 we have derived a functional equation governing the partition function (2.13) as a direct consequence of the reflection algebra (2.9) and the highest/lowest weight property of the vectors |0 and |0 . Some properties of our functional equation have already been discussed in Section 3.2 and here we intend to demonstrate some further properties. More precisely, in this section we shall unveil a set of linear partial differential equations underlying (3.4). This type of hidden structure was first presented in [28] for a similar type of equation and subsequently developed in [27,47]. The first step towards that description is to recast (3.4) in an operatorial form. This can be achieved with the help of the operator D α i defined as follows. Definition 2. Let n ∈ Z >0 and α ∈ Z\{1, 2, . . . , n}. As before we write C[z ±1 1 , . . . , z ±1 n ] for the space of meromorphic functions on C n . Now consider the following operator . . , z n ) := f (z 1 , . . . , z α , . . . , z n ) .
In (3.19) we have made the dependence of L on λ 0 explicit to stress that this reformulation concentrates the whole dependence of our functional equation on λ 0 in the operator L. In particular, this property will allow us to extract a set of partial differential equations from (3.4) due to the fact that there exists a differential realization of (3.18) when we restrict the action of the operator D α i to a particular function space. In order to describe this differential realization we first need to introduce some extra definitions and conventions. Definition 3. Let K[z 1 , . . . , z n ] denote the multivariate polynomial ring in the variables z 1 , . . . , z n with coefficients in an arbitrary field K. We will also use the abbreviation K[z] := K[z 1 , . . . , z n ]. Using this shorthand notation, we define K m [z] ⊆ K[z] to be the subspace of K[z] formed by polynomials of degree m in each variable z i . Lemma 5. The differential operator is a realization of (3.18) on the space K m [z].
Proof. The proof follows from the series expansion of functions in K m [z]. The details of this analysis can be found in [28,27].
The realization (3.20) can not be directly substituted in (3.19) since the function Z we are interested in does not belong to K m [z]. However, as far as the functionZ defined in Lemma 1 is concerned, we have thatZ( 1 , . . . , y ±1 L , q ±1 , t ±1 ] and thus (3.20) can be employed. Here we use the notation of Lemma 4 where, in particular, x j = e 2λ j . We then define the rescaled coefficients and X i,j = {x k : i ≤ k ≤ j} as in Remark 3. Now we can substitute (3.20) in (3.22), and the next step of our analysis is to look at the analytical properties ofL(x 0 ) as function of x 0 or equivalently λ 0 . The explicit expressions for the coefficientsM 0 andM i are obtained from (3.5) and we can readily see thatL contains simple poles at the zeroes of a(2λ 0 ), b(λ 0 − λ i ) and a(λ 0 + λ i ). The residues ofL at the poles a(2λ 0 ) = 0 and b(λ 0 − λ i ) = 0 vanish but the same is not true for the poles at a(λ 0 + λ i ) = 0. Thus (3.22) is of the form whereL R (x 0 ) has no poles for x 0 ∈ C\{0}. Moreover, the direct inspection ofL R reveals that it is indeed a polynomial of the form with differential-operator valued coefficients Ω k . Now sinceL R (x 0 ) is a polynomial, the equationL R (x 0 )Z(X 1,L ) = 0 must be satisfied by each power of x 0 separately. In this way we are left with a total of 2L + 1 partial differential equations formally reading The operator Ω 2L . Due to (3.20) and the fact thatZ(X 1,L ) ∈ K 2L [x 1 , . . . , x L ], the differential operators Ω k are linear and contain partial derivatives with respect to the variables x i of order ranging from 1 to 2L. Although the explicit form of the operators Ω k for a given value of L can be computed from (3.22), (3.23) and (3.24), they mostly lead to cumbersome expressions which are not very enlightening. Fortunately, the situation for the leading term operator Ω 2L is more interesting and we find the following compact expression, The functions U and Y i in (3.26) explicitly read, whereā ω (x, y) := xω − y −1 ω −1 . Some comments are appropriate at this stage. To start with, the direct inspection of (3.26) for small values of the lattice length L reveals that our partial differential equation is fully able to determine the desired polynomial solution up to an overall constant factor that is fixed by Lemma 4. Moreover, the structure of (3.26) resembles that of a quantum many-body hamiltonian with higher derivatives and we can regard the partition function Z as the null-eigenvalue wave-function associated to Ω 2L . It is worth remarking here that a similar structure appeared previously for the standard six-vertex model with domainwall boundaries in [27]. In particular, the structure of Ω 2L is also shared by higher conserved quantities of the six-vertex model as demonstrated in [47]. To conclude we remark that although (3.26) results in a differential equation of order 2L, it can still be recasted as a system of first-order equations using the reduction of order procedure. This analysis is explicitly performed in Appendix G.

Concluding remarks
This work is mainly concerned with the interplay between functional equations and the reflection algebra in the framework developed in [30,33,27]. More precisely, here we have investigated the partition function of the six-vertex model with one reflecting end and domain-wall boundaries through this algebraic-functional approach. This methodology has been previously considered for the dynamical counterpart of the Yang-Baxter algebra in [28,31], and here we demonstrate the feasibility of the reflection algebra for that approach. From this analysis we obtain functional relations satisfied by the partition function of the six-vertex model with both domain-wall and reflecting boundaries. Interestingly, the equation presented here exhibits the same structure as the one obtained in [31,27] for a partition function with simpler boundary conditions. Although [31,27] and our present work consider domain-wall boundary conditions, here we have also included a reflecting end, which makes this algebraic-functional analysis significantly more involved. However, the difference between the functional equations in those works and the present one is restricted to the explicit form of their coefficients.
The starting point for the derivation of (3.4) is the element (3.6) and the corresponding reflection relation of higher degree (3.7). This choice is arbitrary and we would have obtained a different equation if we had started with a different element of W n . For instance, the elementD(λ 0 ) − → 1≤j≤n B(λ j ) would have resulted in an equally simple functional equation. Here we have restricted our attention to the analysis of (3.4) since this equation is already enough to determine the partition function. The solution of our equation is presented in Section 3.2 and is given in terms of a multiple-contour integral over L auxiliary variables. In contrast to the determinant representation obtained in [38], our integral formula offers the possibility of studying the homogeneous limit λ i → λ and µ i → µ straightforwardly. This feature seems to be of relevance for the analysis of the surface free energy of the XXZ model as discussed in [40]. It is also important to remark here that the multiple integral formula given in Theorem 2 can also be shown to satisfy the recurrence relations derived in [38]. Those recurrence relations, in addition to extra properties, are able to uniquely characterize the model partition function and thus can also be used to prove Theorem 2. However, finding an explicit representation would still demand a very non-trivial guess which is not required in our framework. In this sense the approach described here also offers a systematic way of building explicit representations.
The structure of our functional equation is further studied in Section 3.3 and we find interesting properties which are not apparent at first sight. For instance, we shown that our equation actually encodes a set of linear partial differential equations. Any single equation from this set is already able to determine the model's partition function, and thus this set is simultaneously integrated. It is worth remarking here that this property is a common feature exhibited by integrable hierarchies of differential equations.
In this work we have not analyzed the integrability of our partial differential equations in the classical sense, but that direction certainly deserves further investigation. Our construction yields a total of 2L + 1 equations, involving among others the differential operator (3.26), whose structure resembles that of a quantum many-body hamiltonian with higher-order derivatives. Although the order of the corresponding differential equation depends on L, this equation can still be reformulated as a system of first-order equations due to its linearity.
To conclude we remark here that partition functions with domain-wall boundaries and reflecting ends can also be formulated for Solid-on-Solid models as described in [48,49]. In that case the governing algebra is a dynamical version of the reflection algebra and it would be interesting to investigate if our approach can be extended to those cases.

A Properties of |0 and |0
This appendix is devoted to the derivation of formulae (2.12) arising as the eigenvalues of the operators A(λ) andD(λ) with respect to the vectors |0 and |0 defined in (2.11). For that we shall make use of (2.8) keeping in mind the representations (2.4) and (2.7). In this way the entries of (2.10) can be expressed as, recalling that κ ± (λ) = b(h ± λ). Here we are only interested in the operators A(λ) and D(λ), and from (A.1) we can see that (2.12) can be computed from the action of (2.4) on the vectors (2.11). Due to the structure of (2.2) and (2.3) we readily find the following while an analogous computation yields In their turn, the action of B(λ) andB(λ) on the vectors |0 and 0 | does not vanish but they do not correspond to eigenvectors either. Now turning our attention to the functions Λ A , ΛD andΛ A described in Section 2, we can see that Λ A can be directly read off from (A.1) and (A.2). On the other hand, the evaluation of ΛD is more involved as it corresponds to the eigenvalue of the oper-atorD(λ) = D(λ) − c(2λ) a(2λ) A(λ) with respect to the vector |0 . The latter would then require the evaluation of C(λ)B(λ) |0 as we can see from (A.1). Fortunately the Yang-Baxter algebra (2.5) can help us with that computation. Due to the unitarity property R(λ)R(−λ) = a(λ)a(−λ)1 we find the following algebraic relation

B Polynomial structure
In this appendix we prove that the partition function defined in (2.13) has the form stated in Lemma 1. More precisely, here we show that Z(λ 1 , . . . , λ L ) =Z(x 1 , . . . , whereZ(x 1 , . . . , x L ) is a polynomial of degree 2L in each one of the variables x i = e 2λ i . For that it suffices to show that B has the form . . , y ±1 L , q ±1 , t ±1 ] in the notation of Definition 3. In other words, f (2L) B (x) is a polynomial of degree 2L in the variable x, whose coefficients are products of meromorphic functions of y 1 , . . . , y L , q, t and operators on V Q . Throughout this appendix we keep track of the degree of the polynomials by indicating it in superscript as in (B.1).
The expression for B given in (A.1) reduces our task to the analysis of the dependence of κ ± , A,B, B andD with x. From (2.7) it is clear that κ ) and, therefore, it is enough to demonstrate that for a given L we have Here we consider f provide a two-dimensional representation of the U q [sl (2)] algebra obeying the commutation rules Moreover, for L = 1 the monodromy matrices (2.3) consist of a single R-matrix. Thus by writing (2.2) in the auxiliary space as we have taking into account (B.4). We can readily see from (B.7) that (B.2) holds for the case L = 1. Next we use (2.3) and (2.4) to write the following recurrence relations, Thus, if (B.2) holds for L − 1, it follows that (B.2) is true for arbitrary L. This completes the proof.

C Symmetric solutions
Here we demonstrate that any analytic solution of the functional equation (3.4) is a symmetric function. Our argument closely follows the one used in [33], although here we shall need only the first part of that argument.
As the first step of our proof we recall that the symmetric group of order L is generated by any single transposition, in addition to any cycle of length L. Thus it is enough to show that Z is invariant under cyclically permutations of λ 1 , . . . , λ k for all 1 ≤ k ≤ L in order to prove Lemma 2.

D Special zeroes
The strategy employed in Section 3.2 for solving Eq. (3.4) relies on the determination of particular zeroes of the desired solution. The location of these zeroes are stated in Lemma 3 and they are as follows: (λ 1 = µ 1 −γ, λ 2 = µ 1 ) and (λ 1 = µ 1 −γ, λ 2 = −µ 1 −γ). These specialisations of variables are given in terms of the parameter µ 1 but we could have considered any other parameter µ j instead, as will become clear from our proof.
Here we shall focus only on the first specialisation of variables, i.e. (λ 1 = µ 1 −γ, λ 2 = µ 1 ), since the same properties can be used for showing the second case. We start by noticing that the coefficients M i−1 and M i vanish for the specialisation (λ i−1 = µ 1 − γ, λ i = µ 1 ) for 2 ≤ i ≤ L, as can be seen from (3.5) and (2.12). This property is of fundamental importance for our proof. We shall first examine the cases L = 2 and L = 3 for illustrative purposes before considering the general case. L = 2. For L = 2 the functional equation (3.4) consists of three terms and it involves the spectral parameters λ 0 , λ 1 and λ 2 . Upon setting λ 1 = µ 1 − γ and λ 2 = µ 1 , two of the coefficients vanish and we are left with Here we have written · | 1,2 to denote the prescribed specialization of λ 1 and λ 2 . The remaining coefficient is nonzero for generic values of the inhomogeneities µ j and parameters γ and h. Thus we can conclude that Z(µ 1 − γ, µ 1 ) = 0.
Now we can produce L−2 additional equations by switching λ 0 ↔ λ j for 2 ≤ j ≤ L−2 as discussed in Remark 2. These equations can be written in the form

E Asymptotic behavior
The functional equation (3.4) is only able to determine the desired partition function (2.13) up to an overall multiplicative factor. In this way the full determination of Z, as defined in (2.13), requires we are able to compute it for a particular value of its variables. The asymptotic behavior stated in Lemma 4 provides us with that information and here we intend to present its proof.
Using (B.7) and writing x = e 2λ , y i = e 2µ i , q = e γ and t = e h , we find the following asymptotic behavior as x tends to infinity: The operators K and X − appearing in (E.1) were previously defined in (B.4). Also, we can see from (??) that κ ± (x) ∼ ±2 −1 t ±1 x 1 2 as x → ∞. This result combined with (A.1) and (E.1) yields the following asymptotic expansion of the operator B, where we have set From (B.5) it follows that the operators P ± j satisfy the following commutation rules: The behavior of (2.13) in the limit x i → ∞ for 1 ≤ i ≤ L can now be computed using (E.2). To this end it is convenient to introduce the operators Q (n) j := P + j q −2n + P − j q 2n (E.5) such thatZ The operators Q (n) j as defined in (E.5) satisfy the following commutation relations, as a direct consequence of (E.4). Now due to the last relation of (E.7), the summation in the right-hand side of (E.6) reduces to The expression (E.11) can be further simplified by noticing that q n ∆ + n = q −n ∆ − n . This reduces the right-hand side of (E.11) to q − L(L−1)  Here k is a constant that is fixed to be k = sinh (γ) sinh (h − µ 1 ) by the asymptotic behavior discussed in Appendix E. The solution (F.3) can still be recasted as the following contour integral, where the function H is given by (F.5) Here we have already used the explicit form of the constant k. Also, we have used some redundancies in formula (F.5) in order to make the connection with the relation (3.16) more explicit.

G Reduction of order
In Section 3.3 we have unveiled a set of linear partial differential equations underlying the functional relation (3.4). Those equations are given formally by (3.25) and the explicit construction of the set of differential operators {Ω k } was also discussed in Section 3.3. In particular, we found a compact expression for the operator Ω 2L which is given by (3.26) and (3.27). From (3.26) we see that Eq. Ω 2LZ (X 1,L ) = 0 is of order 2L and can be recasted as a system of first-order equations. The resulting system of equations is described in the following lemma. i (x 1 , . . . , x L ) for 1 ≤ i ≤ L and 1 ≤ k ≤ 2L − 1 be multivariate functions. Then the differential equation Ω 2LZ = 0 is equivalent to the following system of equations, where ∂ i := ∂ ∂x i . Proof. The verification is straightforward.
Matricial form. In order to further enhance the structure of (G.1), we finally recast our system of first-order equations as a matrix equation. For that we define the ((2L − 1)L + 1)-component vector   with functions Y i defined in (3.27) and 1 is the L × L identity matrix.