Integral formula for elliptic SOS models with domain walls and a reflecting end

In this paper we extend previous work of Galleas and the author to elliptic SOS models. We demonstrate that the dynamical reflection algebra can be exploited to obtain a functional equation characterizing the partition function of an elliptic SOS model with domain-wall boundaries and one reflecting end. Special attention is paid to the structure of the functional equation. Through this approach we find a novel multiple-integral formula for that partition function.


Introduction
Solid-on-solid (sos) models are Ising-type models for the growth of crystal surfaces. The 'solidon-solid condition' [1] forbids the surface to have overhangs, so that the shape of the surface can be modelled by associating discrete height variables to the sites of a square lattice. The surface tension is taken into account by interactions between these height variables. In this work we are only interested in the exactly-solvable subclass consisting of models that are of 'interaction-round-a-face' type, i.e. the interactions take place between the four vertices sharing a common face, and for which the height difference between all adjacent vertices is precisely one. Accordingly we take 'sos models' to mean exclusively these models, as is standard in the literature on quantum integrability. The matrix containing the face weights satisfies a 'face' (quantum) Yang-Baxter equation (ybe) guaranteeing the commutativity of the corresponding transfer matrices. The partition function is the product of the local weights, summed over all height configurations for the lattice.
Alternatively, sos models can be interpreted as 'generalized six-vertex models'. The sos weight-matrix corresponds to a 'generalized' or dynamical R-matrix that depends on a dynamical parameter keeping track of the original heights. Importantly, the dynamical R-matrix satisfies the ice rule. In addition it obeys the dynamical Yang-Baxter equation, which was first encountered by Gervais and Neveu [2] in the context of (Liouville) conformal field theory and later independently obtained by Felder [3] as the quantization of the modified classical ybe. The appropriate quantum-algebraic setting is that of elliptic quantum groups [3,4]. The ordinary six-vertex model may be recovered in the limit of infinite heights.
In the context of integrability, sos models were introduced by Baxter [5] as a convenient reformulation of the symmetric eight-vertex model. In algebraic terms [6,7] the vertex-face transformation relating the two models is a local (gauge) transformation mapping the eightvertex R-matrix to the dynamical R-matrix. By the ice rule the latter allows for the construction of a highest-weight vector (pseudovacuum), opening the way for a Bethe-ansatz analysis. Like the eight-vertex R-matrix, the associated dynamical R-matrix involves elliptic functions; this elliptic sos model is sometimes referred to as the 'eight-vertex sos' (8vsos) model. The dynamical R-matrix also appears in the quantization of the elliptic Ruijsenaars-Schneider model [8], explaining the adjective 'dynamical': here the dynamical parameter is a coordinate on the phase space of the classical integrable model.
Besides the relation with various models in statistical physics, sos models have several applications in physics and mathematics, depending on the choice of boundary conditions. The case with fixed domain-wall boundary conditions, originating in the study of scalar products of Bethe vectors [9], is of particular interest for enumerative combinatorics in the spirit of Kuperberg [10], whose work was generalized to elliptic sos models at root of unity by Rosengren [11]. Models with domain-wall boundaries further allow for concise expressions for the partition function as certain determinants. For the six-vertex model the partition function takes the form of a single determinant [12], while for sos models it can be written as a sum of 2 L determinants for models on an L × L lattice [11]. Interestingly, the presence of one reflecting (open) boundary makes it possible get a single determinant for six-vertex models [13] as well as for sos models [14,15]. Such determinant-formulas are elegant and helped uncovering striking physics like the possible dependence of thermodynamic quantities on the choice of boundary conditions imposed at finite system size in the computation [16]. Yet many applications of partition functions reside in the homogeneous limit, whose evaluation can be rather nontrivial [17], especially when elliptic functions are involved. This is one issue that makes it interesting to find alternative ways to express partition function. For sos models with domain walls at all boundaries such alternatives are available [18][19][20], and also for the six-vertex model with one reflecting end [21]. In [19][20][21] multiple-integral formulas were obtained as the solution of certain functional equations. Such repeated contour integrals might allow one to access the thermodynamic limit more straightforwardly. In the present work we extend [21] to the case of elliptic sos models with domain-wall boundaries and one reflecting end, complementing the result of [15].
Functional equations are at the core of exactly solvable models of statistical mechanics. The functional equations derived in [19][20][21] have particular origins: the dynamical Yang-Baxter algebra for [19,20] and the reflection algebra for [21]. The use of such algebras as a source of functional equations was put forward by Galleas in [22] for spectral problems and in [23] for vertex-model partition functions. Compact ways to write solutions of such functional equations were only obtained after several refinements of the method [19-21, 24, 25]. At the same time this algebraic-functional method is a beautiful example of the rigid structure imposed by the (quantum) algebra underlying exactly solvable models in statistical physics, which is reflected in the many nice properties of the functional equations obtained in this way.
Outline. This paper is organized as follows. To fix our notation and conventions we describe the elliptic sos model with domain-wall boundary conditions and one reflecting end in Section 2, setting up a graphical notation for generalized six-vertex models along the way. Our main results are stated at the start of Section 3. The first result is a linear functional equation for the model's partition function, which we derive using the algebraic-functional method in Section 3.1. Secondly we show that, despite being linear, this functional equation has a unique analytic solution, up to a constant overall factor. This is a consequence of a recursive relation between the functional equations for successive system sizes, established in Section 3.2. Our third result is a new multiple-integral formula for the partition function, which is found by solving our functional equation by recursion in Section 3.3. We conclude with some remarks and open questions in Section 4. Several technical details are left for the appendices.
The central quantity in this work is the partition function Z of an elliptic sos model on a rectangular lattice with 2L rows and L columns, and specific boundary conditions: one reflecting end and domain walls at the three other ends, as depicted in Figure 3. In this section we set up our notation and conventions, which mostly follow [19][20][21], culminating in an algebraic expression for Z. The set-up is a generalization of that in [21] including (shifts of) the dynamical parameter θ in the appropriate places. Figure 3: The lattice of an sos model with domain-wall boundaries and one reflecting end (left) and the dual lattice with the corresponding generalized six-vertex model (right).

Bulk: dynamical Yang-Baxter algebra
Let V := C e + ⊕C e − be a complex two-dimensional vector space with standard basis vectors e ± , to which we assign weights wt(e ± ) = ±1. This grading of V by weights is conveniently accounted for via the sl 2 Cartan generator on V , which we denote by h ∈ End(V ) and acts as h e ± = wt(e ± ) e ± . Write h ⊆ sl 2 for the corresponding Cartan subalgebra.
Let λ, θ, γ ∈ C be generic complex numbers. We think of the spectral parameter λ and the dynamical parameter θ as variables, and of the crossing parameter γ as fixed.
Dynamical R-matrix. The statistical weights of the height configurations in an sos model are encoded in the dynamical R-matrix R(λ, θ) ∈ End(V ⊗ V ), which must be invertible for generic λ and θ and is required to satisfy the dynamical Yang-Baxter equation on V ⊗ V ⊗ V : (2.1) Here we employ the usual tensor-leg notation, where subscripts indicate on which copies of V in the tensor product the operators act nontrivially; for example, 2) or equivalently, as a matrix acting on the third copy in V ⊗ V ⊗ V , Through the dynamical parameter each R-matrix in (2.1) is sensitive to the weights in any copy of V present on the left of the V ⊗ V on which that R-matrix acts. This can also be seen in terms of a graphical notation in which one depicts the dynamical R-matrix and its entries, labelled by indices α, β, α , β ∈ {±}, by including θ in the usual notation: The arrows serve to keep track of the orientation as we will need to rotate this picture. The dynamical Yang-Baxter equation (2.1) now acquires the usual graphical form as for the sixvertex model, but with fixed height θ in the left-most face of both sides of the equation: weight-preserving, operators on V ⊗ V . In other words, an element of End h (V ⊗ V ) satisfies the ice rule: it commutes with h 1 + h 2 ∈ V ⊗ V . The sos model that we are interested in corresponds to the elliptic solution of (2.1) lying in End h (V ⊗ V ), The entries of this dynamical R-matrix, cf. Figure 1, are given by where in turn f (λ) := −i e −iπτ /4 ϑ 1 (iλ|τ )/2 is essentially the odd Jacobi theta function with elliptic nome e iπτ ∈ C satisfying Im(τ ) > 0. The properties of f required for this work are collected in Appendix A.1. The algebraic structure underlying (2.6)-(2.7) is the elliptic quantum group E τ,γ (sl 2 ) whose properties are discussed in [3,4,28]. To make contact with the ordinary six-vertex model one takes the trigonometric limit, lim τ →i∞ f (λ) = sinh(λ), and subsequently lets θ → ∞, which the yields the [U q ( sl 2 )-invariant form of the] six-vertex model's R-matrix considered in [19,21]. We will also need the following unitarity property of the R-matrix: Notice that the proportionality factor in (2.8) does not depend on the dynamical parameter. Finally we record that transposing the first leg of the R-matrix yields an operator satisfying the ice rule [−h 1 + h 2 , R t 1 12 (λ, θ)] = 0 and the crossing symmetry Here σ y ∈ End(V ) is the second Pauli matrix, and the colons indicates that the h 1 is taken to act after (on the left of) the transposed R-matrix.
Monodromy matrices. A central role in the quantum inverse-scattering method (qism) [6] is played by the monodromy matrix. We are interested in its dynamical version, introduced in [7], which is built out of solutions of (2.1) lying in End h (V ⊗ V ). Focus on any single row of the lattice, with associated 'auxiliary space' V 0 ∼ = V and spectral parameter λ ∈ C. Every column 1 ≤ j ≤ L comes with a 'local quantum space' V j ∼ = V . The 'global quantum space' is defined as the L-fold tensor product W := V 1 ⊗ V 2 ⊗ · · · ⊗ V L . Since (2.1) depends only on differences of spectral parameters we can easily incorporate a (fixed) inhomogeneity parameter µ j ∈ C for each V j , 1 ≤ j ≤ L.
There are two monodromy matrices T 0 (λ, θ),T 0 (λ, θ) ∈ End(V 0 ⊗ W ), defined as oppositely ordered products (2.10) The operator H := L j=1 h j , representing the generator of the Cartan subalgebra h ⊆ sl 2 on W , endows W with a grading by weights: Usually only one of the monodromy matrices in (2.10) is considered since by (2.8) the two are essentially inverse to each other: The two monodromy matrices are also related by crossing symmetry, following from the ice rule for the transposed R-matrix and (2.9), For our purposes, however, it will be convenient to include both T andT .
Dynamical Yang-Baxter algebra. The monodromy matrices (2.10) can be viewed as matrices in auxiliary space with entries acting on W : Together with the ice rule, the dynamical Yang-Baxter equation (2.1) implies that T obeys the quadratic relation 14) It follows that the entries of T as in (2.13) furnish W with a representation of the operator algebra associated with E τ,γ (sl 2 ), known as the (dynamical) Yang-Baxter algebra, which we denote by A = A τ,γ (sl 2 ). Since T 0 (λ, θ) ∈ End h (V 0 ⊗ W ) preserves the grading by weights of V 0 ⊗ W , and not that of W , the generators on the anti-diagonal of (2.13) are not h-invariant: Nevertheless the dynamical Yang-Baxter algebra is compatible with the grading of W when A itself is viewed as a graded algebra: B and C have weights −2 and 2 respectively, while A and D have weight zero. The generators contained inT are analogously seen to have weights wt(B) = −2, wt(C) = 2 and wt(Ā) = wt(D) = 0. The monodromy matrixT obeys This relation is equivalent to (2.14) due to the unitarity relations (2.8) and (2.11). In particular, the entries ofT as in (2.13) subject to the relations encoded in (2.16) also provide a representation of A. From (2.11) and (2.14) it follows that the generators A, B, C, D andĀ,B,C,D are connected by (2.11) as well as by the relations contained in

Reflecting end: dynamical reflection algebra
One can regard the Yang-Baxter algebra, or its dynamical version, as governing the bulk of the lattice. The reflecting boundary is described by an operator K 0 (λ, θ) ∈ End(V 0 ) satisfying the reflection equation, also known as boundary Yang-Baxter equation, Although we consider a dynamical R-matrix, notice that there are no shifts in the dynamical parameter: this is just the ordinary reflection equation. We are interested in the solution of (2.18) that lies in End h (V 0 ), which was obtained in [29] and reads depending on a fixed boundary parameter ζ ∈ C. These quantities can be depicted as The dynamical reflection equation (2.18) then becomes (2.23) As in [30] one can use the relations (2.14)-(2.17) for A and the reflection equation (2.18) to show that this operator satisfies the following quadratic relation: Note that, unlike (2.14) and (2.16)-(2.17), the relation (2.24) does not involve shifts in the dynamical parameters of any monodromy matrix. This permits us to think of θ as a constant and suppress the dependence on this parameter when writing (2.23) as a matrix in auxiliary space with entries in End(W ): Dynamical reflection algebra. The entries of (2.25) together with H, subject to the relations contained in (2.24) and the ice rule for T , endow W with a representation of the (dynamical) reflection algebra, which we denote by B. For our purposes it is convenient to follow [30] and replace the generator D of B by the combinatioñ Just like for A, the ice rule for the double-row monodromy matrix endows B with a grading by weights: wt(B) = −2, wt(C) = 2 and wt(A) = wt(D) = wt(D) = 0. That is, (2.15) also holds for A, B, C,D and it follows that [Thus, instead of including H as one of the generators of B, one may equivalently view B as a graded algebra, with grading by weights.] Finally, with the help of (2.17) and (2.27) one can check that (2.12) implies that B enjoys the crossing symmetry

Domain walls and the partition function
The final ingredient for the partition function are the domain-wall boundary conditions on the three remaining boundaries, which need to be chosen such that Z is nontrivial. We consider the case first introduced for vertex models in [9] and extended to sos models in [16]. These boundary conditions are included in the qism with the help of the following vectors in W .
Pseudovacua. The highest-and lowest-weight vectors |0 , |0 ∈ W are defined as These are also eigenvectors for the entries on the diagonal of the double-row monodromy matrix (2.25). We will need the following 'vacuum' eigenvalues: In Appendix B we show that these eigenvalues are given by Partition function. For models with toroidal boundary conditions one is usually able to express the partition function as a trace of a product of transfer matrices, so that the evaluation of the partition function becomes an eigenvalue problem. This is not the case when other types of boundary conditions are considered. Vertex models with domain-wall boundary conditions [9] still admit an algebraic formulation, yet evaluating the partition function is not a priori an eigenvalue problem. In the case of domain-wall boundaries, possibly with a reflecting end, the partition function can be written as the expectation value of a product of operators. For sos models the situation is analogous, and the partition function of the inhomogeneous sos model with domain walls and one reflecting end, as in Figure 3, is given by It depends on all parameters of the model: L spectral parameters λ = (λ 1 , · · · , λ L ), the dynamical parameter θ, L inhomogeneities µ j , the crossing parameter γ, the boundary parameter ζ, and the elliptic nome e iπτ . (Throughout this work vectors in C L are printed in boldface.) In the graphical notation (2.33) is represented by the diagram on the right in Figure 3.
To conclude this section we collect the properties of Z that are needed to characterize the partition function in our approach. Tracing back the dependence on the spectral and dynamical parameters it is clear that Z, being polynomial in the statistical weight (2.7) and boundary weights (2.20), is meromorphic. Moreover, up to an overall factor originating the denominator in (2.20), it is entire in the spectral parameters: Lemma 1 (Polynomial structure). The partition function (2.33) can be written as where the normalized partition functionZ, viewed as a function of any single λ j , is a theta function of order 2(L + 1) and norm (L − 1)γ.
The basics of higher-order theta functions, which are also known as elliptic polynomials, are summarized in Appendix A.2. For a proof of Lemma 1 we refer to [15], in which a factorizing Drinfel'd twist [31] is used to determine the dependence of the partition function on the spectral parameters. In fact, in [15] it is shown thatZ can be further written as L j=1 [2λ j ] times a higherorder theta function of order 2(L − 1) and norm (L − 1)γ in each variable, but (2.34) suffices for our purposes.
Lemma 2 (Symmetric function). The partition function (2.33) is symmetric in the L spectral parameters contained in λ.
Proof. The reflection algebra (2.24) includes the following relation, obtained using (2.27): The symmetry of (2.33) is thus obvious.
Although we will not need it let us further mention that the partition function is also symmetric in the inhomogeneity parameters µ j [15].
Remark 1 (Crossing symmetry). Let us finally record that the partition function (2.33) inherits crossing symmetry (2.28) of B: for any 1 ≤ j ≤ L it satisfies This property is not a necessary ingredient for our approach, yet it will be useful to get a better understanding of certain issues that we will encounter along the way.

Functional equations
In this section we use the dynamical reflection algebra to obtain a functional equation that determines the partition function (2.33). Our main results are as follows.
Theorem 1 (Functional equation). The partition function (2.33) of the sos model with one reflecting end and domain-wall boundaries obeys the linear functional equation where the caret indicates that the ν'th spectral parameter is omitted, and the coefficients M ν feature the vacuum eigenvalues from (2.31): In the terminology of [32], (3.1) is a cyclic linear functional equation. Observe that it features L + 1 variables whilst the partition function (2.33) depends on only L spectral parameters. In analogy with the conventions from field theory in this work we use the Latin alphabet for indices ranging through {1, 2, · · · , L} and Greek for indices in {0, 1, 2, · · · , L}.
Theorem 1 asserts the existence of solutions to the functional equation (3.1). In view of the linearity of (3.1) the uniqueness of its solutions has to be investigated as well. This issue is addressed by the following result.
Theorem 2 (Uniqueness). Up to an overall λ-independent normalization factor, (3.1) has a unique solution within the class of analytic functions on C L .
In fact, this uniqueness is established by exploiting the presence of the extra variable λ 0 to find a recursion relation between the functional equations for length L and L−1. This recursion moreover allows us to solve (3.1) to obtain a multiple-integral formula for the partition function: Theorem 3 (Multiple-integral formula). The partition function (2.33) can be written as where Ω L is a constant (i.e. λ-independent) normalization factor, the integration contours Γ λ consist of small counter-clockwise loops enclosing only the poles of (3.3) at λ j , 1 ≤ j ≤ L, and the functions m l read (3.4) Direct inspection further suggests that the normalization constant in (3.3) is given by where L/2 is the integer part of L/2. Let us pause for a moment to look at the analytic structure of (3.3)-(3.4). Although (3.3) has apparent poles at z l = −γ/2, these simple poles are removable as the residues of m l (z 1 , · · · , z l )/[2z l + γ] at these points vanish. The simple poles in (3.4) at z l = z j cancel against the corresponding zeroes in (3.3), while the poles at z l = −z j − γ are related to those at z l = z j by crossing symmetry, see Remark 1 above and Lemma 3 below. Finally the simple poles of (3.4) at z l = −θ − ζ matches Lemma 1.
Trigonometric degeneration. Before turning to the proofs let us briefly discuss the limit τ → i∞ in which f (λ) ∝ ϑ 1 (iλ|τ ) degenerates into a trigonometric function: f (λ) → sinh(λ), cf. Appendix A.1. The trigonometric sos model with a diagonal reflecting end arises naturally by applying a vertex-face transformation to the six-vertex model with nondiagonal reflection [33]. The partition function of this model was studied in [14], where it was expressed as a determinant. The set-up of Section 2 and all our results carry through upon reinterpreting [λ] = sinh(λ). In Lemma 1 the notion of elliptic polynomials has to be replaced by that of trigonometric (or Fourier) polynomials, which may equivalently be described as polynomials in x = e 2λ up to an overall factor of e −λ . The degree of the limiting trigonometric polynomial is at most equal to the order of an elliptic polynomial, though it can be less as the example [2λ] illustrates.
When we further take the non-dynamical limit θ → ∞ our functional equation (3.1) and its solution (3.3)-(3.4), as well as the normalization (3.5), boil down to the results of [21]. Indeed, with our conventions all factors involving θ come in ratios, and the expressions of [21] are heuristically recovered from those presented here by simply dropping all such ratios. [To match (3.4) see Remark 3 at the end of Section 3.3 below.]
To keep track of the number of generators A, B, C,D that constitute any element of the dynamical reflection algebra we consider a grading B = n≥0 B (n) in addition to the grading by weights described in Section 2.2. This time we simply declare the generators on the righthand side of (2.25) to have degree one. Thus, the degree-n part B (n) ⊆ B is the linear span of products of precisely n generators A, B, C,D. Note that this grading by number is compatible with the defining relations in (2.24) since those only involve terms of degree two.
Proof of Theorem 1. The idea for obtaining our functional equation for the partition function (2.33) is to use (2.30) to insert an A on one side of the product of B's in (2.33) and use the reflection-algebra relations to move it to the other side.
To save space let us indicate the arguments of the generators of B as subscripts, e.g. A 0 := A(λ 0 ). We start by finding the appropriate algebraic relations in the dynamical reflection algebra. For this we need (2.35) together with the relations in B (2) ⊆ B obtained from (2.24) using (2.27):

6)
whereD was defined in (2.26) such that no term proportional to B 1 A 0 appears on the righthand side of (3.7). By repeated application (2.24) gives rise to relations in B (n) . In particular, (2.35) allows us to suppress the harpoon in (2.33), and from (3.6)-(3.7) it follows that in degree n = L + 1 we have the following relation: Now we multiply this relation from the left by 0 | and from the right by |0 . This may be described more formally as follows. Consider the linear map π : B −→ C[λ ±1 0 , · · · , λ ±1 L ] assigning multivariate meromorphic functions to elements of the reflection algebra defined by for the restriction of π to the degree-n part of B. Apply π L+1 to (3.8) and use (2.29)-(2.30) to express the result in terms of π L . This yields the desired functional equation with coefficients (3.2).
A few comments are in order. Firstly, rather than inserting an A(λ 0 ) in (2.33), one could equally well useD(λ 0 ). In [25] the corresponding functional equations are referred to as being of 'type a' or 'type d'. From the structural similarity between (3.6) and (3.7) it is clear that the equation of type d is also of the form (3.1) but with different coefficients. In the following section we will see that, like in [21], the functional equation (3.1)-(3.2) already suffices to characterize the partition function, and we do not need the explicit form of the functional equation of type d.

Reduction and uniqueness
In this and the next section we solve our functional equation. Our approach elaborates on the analysis developed in [19][20][21]. As a byproduct we will be able prove Theorem 2. Let us reserve the symbol 'Z' for the partition function (2.33) and study the linear functional equation L ν=0 M ν (λ 0 ; λ) Z(λ 0 , · · · , λ ν , · · · , λ L ) = 0 (3.12) with coefficients (3.2). In view of Lemmas 1 and 2 we are looking for solutions Z to this equation in the space of symmetric higher-order theta functions (up to an overall factor) on C L . We begin by collecting some useful properties of the functional equation (3.12). First, the properties from Lemmas 1 and 2 (and Remark 1) are in fact automatic when we look for sufficiently nice solutions: Lemma 4 (Properties of analytic solutions). Any analytic solution Z of (3.12) has the following properties: i) The normalized solutionZ is a theta function of order 2(L + 1) and norm (L − 1)γ with respect to each λ j .
ii) Z is symmetric in the L spectral parameters; iii) Z satisfies the crossing symmetry (3.14) The proofs can be found in Appendices C.1-C.3. Next, since (3.12) involves one more spectral parameter than Z depends on, we can specialize any single spectral parameter to any value we like. Together, (2.31) and (3.2) show that the greatest simplification of (3.12) occurs when we set λ 0 = ±µ k − γ for some 1 ≤ k ≤ L. Indeed, M 0 (±µ k − γ; λ) =Λ A (±µ k − γ) consists of a single λ-independent product that is nonzero for generic values of the parameters. This allows us to express Z(λ) as a linear combination of Z's that each depend on only L − 1 free spectral parameters. In particular, when L = 1 this procedure completely determines the solution up to an overall scale: Lemma 5 (Case L = 1). For L = 1 any solution to the functional equation (3.12) can be written as for some normalization Ω 1 not involving λ, where m 1 was defined in (3.4).
[The second equality is recorded anticipating Section 3.3.] For L ≥ 2 we can further reduce the expression for Z(λ) by exploiting the following vanishing property for Z.
Following [19][20][21] we refer to these 4L combinations as 'special zeroes' of the solution. Note that for given k the four special zeroes are precisely related by crossing one or both values, µ → −µ − γ, cf. Lemma 4 (iii). Lemma 6 is substantiated in Appendix C.4. Of course, by part (ii) of Lemma 4, an analytic solution thus vanishes when any two of its arguments are specialized to λ + and λ − . Together with parts (i) and (iii) of Lemma 4 this has the following consequence.
Corollary 1. In the setting of Lemma 6 suppose that Z is an analytic solution of (3.12). Then we have for any i ∈ {1, · · · , L} that where theZ ± are analytic functions in L − 1 spectral parameters with the properties from Lemma 4, now with L − 1 instead of L.
Although the preceding holds for any choice λ i ∈ {±µ k − γ, ∓µ k } with any 1 ≤ k ≤ L, only when k = 1 or k = L it is possible to reduce the problem of solving (3.12) for length L to that for L − 1. The results for these two choices for k turn out to be the same, cf. the proof of Theorem 2 and Remarks 2 and 3 below, yet the case k = L is easier to treat. The four resulting special zeroes can be treated simultaneously, yet for now we consider the two choices (µ L − γ, λ * ) with λ * ∈ {−µ L − γ, µ L } to avoid cumbersome notation. This is consistent with crossing symmetry, implying that the result should not depend on this choice at any rate. We obtain the following recursive relation.

Proposition 1 (Reduction).
If Z is an analytic solution of the functional equation (3.12) for length L then we can write (3.17) whereZ is an analytic solution of (3.12) for length L − 1, and the proportionality constant Ω is independent of λ.
Proof. Let us employ the following shorthand for the arguments of various functions. We write e.g. λ 0,1,··· for λ 0 , λ 1 , · · · . The omission of a spectral parameter is indicated by a caret on the corresponding index; the presence of λ * is shown by an asterisk. Suppose that Z satisfies (3.12) for L. The idea is to use suitable specializations of λ 0 and λ L . In view of the discussion preceding Lemma 5 we first specialize λ 0 = µ L − γ in (3.12) and solve for Z(λ) to get (3.20) One can check that the two possible choices of λ * only lead to a constant (λ-independent) overall (ν-independent) factor for the coefficients (3.20). Since this factor can be dropped in (3.19) we have not indicated the dependence of theM ν on this choice. Clearly the reduced functional equation (3.19) has the same structure as our original functional equation when we take L − 1 instead of L in (3.12); the coefficients (3.20) also exhibit the symmetries from Lemma 3 for length L − 1. Although (3.19) for L = 2 happens to be proportional to (3.2) with L = 1, unlike in earlier works such as [19] (but like in [21]), this simple relation does unfortunately not persist for L ≥ 3.
At the end of Section 3.1 we noticed that interchanging λ 0 ↔ λ j in (3.12) for some 1 ≤ j ≤ L yields another functional equation. Thus, for length L − 1, in fact we have L different functional equations: (3.12) together with the L−1 equations obtained by switching λ 0 ↔ λ j . (Recall that not all of these equations are linearly independent.) We claim that for any length L the left-hand side of the reduced functional equation (3.19) can be written as some linear combination of the left-hand sides of these L functional equations corresponding to length L − 1. To see that this is indeed the case we form the L × L matrix with entries given by (3.11) for length L − 1, except that the last row is replaced by the coefficients (3.20) of the reduced equation. Our claim is true if this matrix has zero determinant. Whilst we have not managed to prove this rigorously, analytic and numerical investigations for L ≤ 7 confirm that this is indeed the case for these system sizes, and we see no obstruction for this pattern to continue for larger L too.
Due to crossing symmetry, cf. Lemma 4 (iii), the other special zeroes, (−µ L − γ, λ * ) with λ * ∈ {µ L − γ, −µ L }, yield the same result. Another way to understand this is the fact that the proof of Proposition 1 does not involve any algebraic choices, unlike e.g. for the equation λ 2 = 1, so that the procedure only has one possible result. (It may be instructive to convince oneself of this for L = 1 and L = 2.) Correspondingly it is now easy to prove the uniqueness, up to a constant overall factor, of analytic solutions of (3.12).
Proof of Theorem 2. We use induction on L. The base case is furnished by Lemma 5. Suppose that Z is an analytic solution of (3.12) for length L. Then by Proposition 1 the functionZ in (3.17) is analytic and solves the equation for length L − 1. Hence, according to the induction hypothesis,Z is unique up to a constant normalization factor. But (3.17) determines Z in terms ofZ, again up to a constant factor. This proves the theorem.
Note that the proof crucially depends on the recursion between the functional equation for successive lengths. In particular it also applies to the equations derived in [19][20][21].

Remark 2.
In view of the preceding proof one may wonder about the dependence on the choice k = L in Proposition 1. The only other choice of k that allows one to relate the reduced functional equation (3.19) to the equations for length one lower is k = 1. In that case one finds that any analytic solution of (3.12) for length L can be written in terms of a solutionŽ of (3.12) for length L − 1 and with inhomogeneities (µ 2 , · · · , µ L ) instead of (µ 1 , · · · , µ L−1 ): (3.21) At the end of the next section we will see that the result of this alternative approach is the same as that for k = L, in accordance with Theorem 2.

Multiple-integral formula
In this section we solve our functional equation (3.12) by induction on L, based on the analysis of the previous section and Proposition 1 in particular, culminating in the multiple-integral formula from Theorem 3. As in the proof of Theorem 2 we proceed by induction on L. Iterating the relation from Proposition 1 we can find a closed expression for the solution to (3.12), which by Theorem 2 is unique up to an overall constant factor. Proposition 2 (Solution as symmetrized sum). Up to a constant overall normalization factor, the solution to the functional equation (3.12) for L ≥ 2 can be written as the following symmetrized sum:

22)
where S L denotes the symmetric group in L symbols, Z(λ) is given by (3.15), the factor M l is to be understood as given by (3.2) for length l, and m l was defined in (3.4).
Proof. The second equality follows directly from (2.31), (3.2) and (3.15). The proof of the first equality is by induction on L. For L = 2 the statement follows directly from Proposition 1. The inductive step is straightforward, using Proposition 1, the properties of the coefficients M l recorded at the end of Section 3.1, and the bijection of labelling sets {1, · · · , L} × S L−1 ∼ − − → S L given by (i, σ) −→ (i, i + 1, · · · , L) • σ , where σ ∈ S L is the extension of σ fixing L.
The remaining step in the proof of Theorem 3 is based on the following trick, which appears to be common lore.
Lemma 7 (Multiple-integral formula for symmetrized sums). Consider L ≥ 1 and let λ ∈ C L be generic. Suppose that g : C L → C is a meromorphic function that is regular in a neighbourhood of z i = λ j for all i, j ∈ {1, · · · , L}, and that f : C → C is analytic in a neighbourhood of the origin and satisfies f (0) = 0 = f (0). Then we can write

23)
where each z i is integrated over the contour Γ λ consisting of small counter-clockwise oriented loops around all the λ j , 1 ≤ j ≤ L.
Proof. Again one proceeds by induction on L. For L = 1 the statement follows immediately from Cauchy's residue theorem for the single pole at z = λ, with residue 1/f (0). For L ≥ 2 we assume that λ ∈ C L is such that f (λ i − λ j ) = 0 for all i = j; in particular this means that the components λ j should all be distinct. The inductive step entails applying the residue theorem to integrate over z L , then employing the induction hypothesis (3.23) to the L functions g(z)| z L =λ i L−1 j=1 f (λ i − z j ), and finally using {1, · · · , L} × S L−1 ∼ − − → S L as in the proof of Proposition 2 above.
It is now easy to prove our final result.
Proof of Theorem 3. The function f (λ) ∝ ϑ 1 (iλ|τ ) and the function g(λ) determined by the second line in (3.22) satisfy the assumptions of Lemma 7, allowing us to pass to repeated contour integrals.
Remark 3. If one proceeds along the lines of Remark 2 rather than Proposition 1 one finds an equivalent symmetrized sum: Here the lth factor of M 1 is understood to be given by (3.2) for length L − l + 1 and with all inhomogeneity parameters shifted as indicated. One can verify that this expression is equal to (3.22). The resulting representation of the partition function as a multiple contour integral does look slightly different from (3.3), and is immediately seen to generalize the result of [21].

Conclusion
In this work we consider the elliptic sos model with domain-wall boundaries and one reflecting end. Using the dynamical reflection algebra we derive a cyclic linear functional equation for the model's partition function (Theorem 1). By construction this functional equation has a solution.
We also prove its uniqueness, up to normalization, within the class of analytic functions on C L (Theorem 2). This allows us to obtain a closed expression for the partition function as a symmetrized sum (Proposition 2), which can be concisely rewritten as a multiple contour integral (Theorem 3). By taking the trigonometric limit analogous results are obtained for the trigonometric sos model on the same lattice. Along the way special attention is paid to the properties of our equation. Our results generalize those of Galleas and the author [21]. The equation that we obtain is by no means the only one that can be found in this way. Indeed, as we mention in Section 3.1 our functional equation may be said to be of 'type a'. One can similarly derive functional equations of 'type c', which we expect to be significantly more involved (cf. [19]), or of 'type d'. In each case the form is as in (3.1) but with coefficients different from (3.2). Since as in [21], but unlike in e.g. [25], our functional equation alone already determines the partition function (Theorem 2) we do not pursue these other possibilities.
A common feature of our functional equation and others obtained via the algebraicfunctional method [19][20][21][23][24][25] is its structure (3.1): it can be described as a cyclic linear functional equation [32]. Of course the coefficients, presently (3.2), differ from case to case. Another notable difference with the functional equation from [20] is that, due to the reflecting boundary, our equation does not determine the behaviour of the partition function on the dynamical parameter θ. This can already be anticipated at the algebraic level by comparing the θ-dependence of the monodromy matrices in the defining relations (2.14) of the dynamical Yang-Baxter algebra A = A τ,γ (sl 2 ) and (2.24) of the dynamical reflection algebra B ⊆ A.
There are several problems left that we do not address in this work. These include the homogeneous limit, in which our expressions for the partition function might be easier to handle than the determinant formulas obtained by Filali and Kitanine [14,15], the thermodynamic limit, which has been studied for the ordinary six-vertex model with the same boundary conditions in [34], and the relation to combinatorial problems at the ice-like point γ = iπ/3 [10,11]. For the trigonometric sos model it should also be quite straightforward to convert our functional equation into a family of partial differential equations as in [21], which might give new insights into the symmetries of the partition function. Another possible direction is to investigate if this method can also be applied to models associated with Lie algebras sl N of higher rank, or even Lie superalgebras. It would be interesting to return to such issues in the future. The series in (A.1) converge absolutely when the elliptic nome e iπτ satisfies | e iπτ | < 1. When e iπτ tends to zero the theta function degenerates into a trigonometric function, and in fact lim τ →i∞ f (λ) = sinh(λ).
In this work we often use the short-hand [λ] := f (λ) and [λ 1 , λ 2 , · · · ] := [λ 1 ] [λ 2 ] · · · . The properties that we need to work with the special function f are as follows: i) Analytic structure: f is entire (holomorphic on all of C) and only has simple zeroes; ii) Double quasiperiodicity: iii) Oddness: iv) Addition rule: Properties (i)-(iii) readily follow from the series (A.1). The zeroes of f form a lattice, iπZ + iπτ Z ⊆ C. Using Liouville's theorem in complex analysis it is easy to see that (i)-(iii) imply property (iv). Reversely, property (iv) implies (iii) and, up to a constant factor, (i) and (iv) uniquely characterize f along with its trigonometric and rational limiting cases.

A.2 Higher-order theta functions
Fix τ ∈ C with Im(τ ) > 0 and consider the function f from (A.1). For N ∈ N 0 and t ∈ C one defines a theta function of order N and norm t to be a complex function F (λ) for which there exist numbers Ω, t 1 , · · · , t N ∈ C with N n=1 t n = t such that F can be written in the factorized form Let Θ N,t be the set of theta functions of order N and norm t with respect to the variable λ. A classic result [36, §15] is that F (λ) ∈ Θ N,t if and only if F (λ) is entire and doubly quasiperiodic with quasiperiods iπ and iπτ such that For example one can verify that F (λ) := [nλ + γ] lies in Θ n 2 ,nγ with respect to λ. As a corollary of (A.6) we see that Θ N,t is a vector space: any linear combination of functions in Θ N,t also satisfies (A.6). This factorization property for higher-order theta functions is very useful.
For completeness we also mention that when N ≥ 2 the dimension of Θ N,t is equal to N , and that there is an interpolation formula expressing F (λ) ∈ Θ N,t in terms of its values at N generic points λ n ∈ C: Further details can be found in [18,37].

B Vacuum eigenvalues
In this appendix we ascertain that the pseudovacua ( Thus we first compute the action of the generators of the dynamical Yang-Baxter algebra A on the vectors (2.29). Due to (2.6) and (2.7), |0 is a simultaneous eigenvector of A, C, D and A,C,D, with corresponding eigenvalues Likewise 0 | is an eigenvector of these operators, with eigenvalues Combining these with (B.1) we directly find that Λ A is given by the expression in (2.31). In contrast, neither pseudovacua is an eigenvector of B orB. This prevents a simple evaluation of CB|0 needed for Λ D in view of (B.1). This issue can be circumvented by using the relation (2.17), and (2.27), to rewrite the problematic term as C(λ, θ)B(λ, θ) =B(λ, θ + γ) C(λ, θ + γ) Together with (B.1) and (B.2) this yields the result for ΛD from (2.31), whereD was defined in (2.26) and we also used the addition rule (A.4) to rewrite the prefactor. ForΛ A we proceed analogously. The evaluation of 0 |BC is avoided by exploiting the following relation contained in (2.17):

C Properties of solutions
In this appendix we derive various properties of solutions of our functional equation.

C.1 Polynomial structure
Here we investigate the analytic structure of any given solution of the functional equation (3.12).
In view of Lemma 1 we expect a simple pole at −θ − ζ. Part (i) of Lemma 4 states that up to these poles, the solution is a higher-order theta function whose order and norm in the spectral parameters are determined by the functional equation.

(C.2)
First we focus on the dependence on λ 0 . The explicit form of the coefficients, see (2.31) and (3.2), show thatM 0 has no poles in λ 0 whilst (for generic λ) theM i only have simple zeroes in λ 0 . ThereforeZ is entire in λ 0 . Next one checks thatM 0 is a theta function of order 4L + 6 and norm (L + 2)γ − θ whereas eachM i is a theta function of order 2L + 4 and norm 3γ − θ. Thus comparing the λ 0 -dependence of the terms in (C.2) we conclude thatZ must be a theta function of order 2(L + 1) and norm (L − 1)γ in λ 0 .
One proceeds likewise for the dependence on λ i for any fixed 1 ≤ i ≤ L. With respect to this variableM i is a theta function of order 4L + 6 and norm (2L + 1)γ while eachM ν with ν = i is a theta function of order 2L + 4 and norm (L + 2)γ. We conclude that Z is also a theta function of order 2(L + 1) and norm (L − 1)γ in the λ i .
It is worth pointing out that, unlike for the functional equation in [19,20], similarly focussing on the dependence on the dynamical parameter θ in (C.2) does not give us more information: theM ν are theta functions of the same order and norm in θ for each ν ∈ {0, 1, · · · , L}. This is in accordance with our notation, in which we do not indicate the dependence on θ of the coefficients or the solution (or partition function).

C.2 Symmetric solutions
In the proof of Theorem 1 we used the symmetry of the partition function Z from (2.33), cf. Lemma 2, to derive our functional equation. Here we show that although (3.12) is not manifestly symmetric in all spectral parameters, cf. the end of Section 3.1, reversely any analytic solution to (3.12) is symmetric.
Proof of Lemma 4 (ii). Consider the functional equation (3.12) in the limit λ 0 → λ i for some 1 ≤ i ≤ L. When Z is analytic so that part (i) of Lemma 4 applies, (3.2) shows that the only singularities in (3.12) are the simple poles in M 0 and M i , with opposite residues Computing the residue of (3.12) as λ 0 → λ i we thus obtain where λ σ i = (λ i , λ 1 , · · · , λ i−1 , λ i+1 , · · · , λ L ) is the image of λ under the (right) action of the permutation σ i = (i, i − 1, · · · , 2, 1) ∈ S L : its jth component is (λ σ i ) j = λ σ i (j) . Since (C.3) is generically nonzero and Z is analytic, (C.4) implies that Z(λ) = Z(λ σ i ). Hence Z is invariant under the cycle σ i of length i. Now any two cycles of length 2 and L already generate the full permutation group S L . Using the above argument for i = 2 and i = L we may thus conclude that Z is invariant under the action of S L permuting the variables λ j , as we wanted to show.

C.3 Crossing symmetry
The final part of Lemma 4 states that all solutions of the functional equation enjoy crossing symmetry. The proof is analogous to that of part (ii) above.
Proof of Lemma 4 (iii). We investigate (3.12) as λ 0 → −λ i − γ for some 1 ≤ i ≤ L. If the solution Z is analytic, by part (i) of Lemma 4 the sole singularities in (3.12) are again the simple poles in M 0 and M i , with residues given by (C.5) Here we used (2.32). Using the symmetry of Z from part (ii) of Lemma 4 the residue of (3.12) As (C.5) is generically nonzero and Z is analytic, (C.6) implies that Z is invariant under crossing λ i → −λ i − γ, which is what we had to demonstrate.
Our task is to show that the determinant of this matrix is nonzero. Although we have not been able to prove this rigorously, analytic and numerical inspection for L ≤ 9 reveal that this is indeed the case. We see no reason to doubt that this is the case for larger L as well.