Density matrices in integrable face models

Using the properties of the local Boltzmann weights of integrable interaction-round-a-face (IRF or face) models we express local operators in terms of generalized transfer matrices. This allows for the derivation of discrete functional equations for the reduced density matrices in inhomogeneous generalizations of these models. We apply these equations to study the density matrices for IRF models of various solid-on-solid type and quantum chains of non-Abelian $\bm{su(2)_3}$ or Fibonacci anyons. Similar as in the six vertex model we find that reduced density matrices for a sequence of consecutive sites can be 'factorized', i.e.\ expressed in terms of nearest-neighbour correlators with coefficients which are independent of the model parameters. Explicit expressions are provided for correlation functions on up to three neighbouring sites.


Introduction
Exact Bethe ansatz solutions of integrable lattice models provide valuable insights into properties which can be related to their spectrum such as thermodynamic properties and the nature of low energy excitations. The computation of general correlation functions in this framework is much more involved. For certain integrable vertex models and in particular the spin-1/2 Heisenberg model, however, manageable expressions for correlation functions have been obtained using methods based on the representation theory of quantum algebras, functional equations of q-Knizhnik-Zamolodchikov (qKZ) type, or the algebraic Bethe ansatz [1][2][3][4][5].
Considering inhomogeneous generalizations of these models a remarkable property of the corresponding reduced density matrices has been established: in Refs. [6][7][8] it was found that correlation functions of spins on N consecutive sites in the ground state of the infinite length antiferromagnetic Heisenberg chain are solutions of the qKZ equation (or a reduced version thereof) and can be expressed as sums of terms factorizing into products of nearest neighbour (two-point) functions of the generalized models. Their coefficients are recursively defined elementary functions of N spectral parameters and do not depend on model parameters such as the system size or choice of inhomogeneities. While this approach based on the qKZ equations is limited to infinite chains the factorization was also observed in studies of the three-site reduced density matrix for isotropic chains of finite length (or at finite temperature) [9,10]. Later, using the fermionic structure in the space of operators of the XXZ model [11,12], this property has been proven to hold for arbitrary correlation functions of the XXZ Heisenberg chain in an external magnetic field and at finite temperature [13].
In another approach to the computation of correlation functions of the Heisenberg chain discrete functional equations of reduced qKZ-type have been derived starting from the local properties of the six vertex model [14,15]. These equations can be shown to characterize the generalized N-site reduced density matrix D N of the model uniquely when complemented with an asymptotic reduction relating the latter to D N−1 when one of the spectral parameters is sent to infinity [15]. Taking the factorized form of the reduced density matrix from Ref. [7] as an ansatz it is found that the latter does indeed satisfy these relations. Therefore, this approach provides an alternative, though not constructive proof of the factorization property.
Less is known for integrable interaction-round-a-face (IRF) models. For an important class of IRF models, the Andrews-Baxter-Forrester (ABF) series of solid-on-solid (SOS) models, Baxter's corner transfer matrix [16] has been used to compute local height probabilities in the infinite lattice [17]. Vertex operators for the ABF models introduced in this approach can be related to representations of quantum group symmetries present in the infinite system and their correlation functions satisfy qKZ equations [18] and the algebra formed by the vertex operators has been bosonized to obtain integral representations for multi-point local height probabilities in restricted SOS (or RSOS) models in the thermodynamic limit [19]. As for the vertex models the integrability of IRF models is based on a Yang-Baxter equation satisfied by the Boltzmann weights. For the SOS models these weights can be arranged in an R-matrix depending on an additional dynamical parameter. Due to this formulation certain aspects of the Quantum Inverse Scattering Method can be employed for the analysis of these models: the corresponding modified (dynamical) Yang-Baxter algebra allows for the solution of the spectral problem by means of an algebraic Bethe ansatz [20] or an adaption [21] of the framework of Sklyanin's Separation of Variables [22]. Similarly, the inverse problem relating local spin operators to elements of the Yang-Baxter algebra [23,24] has been solved for the dynamical vertex models and also allows to express local heights in this context [25][26][27]. A remaining difficulty for the application of the algebraic Bethe ansatz approach to the calculation of local height probabilities arises from the expressions of matrix elements of local height operators. Even for the simple case of the cyclic SOS model with rational crossing parameter these appear to be more complicated than in the (non-dynamical) six-vertex model [26,27]. Factorization properties of the reduced density matrices in integrable face models have, to our knowledge, not been studied so far.
In this paper we address some of the issues appearing in the calculation of correlation functions in generic face models in particular on finite lattices without refering to a possible formulation as a dynamical vertex model. In the following section we specify the possible configurations of an IRF model and the Hilbert space of the related quantum chains of non-Abelian anyons. The Boltzmann weights of local configurations for these models are used to define a family of transfer matrices with generalized boundary conditions. In Section 3 the solution of the inverse problem is presented for a class of local operators for an inhomogeneous IRF model with certain local properties of the Boltzmann weights, in particular unitarity, together with an initial condition for their dependence on the spectral parameter. Expectation values of these operators in eigenstates of the transfer matrix are encoded in generalized N-point density matrices D N (λ 1 , . . . , λ N ) with independently chosen spectral parameters λ n . In Section 4 we derive a set of linear functional equations of reduced qKZ-type for D N which holds for a discrete set of values of the spectral parameter λ N .
In a final section we consider several SOS models and the related chain of Fibonacci anyons where we find that these functional equations together with the analytical properties of the density matrices derived from those of the local Boltzmann weights do in fact uniquely determine the D N in these models. Assuming that the factorization property for the N-point density matrices mentioned above also holds for integrable IRF models we propose an algorithm for the efficient computation of the structure coefficients.

Integrable face models
Interaction-round-a-face (IRF or face) models are classical statistical models defined on a square lattice with a spin (or height) a assigned to each site . The heights take values from a set S possibly subject to adjacency rules constraining their values on neighbouring vertices. These rules are conveniently presented in the form of a graph with nodes a ∈ S and adjacency matrix A ab ≡ 1, spins a and b are allowed to be adjacent 0, spins a and b are not allowed to be adjacent . (2.1) The energy of the face model for a given height configuration is determined by local Boltzmann weights depending on the spins on the vertices surrounding an elementary face [17]. These weights are allowed to depend on an arbitrary (spectral) parameter u and are depicted graphically as Related quantum models describing interacting non-Abelian anyons in one spatial dimension can be obtained from face models in their Hamiltonian limit, see e.g. [28][29][30][31]. Mathematically these anyon models can be described by braided tensor categories [32] consisting of a collection of objects {ψ a } (including an identity). They are equipped with a set of fusion rules This rule for the fusion of objects ψ a and ψ b into ψ c can be represented graphically, where the vertex (to be read from top-left to bottom-right) b a c , is allowed provided that N c ab = 0. The fusion rules allow to construct the Hilbert space of a chain of L interacting anyons with topological charge ψ a * and their possible local interactions. In this paper we consider tensor categories which are free of multiplicities, i.e. N c ab ∈ {0, 1}, and starting with an auxiliary anyon ψ a 0 an orthogonal basis of 'fusion path' states is constructed by fusing ψ a and ψ a * into ψ a +1 resulting in or, using the graphical representation of these consecutive fusion processes, Note that the sequence {a} = (a 0 a 1 . . . a L ) coincides with a possible height configuration along a horizontal line of vertices in the face model on a lattice with L × N faces provided that the possible topological charges are labelled by the elements of S and N a +1 a a * = A a a +1 .
Clearly the Hilbert space H L spanned by these fusion paths can be decomposed into sectors H L αβ labeled by the auxiliary spins α = a 0 and β = a L . Below we consider anyon chains (face models) with periodic boundary conditions (in the horizontal direction). Therefore we have to identify a 0 and a L which allows to remove the label a 0 from the basis states (2.4). As a result the model is defined on the Hilbert space ('quantum space') Note that states in H L per correspond to periodic paths of length L on the adjacency graph. Similarly, the sequences {α} = (α 0 α 1 . . . α N ) with N α n+1 α n a * = A α n α n+1 = 1 of heights on vertical lines can be identified with fusion paths spanning an 'auxiliary' Hilbert space V N of anyons. Here and below we use greek indices for the auxiliary height variables on vertical lines and latin indices to label the height variables corresponding to an anyonic quantum state on horizontal lines.
We represent the matrix elements of generic linear operators B on the space V N⊗ H L as 1 (2.6) The matrix elements of B in V N are linear operators on H L (and vice versa): (2.7) As an example, we define an operator T(u) on V 1⊗ H L from a single row of the the Boltzmann weights (2.2): 1 The symbol⊗ indicates that the index of the joint vertex of the two factors coincides.
Here, the complex numbers {u i } L i=1 parameterize local variations in the interactions around a face. Taking the trace in V 1 this gives the transfer matrix of the inhomogeneous face model with periodic boundary conditions in the horizontal direction For later use we also define the following generalized transfer operators |b = δ a n ,α δ b n ,β j =n δ a j b j = a n−1 a n+1 α = a n β = b n a|E α n 1 ...α n 2 β n 1 ...β n 2 |b = n 2 k=n 1 δ a k ,α k δ b k ,β k j/ ∈{n 1 ...n 2 } δ a j b j = a n 1 −1 a n 2 +1 α n 1 = a n 1 acting locally on a single site n or on sequences of n 2 − n 1 + 1 of neighbouring sites (subject to constraints for the heights on the neighbouring sites as a consequence of the adjacency conditions). A face model is said to be integrable if its transfer matrix commutes for different values of the spectral parameter, i.e. [t(u), t(v)] = 0. This is guaranteed by a local condition on the Boltzmann weights: the face version of the Yang-Baxter equation reads 12) or, in the graphical notation, where heights on nodes with a solid circle are summed over and heights connected by dotted lines are taken to be equal. We also assume a unitarity condition crossing symmetry 15) and the initial condition Here λ is the crossing parameter and ρ(u) a function, both are model-dependent. We assume ρ(0) 2 = 1 which can always be reached by rescaling the Boltzmann weights.

Reduced density matrices
A complete characterization of the models introduced above requires the computation of generic correlation functions, which can in turn be expressed through reduced density matrices, i.e. matrix elements of the operators (2.11) Here |Φ 0 ∈ H L per is the ground state of the model. More generally one may consider the right (left) eigenvectors |Φ ( Φ|) of the transfer matrix corresponding to a particular eigenvalue Λ(u), i.e. t(u)|Φ = Λ(u)|Φ ( Φ|t(u) = Φ|Λ(u)).
An important step towards the calculation of these quantities in integrable (vertex) models has been the solution of the 'inverse problem', i.e. expressing local operators through elements of the Yang-Baxter algebra. This has been achieved for the six-vertex and related models in [23,24,33]. In Refs. [25,26] this construction has been generalized to local spin and local height operators in face models allowing for a formulation as a dynamical vertex model [20,34,35]. Here we formulate the solution to the inverse problem for a general integrable face model, i.e. without using the existence of an R-matrix satisfying a dynamical Yang-Baxter equation, by expressing the elementary height operators (2.11) in terms of the single-row operators (2.8), (2.10) introduced above: We use the convention that empty products are 1.
Proof. Let us first proof the statement for an easy example where the chain length is L = 2, i.e. two faces per row. For |a , |b ∈ H 2 per we calculate: For general L the procedure works similarly, see Appendix A.
Unitarity of the Boltzmann weights implies This allows to reformulate the theorem as: We emphasize that this construction and proof is valid for all face models whose Boltzmann weights satisfy the unitarity conditions.
for the example N = 2 displayed above.
Comparing to (3.5) we observe that the reduced density matrix of the face model for N consecutive edges (or segments of the fusion path) in an eigenstate |Φ of the transfer matrix is obtained from D N by proper choice of the arguments λ k : In a slight misuse of notation we shall denote D N as N-site density matrix below. D N is normalized such that tr V N D N (λ n 1 , . . . λ n 2 +1 ) = 1, which gives a constraint on the diagonal elements of D N . Taking partial traces, i.e. summing over pairs (α , β ) of auxiliary indices, any n-point function with n ≤ N can be computed from D N .

Functional equations
For the next step towards the calculation of correlation functions in an IRF model the functional dependence of the density matrices D N of the generalized problem on the spectral parameters λ n 1 , . . . , λ n 2 +1 has to be found. In integrable models such correlation functions are often closely related to solutions of functional equations of quantum Knizhnik-Zamolodchikov type [36]. In the context of face models from the Andrews-Baxter-Forrester (ABF) series [17] such difference equations have been obtained by Foda et al. for the infinite lattice using corner transfer matrix and vertex operator techniques [18]. Below we show that the density matrices D N satisfy a discrete version of such equations using the local properties of the Boltzmann weights of an integrable model. Our approach resembles that of Ref. [15] for the Heisenberg model. Other than the approach of Ref. [18], it is also applicable for finite size face models.
To derive the functional equations for the density operators (3.7) we introduce the linear operator . Let B be an arbitrary operator acting on V N as defined in (2.7). The action of A N on B graphically as For models with crossing symmetry as in (2.15) the operator P − ∈ End(V 1 ) is obtained by evaluation of the Boltzmann weight (2.2) at u = λ. For more complicated cases this expression needs to be modified, see (5.16) below. Note the extra Kronecker δ's enforcing that the image of B has elements acting on H L per only. As an example consider the action of A 2 on the density matrix D 2 , here shown for a system of length L = 2: We can now formulate the main theorem of this chapter.
For general λ N the functional equation (4.2) is a difference equation for the elements of the density operator resembling the reduced quantum Knizhnik-Zamolodchikov (qKZ) equation for correlation functions of the six-vertex model in the thermodynamic limit [7,36,37]. In general Eq. (4.2) is valid only for a discrete set of values, namely λ N ∈ {u 1 , . . . , u N } -as in the finite temperature case for the spin-1/2 Heisenberg model [15]. For the face models considered here, however, it is straightforward to show that this restriction can be dropped for matrix elements of (4.2) where α N−1 is a leaf node on the adjacency graph G, i.e. if it has exactly one neighbour.
Another functional equation satisfied by the density matrices follows directly from the Yang-Baxter equation. Introducing the operator we have for 1 ≤ i < N:

Applications
In this section we will use the functional equation (4.2) to compute the density matrices of face models of solid-on-solid (SOS) type and the related anyon chains. This class of face models has been introduced by Andrews, Baxter and Forrester as an auxiliary tool to solve the 8-vertex model. Specifically, we shall consider two critical models with a finite set S of height variables, i.e. the cyclic solid-on-solid (CSOS) model [38,39] and the restricted solid-on solid (RSOS) model [17].

The cyclic solid-on-solid model
The height variables of the CSOS model take integer values 0 ≤ a ≤ r − 1 for a positive integer r. Heights on adjacent sites are required to differ by ±1 modulo r, hence a configuration of neighbouring spins 0 and r − 1 is allowed. As a consequence the adjacency graph for this model corresponds to the Dynkin diagram of the affine Lie algebraÃ r−1 , see Figure 1(a). The Boltzmann weights of the critical CSOS model are (heights in the arguments of W are taken modulo r) Here the crossing parameter is λ = πm/r where 1 ≤ m ≤ r − 1 is coprime to r (the Boltzmann weights of the general CSOS model are elliptic functions of the spectral parameter u depending explicitly on the height variable a through an additional phase angle which, however, plays no role in the critical case). Note that the weights (5.1) coincide with the nonzero vertex weights in the R-matrix of the six-vertex model. In fact, this relation has been used extensively, e.g. to identify the operator content of the low energy effective theory of the lattice model in the thermodynamic limit [40]. Furthermore, and unlike most other face models, the transfer matrix of the CSOS model has a simple eigenstate which allows for its solution by means of the algebraic Bethe ansatz method based on an R-matrix depending on a dynamical parameter related to the height variables. This property has already been used to compute form factors in the basis of Bethe eigenstates of this model [25].
Here we will utilize the existence of a particularly simple eigenstate of the CSOS transfer matrix to illustrate the approach to compute correlation functions based on the functional equation (4.2). To be specific we choose the CSOS model with r = 3 and crossing parameter λ = 2π/3. Considering a lattice of length L = 3k with integer k it is easy to verify that is an eigenstate of the transfer matrix with eigenvalue From the definition of the Boltzmann weights it follows that the action of the single row operators (2.8) on this state is This allows to analyze the density matrices D N in the state (5.2). The simplest case is N = 1 where periodic boundary conditions imply that α 0 = β 0 and α 1 = β 1 . As a consequence D 1 (λ) is a diagonal operator on V 1 whose diagonal elements can be directly read off from Eqs. (5.5), e.g.

01|D
(Ω) ) and respectively ford. Note that the trace condition . Similarly, the diagonal elements of the two-site density matrix D 2 (λ 1 , λ 2 ) in the reference state are obtained from (5.5). The functional equation (4.2) allows for the direct computation of all off-diagonal elements: choosing as a basis for the auxiliary space V 2 and using the fact that the Boltzmann weights are invariant under the shift of all heights by an integer we find that the density matrix has a structure of three identical 4 × 4 blocks D (Ω) 2 (λ 1 , λ 2 ). Restricting ourselves to the first of these blocks we find for the reference state (5.2) (5.8) or, using the notation introduced in Eq. (3.9), With this ansatz we obtain from the functional equations (4.2) and (4.4) after some algebra an explicit expression for the off-diagonal element ( Hence, as a consequence of the simple form of the reference state (5.2) of the r = 3 CSOS model, the two-site density matrix is completely determined by the one-point functionã(λ). This is also true for the three-site density matrix D for the [0, a]-blocks in the auxiliary space V 3 we find: 2 This suffices for the calculation of the nearest and next-nearest neighbour correlation functions in the reference state (5.2). In the homogeneous limit (i.e. all inhomogeneities u k = 0) we haveã(0) = 1/ √ 3,d(0) = 0, and g(0, 0) = 0. Therefore, the two and three-site density matrices for λ i = 0 are diagonal with non-zero elements 012|D (Ω) (5.12) With Eq. (3.10) this yields the expected results for the two-and three-point functions in the reference state |Ω of the r = 3 CSOS model.

The restricted solid-on-solid model
The RSOS model can be treated in a similar way. The state space is obtained by removing 0 from the set of heights allowed in the CSOS model. As a consequence the adjacency graph corresponds to the Dynkin diagram of A r−1 , Fig. 1(b). The Boltzmann weights of the critical RSOS model are again given in terms of trigonometric functions and a crossing parameter λ = π/r. They satisfy the face Yang-Baxter equation (2.12) and the unitarity condition (2.14). As a consequence of the gauge factors g x in the definition of the Boltzmann weights the crossing relation (2.15) is modified to This modification has to be taken into account whenever crossing symmetry is used, in particular in the definition of the A-operator in (4.1). To cancel the additional factors from the Boltzmann weight evaluated at u = λ we have to rescale the corresponding weight giving the operator P − ∈ End(V 1 ): In addition, the third step of the proof in Appendix B needs to be reconsidered. Keeping track of the gauge factors we find, that the A-operator needs to be multiplied by an additional factor of g β N−1 /g α N−1 .
By construction the transfer matrix (2.9) of this model and its eigenvalues Λ(u) are Fourier polynomials of degree L The leading Fourier coefficients are known to take values [41] This allows to decompose the spectrum of the RSOS model into topological sectors with 'quantum dimension' d q (j) = sin(π(2j + 1)/r) sin(π/r) , (5.19) labeled by the quantum number j. In addition there is a discrete symmetry due to the invariance of the Boltzmann weights under a reflection of all heights, i.e. a → r − a. 3 This symmetry is inherited to the transfer matrix and the reduced density operators.
Starting with the single-site density matrix D 1 (λ 1 ) we observe that only its diagonal elements are allowed to be non-zero. We will now prove that D 1 (λ) is independent of the spectral parameter λ in any eigenstate |Φ of the transfer matrix (although the matrix elements may still depend on the choice of inhomogeneities {u i } L i=1 ): to compute D [12] 1 (λ) we note that due to the adjacency condition where we have used the definition of D 1 . Note that the one-site projection operators, defined as a|P (ā) |b = δ a ā δ a j b j , |a , |b ∈ H L , (5.21) are independent of the spectral parameter. Hence, the 1-point function (5.20) is the local height probability for finding a spin a = 1 if |Φ = ( Φ|) † . With the same reasoning, one concludes that D [21] D [21] 1 (λ) + D Given that D [21] 1 (λ) is constant we find that D [23] 1 (λ) is also independent of λ. Repeating this procedure we find that in fact all matrix elements are independent of the spectral parameter and given as sums of the local height probabilities. Generically the latter depend on state |Φ and the inhomogeneities. For the critical RSOS models considered here we find, however, that they are functions only of r and the local spin in the topological sectors with quantum dimension d q = 1. Using the known values for the local height probabilities in the thermodynamic ground state of the homogeneous system [17] we find D [a,a+1] 1 (λ) = sin (πa/r) sin (π(a + 1)/r) r cos (π/r) (5.23) for the non-zero elements of the single-site density matrix in these sectors.
and also, due to reflection symmetry, 232|D 2 (λ 1 , λ 2 )|232 = 1/4. In addition, we have that Hence, we find that the non-zero blocks in D 2 (λ 1 , λ 2 ) are D [11] 2 (λ 1 , λ 2 ) = follow from height reflection a i → r − a i . Generically the two functions f and g are independent. Evaluating equation (4.4) we find that f (u, v) = f (v, u) and g(u, v) = g(v, u), i.e. the two site density operator for r = 4 is symmetric under exchange of the arguments.
Taking (5.28) as an ansatz in the functional equation (4.2) we obtain 2L linear relations for the unknown functions f and g at λ 2 ∈ {u 1 , u 2 , . . . , u L }: For the actual computation of of the density matrices we note that, as a consequence of (3.7) the elements of D N (λ 1 , . . . , λ N ) N k=1 Λ(λ k ) are Fourier polynomials in the spectral parameters λ k . We have checked that, for small N and system sizes, the (L + 1) N unknown Fourier coefficients can be determined uniquely for a given transfer matrix eigenvalue Λ(u) using the discrete functional equations (5.29) for N = 2 and similarly (4.2) for general N once D N−1 is known (cf. the appearance of D 1 in the sum rules (5.26) and (5.27) for D 2 ).
This procedure is simplified when we consider density operators for eigenstates in the sectors with quantum dimension d q (j) = 1, i.e. topological quantum numbers j ∈ {0, 1}: here we find that D 2 is determined by a single function of the spectral parameters f (λ 1 , λ 2 ) ≡ g(λ 1 , λ 2 ) such that equations (5.29) for the elements of the two-site density matrix degenerate to a set of L equations. Another simplification in these sectors is found for spectral parameter λ 2 → i∞: in this limit the functions f and g vanish and D 2 (λ 1 , λ 2 ) becomes the single-site density matrix D 1 (λ 1 ), written as an operator on V 2 using the basis (5.25). In fact, we find that a similar reduction relating D N for large λ N to D N−1 for N ≥ 2 holds in the topological sectors with quantum dimension d q = 1 of the RSOS models where (recall that D 1 is independent of the spectral parameter and diagonal): 4 Using (4.4) one obtains expressions for D N in the limit of large λ k , k < N. Hence, the asymptotics of the N-site density matrix is determined by the (N − 1)-site one, e.g. (recall that f = g in these sectors) for the three-site density matrix of the r = 4 model in the basis Remarkably, it has been shown that the density matrices of the Heisenberg spin chain can be written as in terms of a nearest neighbour two-point function ω and a set of recursively defined elementary functions f N;I,J of the spectral parameters λ j , so-called 'structure functions' [7].
Here I = (i 1 , . . . , i m ) and J = (j 1 , . . . , j m ) such that I ∩ J = ∅, 1 ≤ i p < j p ≤ N and i 1 < · · · < i m . For the density matrices in eigenstates from the topological sectors with quantum dimension d q (j) = 1 (j ∈ {0, 1} for the r = 4 RSOS model) we observe a similar behaviour, e.g. for the three-point density matrix: motivated by Eq. (5.32) we assume that the matrix elements of D 3 (λ 1 , λ 2 , λ 3 ) can be written as where f 0 and the f I,J are rational functions of e 2iλ 12 and e 2iλ 23 (λ k ≡ λ k − λ ), and f (u, v) is the single function from (5.28) which determines the two-site density matrix in these topological sectors.
Most importantly, the model parameters such as the system size L and the inhomogeneities {u k } enter the expressions (5.32) only via the two-point function ω (or f in (5.33)). This fact can be used to implement an efficient algorithm 5 for the numerical calculation of f 0 and f I,J in the ansatz (5.33) for the 3-site density matrix of the r = 4 RSOS model (or the structure functions f N;I,J appearing in an ansatz such as (5.32) for elements of the N-site density matrix D N ): 1. choose a set of spectral parameters Λ = {λ 1 , . . . , λ N }, 2. diagonalize the transfer matrix of a sufficiently small system with randomly chosen inhomogeneities, 3. pick an eigenstate of the transfer matrix (from the topological sector considered) and compute the generalized N-site density matrix D N (λ 1 , . . . , λ N ) and the two-site density matrix D 2 (λ j , λ k ) for pairs (λ j , λ k ) from Λ using their definition (3.7), 4. compare D 2 to (5.28) to obtain numerical values for the corresponding two-point functions f (λ j , λ k ), 5. insert the data from steps 3 and 4 into (5.33) (resp. (5.32)) to get a linear equation relating the structure functions, 6. repeat steps 2 to 5 to build a system of linear equations which can be solved for the structure functions f N;I,J (λ 1 , . . . , λ N ).
Once these functions are known for a range of spectral parameters it is straightforward to find analytical expressions, e.g. by Fourier analysis, which can be checked using (4.2). A slight complication in the present case of the r = 4 RSOS model is that the decomposition (5.33) is not unique. Evaluating the diagonal element α = β = (1, 2, 1, 2) of the functional equation (4.2) for D 3 (x, y, z) we find an additional relation satisfied by the twopoint function f : This identity holds for arbitrary values of λ j , j = 1, 2, 3, as a consequence of α 2 = 1 being a leaf node on the adjacency graph (c.f. the remark in Appendix B). Taking this into account we have used the procedure outlined above to compute the factorized form of the three-point density matrix D N=3 of the r = 4 RSOS model. Remarkably, it turns out to be sufficient to compute the initial data for a system of length L = 2 < N. Moreover, we find that the structure functions are the same for all eigenstates |Φ of the transfer matrix in the topological sectors considered here. As a result we obtain D [12] 3 (λ 1 , λ 2 , λ 3 ) = In addition, D 2 (0, 0) is completely fixed by the two-point function f (0, 0). For the computation of D 3 (0, 0, 0) the singularities for λ 1 = λ 2 and λ 2 = λ 3 in Eq. (5.35) has to be taken care of. We expand the two-point function as f (λ 1 , λ 2 ) (0, 0) + (1, 0)λ 1 + (0, 1)λ 2 + 1 2 (2, 0)λ 2 1 + 2 (1, 1)λ 1 λ 2 + (0, 2)λ 2 2 + . . . (5.37) with (k, ) ≡ ∂ k 1 ∂ 2 f (λ 1 , λ 2 )| λ 1 =λ 2 =0 . Note that (k, ) = ( , k) due to the symmetry of f (λ 1 , λ 2 ). Additional relations between the coefficients of the r = 4 two-point function follow from the identity (5.34), e.g. (1, 0) = 0 and (2, 0) − 2(1, 1) = 4(0, 0). As a result the singularities are removed and the homogeneous limit of D 3 is found to be As a consequence the two-and three-point correlations in a transfer matrix eigenstate |Φ from the topological sectors with d q = 1 of the homogeneous r = 4 RSOS model are given in terms of (0, 0), i.e. the numerical value of the two-point function at spectral parameters λ 1 = λ 2 = 0, alone. The latter is directly related to the corresponding eigenvalue of the RSOS hamiltonian H = J ∂ u ln t(u)| u=0 , i.e.
Hence, explicit expressions for two-and three-point functions in the ground states of the infinite system can be obtained from Eqs. (5.28) and (5.38) using the known results for the energy density of the RSOS model in the thermodynamic limit [45]. We find for the ground state of the RSOS hamiltonian with J = −1 (+1).
The r = 5 RSOS model. As a second example we consider the r = 5 RSOS model with local heights 1 ≤ a ≤ 4. Similarly as above, we can compute the single site density matrix.
In states from the sector with topological quantum number j = 0 (recall that the topological sectors in the odd r RSOS models are labelled by integers 0 ≤ j ≤ (r − 3)/2 = 1) we find that the matrix elements are independent of the system size and the inhomogeneities {u k }, namely The auxiliary space V 2 for the two-site density matrix of the r = 5 model has dimension ten. Due to the adjacency condition the Hilbert space of states splits into two spanned by fusion paths with a 0 even and odd, respectively. The transfer matrix is a map between these two subspaces. Similarly, products of an even number of transfer matrices (or more general single row operators) are therefore block diagonal and may be written as the sum of an even and and odd part. In view of this decomposition of the Hilbert space we chose the following basis for V 2 : The two sets are related by reflection a i → r − a i and hence we may restrict ourselves to the subspace generated by the first. Again the structure of the density operator D 2 (λ 1 , λ 2 ) is constrained by sum rules such as (5.26) and (5.27). We find the non-zero blocks of D 2 in the first subspace with odd a 0 to be (b is a constant) The sum rules immediately imply Furthermore the trace condition tr V 2 D 2 (λ 1 , λ 2 ) = 1 yields d(λ 1 , λ 2 ) = f (λ 1 , λ 2 ) + D [21] 1 . (5.45) Using the relations (4.4) find that the off-diagonal functions are related via c 1 (λ 1 , We find further simplifications for eigenstates of the transfer matrix belonging to the j = 0 topological sector (where b = 1/(5 + √ 5) according to Eq. (5.41)): in this sector the off-diagonal elements of D 2 coincide, i.e. g(λ 1 , λ 2 ) = g(λ 2 , λ 1 ), and therefore g(λ 1 , λ 2 ) = f (λ 1 , λ 2 ) . As a consequence D 2 can again be expressed in terms of a single scalar function f (λ 1 , λ 2 ) satisfying the functional equation for λ 2 ∈ {u 1 . . . . , u L }. As in the r = 4 RSOS model the density matrices D N can be computed recursively for any given transfer matrix eigenvalue using their analytical properties and the functional equations. In particular, we find that the asymptotic behaviour of the N-site density matrices is related to the N − 1-site one as given by (5.30), e.g. lim in the topological sector with quantum dimension d q = 1, i.e. j = 0 for the r = 5 RSOS model. Similar as in (5.35) for r = 4 we have been able to express the three-point density matrix of the r = 5 RSOS model in this topological sector as a sum of terms factorizing into spectral-parameter dependent elementary functions and the two-point function f (λ 1 , λ 2 ) solving the functional equation (5.47). Proceeding as for r = 4 we find the factorization of the D 3 in the one-dimensional block corresponding to the sequence α = (1234) of heights to be D [14] We present the complete list of non-zero matrix elements of D 3 in Appendix C. Now it is straightforward, to calculate D 3 in the homogeneous limit. Expanding the two-point function as in Eq. (5.37) for the case r = 4 we find for the one-dimensional block considered above D [14] 3 (0, 0, 0) = All other matrix elements may be computed using Appendix C.

Fibonacci anyons
As discussed earlier we can relate face models to one-dimensional quantum chains with anyonic degrees of freedom on each lattice site. Considering the Hamiltonian limit of the homogeneous RSOS model with r = 5, i.e. u i ≡ 0 in (2.8), one obtains an integrable model of su(2) 3 or Fibonacci anyons [28,46]. Despite its simplicity this non-Abelian anyon model gives rise to universal quantum computation. It contains only two types of anyons, the trivial anyon 1 and a second one, τ. Here we will use the functional equation (4.2) to compute the two-site density matrix for the chain of τ-anyons. The Hilbert space of fusion paths for these anyons can be shown to be isomorphic to the a 0 odd part of the RSOS Hilbert space H L per for r = 5. A Hamiltonian for a chain of L τanyons with local interaction can be constructed using the operators P (ττ→1) = 1 − P (ττ→τ) projecting on one of the outcomes of fusing neighbouring anyons according to the rule τ ⊗ τ = 1 ⊕ τ. In Appendix D we show that these operators can be expressed in terms of the Boltzmann weights of the RSOS model (5.13). This allows for an embedding of the anyon Hamiltonian into the family of commuting operators generated by the transfer matrix of the RSOS transfer matrix (2.9). By choosing a negative (positive) coupling constant J fusion of two neighbouring anyons to a trivial (τ) anyon is energetically favoured.
We will now use the inverse problem and relate the energies of the anyon model to the density operator of the homogeneous model. Note that in this case t(0) is the translation operator with eigenvalues Λ(0) = exp(2πik/L) for some integer k. Furthermore, we have ρ(0) = −1 for the RSOS model. Assuming that the eigenstates of the transfer matrix are normalized, Φ|Φ = 1, the two-point function (3.10) is Translation invariance, i.e [H, t(0)] = 0, implies that the energy is E Φ = L J Φ|P (ττ→1) 1 |Φ for any eigenstate |Φ . As P (ττ→1) 1 depends only on the first three heights of the chain, we can directly use (5.52) and obtain relating the energy density of the anyon chain to certain correlators of the RSOS model. Plugging in the the explicit expression of P (ττ→1) 1 into (5.53) and using the simplified form of the two-site density matrix (5.43) for eigenstates |Φ in the topological sector j = 0 of the r = 5 RSOS model 6 we finally obtain

Submission
As for the r = 4 RSOS model the ground state energies of the anyon model in the thermodynamic limit are known [45] giving for the corresponding two-point functions f (0, 0). Finally, we show how our results can be used for the computation of 3-point functions. Therefore, we consider the operator P (τττ→1) which projects the fusion product of three consecutive τ-anyons to an anyon of type 1. Using the homogeneous limit of D 3 (again for eigenstates |Φ in the topological sector j = 0 of the r = 5 RSOS model) as discussed in the previous section and (3.10) we find

Conclusion
We have studied correlation functions for generic models with interactions-round-a-face on a square lattice (and the related anyonic quantum chains). To make use of the Yang-Baxter integrability local operators have been expressed in terms of (generalized) transfer matrices for inhomogeneous versions of these models, see Eqs. (3.2) and (3.5). This allowed to encode correlation functions in N-point reduced density matrices D N , Eq. (3.7), depending on a set of N spectral parameters. We have constructed a set of discrete functional equations of reduced quantum Knizhnik-Zamolodchikov type (4.2) which determine the functional dependence of D N on these spectral parameters. This framework has been applied to several critical solid-on-solid models: in the simple 'reference' state (5.2) of the r = 3 CSOS model we have obtained explicit expressions for the generalized reduced density matrices on up to three neighbouring sites in terms of the one-point function depending on the spectral parameter and the choice of inhomogeneities. By contrast, the one-point functions in the RSOS models are independent of the spectral parameter. For the r = 4 and 5 RSOS models we have been able to express the two-site density matrices in terms of two unknown functions similar to the spin-1/2 Heisenberg chains [13,15]. We find that these functions (and, together with the application of the sum rules, the elements of the N > 2-site density matrices) are uniquely determined by the discrete functional equations for any transfer matrix eigenstate, given the analytical properties inherited from the Boltzmann weights.
Additional properties of the density matrices are found in topological sectors with quantum dimension d q = 1 (containing the ground states of the RSOS model): here we have observed a reduction relating D N to D N−1 when one of the spectral parameters is sent to infinity. This resembles the asymptotic condition on the density matrices complementing the discrete qKZ equation for the Heisenberg spin chain guaranteeing the uniquess of its solution [14,15]. Moreover, we have found that the reduced density operators of the r = 4 and 5 RSOS models in these topological sectors can be expressed in terms of a single function determining the nearest-neighbour two-site correlations in the particular transfer matrix eigenstate together with a set of elementary structure functions. This observation and preliminary results for r > 5 lead us to conjecture that this holds for all RSOS models in these topological sectors.
For the remaining tasks of calculation of the nearest-neighbour function (the physical part of the correlation functions) and the structure functions (the algebraic part) in the topological sectors with d q = 1 we have used an ansatz motivated by the 'factorized' form of density matrices for the spin-1/2 Heisenberg chain [7] together with the observation that the density operator depends on the particular realization of the model, i.e. the system size and choice of inhomogeneities, only via the two-site function. This allows for an efficient computation of the structure functions. We have applied this algorithm to reduce the calculation of the density operators on three neighbouring sites for the r = 4 and r = 5 RSOS models and related correlation functions for the Fibonacci anyon chain to that of the physical part. The latter solves discrete difference equations (5.29) and (5.47) resulting from the functional equation for the two-site density matrix. Explicit expressions for the two-site functions solving these equations are limited to RSOS models of sufficiently small length or a special state as for the CSOS model. In the fermionic basis approach for the spin-1/2 Heisenberg model the physical part of the correlation functions at both finite length and finite temperature has been characterized in terms of integral formulae and difference equations [13] or in terms of solutions to linear and non-linear integral equations [47]. For face models such representation of the two-site function is not known. As a first step, one might consider the zero temperature ground state of these models in the thermodynamic limit: choosing a continuous distribution of the inhomogeneities the functional equations (4.2) are expected to hold for arbitrary complex values of the spectral parameter λ N allowing for their solution.
To conclude let us mention just two open problems which may be addressed based on the approach presented here: first of all, in the context of RSOS models a comparison of the density matrices with the corresponding quantities for the related anisotropic Heisenberg chains at roots of unity can provide insights on the boundary contributions to correlation functions resulting from the peculiar fusion path nature of the RSOS Hilbert space. Secondly we want to emphasize that the discrete functional equations (4.2) for the density operators hold for generic integrable IRF models (such equations are also known for vertex models and spin chains related to quantum groups [48]). Together with the algorithm used here for the computation of the structure functions in factorized expressions (5.32) and (5.33) this may well allow to shed some light on the question whether the factorization of correlation functions is a general property of integrable models which extends beyond RSOS models and spin-1/2 chains.

A Proof of Theorem 1
Here we provide the proof of (3.2) for arbitrary values of L and 1 ≤ n < L. Let us first look at its graphical representation (double arrows indicating periodical boundary conditions in the horizontal direction): for |a , |b ∈ H L per we consider the matrix element After using the initial condition (2.16) in each row the Boltzmann weights can be turned into Kronecker δ's by repeated use of unitarity (2.14). To understand the principle we have

SciPost Physics
Submission a closer look at the four rows n − 1, n, n + 1 and n + 2, i.e.

α β
where we used u i,j ≡ u i − u j and ρ i,j ≡ ρ(u i − u j )ρ(u j − u i ). Continuing this procedure throughout (A.1) as on the shaded regions shown above, it is easy to see that . . . a n−1 a n = α a n+1 . . .
is produced which finishes the proof.

B Proof of Theorem 2
Here we will provide the proof for the functional equation (4.2). It consists of three steps. The idea is, to consider the action of A N on a density operator with one row added, i.e. D N+1 (λ 1 , . . . , λ N , λ N + λ). Recall, that for every n ∈ N D n (λ 1 , . . . , λ n ) {α}{β} = 0 if α 0 = β 0 or α n = β n . In those cases the functional equation holds trivially, so we may assume α 0 = β 0 . Also note that (A N (λ 1 , . . . , λ N ) [D N+1 (λ 1 , . . . , λ N , λ N + λ)]) {α}{β} is a matrix element of an operator V N+1 → V N+1 . To keep the presentation legible, the following steps will be shown graphically for N = L = 2, e.g. Section 2 for the definition of the symbol⊗) we perform two '(constrained) partial traces' over the factor V 1 each leading to one side of the functional equation, i.e. operators on V N . In a final step we show that for the special choice of λ N being one of the inhomogeneities {u i } the constraint can be dropped and both summations lead to the same result.
In a first step we note that from the definition (3.7) of the density operators (B.2) Note that this fixes the spin γ to be equal to α 1 in the graphical representation above (or α N−1 for general N). Therefore we have obtained the left-hand side of the functional equation (4.2).
For the second step we sum over α = α N = β N and fix the spin γ to be equal to β N−1 . For N = L = 2 this becomes (thick dotted lines indicate where the constraint δ γβ 1 is used) with α = (α 0 , α 1 , α 3 ) and β accordingly. Here we have used the initial and crossing conditions to evaluate the Boltzmann weight at λ, the Yang-Baxter equation to pull the rotated weight from the right to the left and finally Φ|T α 0 α 0 (λ 2 ) = Λ(λ 2 ) Φ| P α 0 α 0 with the projection operator P αα : H per → H αα . For general N > 2 these operations have to be iterated to move the row T(λ N ) to the top yielding the right-hand side of Eq. (4.2).
The final step consists of showing that for the special choice of λ N the two operations shown above yield the same result. Combining the unitarity condition (2.14) and crossing symmetry (2.15) it follows immediately that Using the initial condition for the rotated weight with spectral parameter λ we obtain that β N = β N+1 . By definition of the operator A N and periodic boundary conditions in quantum space H per we also have α N = α N+1 . The spin γ is not connected to one of the Boltzmann weights any more and therefore the partial traces considered above yield the same result. This proves the theorem. Note that the restriction of λ N ∈ {u i } in the functional equation (4.2) can be dropped for matrix elements where α N−1 is a leaf node of the adjacency graph G: in this case all neighbouring spins are necessarily equal and therefore the lowest two rows can be removed using (B.3) for arbitrary values of λ N .
Finally one should remark that, depending on the definition of the Boltzmann weigths for a particular model, the crossing relation may be modified by height dependent gauge factors, see e.g. (5.15) for the RSOS model. While these factors cancel in calculations where periodic boundary conditions can be imposed they have to be taken care of in the functional equation (4.2) -either by rescaling the operators A and D or by adding constant (i.e. not spectral parameter dependent) prefactors.

D Non-abelian anyons
The sequence of fusion processes defining the fusion path basis (2.4) of states for an anyonic quantum chain can be reordered by so-called F-moves [32]. Up to some gauge freedom the latter are determined by the fusion rules ( where e, f ∈ {1, τ} and φ ≡ 1 2 (1 + √ 5) is the golden ratio (see e.g. Ref. [28]. All others are 1 if they are allowed by the fusion rules and 0 else.
Local projection operators can be expressed in terms of the F-moves as Though depending on the sites i − 1, i and i + 1 they leave the first and the last invariant.
Since the F-moves of the Fibonacci anyons are real valued, we can drop the complex conjugation. A straight forward generalization allows to express local three anyon projection operators as (D.4) The anyons of an su(2) 3 theory can be labeled by generalized spins j = 0, 1 2 , 1, 3 2 . An automorphism of the corresponding fusion algebra allows to identify j = 0, 3 2 with the trivial (x = 1) and j = 1 2 , 1 with the τ-anyon of the Fibonacci chain. Starting from fusion path states |x 0 x 1 . . . x L of Fibonacci anyons with x n ∈ {1, τ} we obtain the Hilbert space of r = 5 RSOS model defined in Section 5.2 by mapping x n → a n ≡            1 for x n = 1, n odd 2 for x n = τ, n even 3 for x n = τ, n odd 4 for x n = 1, n even . (D.5) Note that this mapping gives only half of the basis states of the RSOS model, since a n will be even (odd) on the even (odd) sublattice. The other half is obtained by switching odd and even in (D.5). This also provides a mapping of the anyon Hamiltonian to an operator in the RSOS model. Similarly the projection operators P (ττ→1) i are mapped (up to a factor φ) to local operators e i forming a representation of the Temperley-Lieb algebra [28,50] a|e i |b = δ a i−1 a i+1 g a i g b i g a i−1 g a i+1 k =i δ a k b k (D. 6) with gauge factors g x from Eq. (5.14). Comparing this with the RSOS Boltzmann-weights (5.13) we find: We will now express the Hamiltonian by means of the transfer matrix. Therefore, we observe that where W is the derivative with respect to the spectral parameter u. Furthermore, crossing symmetry and unitarity imply that t −1 (0) = t(λ) for the homogeneous model (i.e. all inhomogeneities u i = 0). In that case, both t(0) and t(λ) are shift operators due to the initial condition. Hence, the logarithmic derivative of the transfer matrix at u = 0 is the sum of local operators is a member of the commuting family of operators generated by the transfer matrix (2.9) of the r = 5 RSOS model. Likewise we can also map the 3-anyon projector (D.4) to an operator acting on the r = 5 RSOS model Hilbert space.