Exceptional solutions to the eight-vertex model and integrability of anisotropic extensions of massive fermionic models

We consider several anisotropic extensions of the Belavin model, and show that integrability holds also for the massive case for some specific relations between the coupling constants. This is done by relating the S-matrix factorization property to the exceptional solutions of the eight-vertex model. The relation of exceptional solutions to the XXZ and six-vertex models is also shown.

The LL model is a classical integrable model and has been extensively investigated using various methods [21,22]. It is the continuous limit of the Heisenberg spin chain, which is much simpler to deal with due to absence of singular interactions. It is, however, an interesting toy model to quantize directly in the continuous limit, i.e., without invoking its discretized version [35]. This activity becomes particularly relevant, because the lattice regularization of the AAF model has not so far been constructed, and, therefore, it is an important open problem to understand the quantization of such continuous systems. Due to the singular nature of the interaction terms in the Lagrangians of both the LL and AAF models, it has been shown that a consistent quantization procedure requires the construction of self-adjoint extensions of the Hamiltonians of both theories [14,23,24,26,34]. The reason is the presence of terms proportional to derivatives of the fields in the interaction Lagrangians, which results in terms containing derivatives of delta-function in the quantum mechanical Hamiltonians. 1 In addition, the AAF model is also non-ultralocal, which further complicates the analysis of its integrable properties [7,15,33,37]. To properly address some of these problems, one has to consider the fields as operator-valued distributions, and use a suitable regularized product of quantum fields. This has been shown to produce consistent and known results for the LL model [26], and, for the AAF model, it has correctly reproduced the free fermion limit [14]. In addition, it was shown that employing distributions and a regularized field product in the quantum theory naturally reproduces Maillet symmetrization procedure in the classical limit [14]. Nonetheless, it is still very hard to proceed with concrete calculations in the case of the AAF model, 2 since the interaction terms go up to the sixth order in the fermion fields and their derivatives.
The LL model has another interesting feature which helps to understand its integrable properties without appealing to its discretized version. Namely, it is the low-energy limit of the Faddeev-Reshetikhin (FR) model [5,9,10,[38][39][40] -integrating out the high-energy modes of the FR model in the QFT description naturally leads to the derivatives of the fields in the interaction terms of the LL model [22]. The main advantage of the FR model is that the interaction term needed to extract the two-particle S-matrix does not contain any derivative of the fields, and, therefore, there are no singularities in the corresponding quantum-mechanical Hamiltonian term. Besides the method proposed by Fadeev and Reshetikhin, one can also solve the model by fermionization techiques [41,42]. In the case of the AAF model, however, there is no such simpler high-energy theory from which it could be reproduced in the low-energy limit. In fact, there are not many known purely fermionic continuous integrable systems to begin with. The most well-known Thirring model does not clearly contain enough fermionic degrees of freedom in order to integrate out the high-energy modes and still obtain as a result a non-trivial fermionic model similar to the AAF model. Thus, motivated by these considerations, the main goal of this paper is to find a new class of purely fermionic integrable models without derivatives of the fields in the interaction terms, and with enough fermionic degrees of freedom so that the resulting model in the low-energy limit is still non-trivial. In addition, we require that such a fermionic model be massive and contain interaction terms similar to those of the AAF model.
The simplest known candidate satisfying the first requirement is the model originally considered by Belavin [43] (see also [44,45]) and subsequently explored in [46][47][48][49][50][51][52][53][54][55], which exhibits asymptotic freedom and dimensional transmutation. It is a su (2) symmetric fourcomponent fermionic model which has been shown to be integrable in the zero-mass limit. Our first result is that Belavin's model is in fact integrable also in the massive case for some 1 For details see, for example, [36] 2 The explicit form of the Lagrangian for the AAF model is given in (3.1). particular cases of the coupling constants. The analysis requires several steps. First, as we discussed above, the fields should be considered as operator-valued distributions and, furthermore, the operator products should be regularized to avoid singular terms. Next, the integrability in the quantum case reduces to an analysis similar to that of the R-matrix factorization condition for the XYZ model conducted by Baxter [56]. Recall, that Baxter has deduced a general condition for such factorization to occur, and has given a general solution in terms of elliptic functions. 3 Indeed, as shown by Belavin [43] and later by Dutyshev [55], it can happen only for the massless case. There are, however, some exceptional cases for which the R-matrix factorization may also occur, corresponding to the situations in which Baxter's condition is not valid due to the degeneracy of the resulting Yang-Baxter equations. Such exceptional cases have been classified and discussed in details in [58] (see also [59,60]). For Belavin's model, we find one such exceptional solution for some particular relation between the coupling constants. 4 In addition, we generalize our analysis by adding two new interaction terms that mimic the interaction terms of the AAF model. The resulting S-matrices are shown not to have the form of the R-matrix in the XYZ model, but rather its generalized inhomogeneous form, which had been considered in [58], and where its various non-trivial solutions have been classified, including the exceptional cases. Here too, we find that there are such exceptional solutions, corresponding to the massive theories.
Our paper is organized as follows. In section 2, we consider the anisotropic Belavin model, introduce Sklyanin's product to deal with singular operator products, and obtain the conditions necessary for the diagonalization of the two-particle sector. In section 3, we extend these considerations for the Lagrangian with new interaction terms resembling those of the AAF model. In section 4, we consider the necessary conditions for the integrability of the model, i.e., the S-matrix factorization condition, and show that it reduces to finding the exceptional solutions of the eight-vertex model with homogeneous and inhomogeneous S-matrix structures. In section 5, we show the relation of our exceptional solutions to the XXZ and six-vertex models. In the last section 6, we outline some possible continuations and open problems.

Anisotropic Belavin model
We start from the classical Lagrangian density introduced by Belavin in [43] and generalized to the anisotropic case in [55] to describe the interacting theory of an isotopic fermion in (1 + 1)-dimensions ψ α i (x, t) with spinorial index i ∈ {1, 2} and isotopic index α ∈ {1, 2}. The isotropic interaction coupling constant g 0 is supplemented by three distinct coupling constants g a , a ∈ {1, 2, 3} for the isotopic interaction. The γand τmatrices are chosen as in terms of the usual representation for the Pauli matrices: The spacetime metric is g µν = diag(1, −1) and the metric in the isotopic space is simply δ ab . The Hamiltonian density corresponding to (2.1) is where we introduced the following matrices for the interaction terms. The Poisson algebra reads:

Quantization and Sklyanin product
To properly quantize and diagonalize a continuous hamiltonian such as (2.5) without resorting to a lattice regularization, it is necessary to deal with the singularities associated with operator products at the same point. As discussed in [14], this requires a two-step regularization procedure involving both the regularization of the fields and the operator product.
First, one has to take care of the singularities related to the canonical anticommutation relations: (2.9) in the limit x → y. One solution to this problem is to formulate the quantum theory in terms of operator valued distributions (2.10) in which the function F µ (x, y) is a symmetric test function which satisfies The c-number ζ can be subsequently fixed to where ζ 0 ∈ R, by requiring the algebra of regularized fields to be finite in the limit x → y. In (2.12), the function δ µ corresponds to the following regularization of the delta function with respect to the regularizing parameter L = 1 /µ associated to the length of the box in which we consider our system, resulting, as a consequence, in asymptotic Bethe Ansatz equations.
It is important to emphasize that treating the fields as operator valued distributions has fundamental consequences in the context of integrable models. For instance, it is in general only possible to obtain a meaningful Yang-Baxter equation from which one can subsequently derive the trace identities in terms of fields regularized as (2.10). One notable example of this situation is given by the LL model, in which the afore mentioned procedure allows to avoid ill-defined expressions such as ∂ 2 x δ(0). Furthermore, only by treating the fields as operator valued distributions one can construct the appropriate self-adjoint extensions and derive the corresponding spectra [26].
There remains, however, the singularities arising from possible discontinuities of the wave-function and its derivatives during the diagonalization of the quantum Hamiltonian. This issue can be sorted out by also regularizing the operator product of fields. Let A i , i ∈ Z + denote an arbitrary set of fields and consider the operator corresponding to (2.10), i.e., (2.14) We define the •-product as [14]: 5 Here the integration is taken over a k-dimensional hypercube of side ∆ ∈ R around the point x where all possible singular (k − 1)-dimensional hyperplanes are removed from the integration domain. Hence, there are k! disjoint regions of volume ∆V i , 1 ≤ i ≤ k! corresponding to all possible orderings of the variables ζ j , 1 ≤ j ≤ k separated by the length of a regularization parameter ǫ ∈ R. As argued in [14], the •-product (2.15) correctly reproduces the Maillet symmetrization procedure for non-ultralocal models [7] in the classical limit. To further clarify the meaning of the •-product (2.15) and prepare for the diagonalization of the quantum hamiltonian, in the following, we evaluate in details the •-product of two and four F -regularized fields (2.14). The •-product of two operator valued distributions corresponds to: b . 5 The •-product is a modified Sklyanin product of [35].
In (2.16), the integration is taken over a square of side ∆ minus a strip of size 2ǫ around the diagonal ξ 1 = ξ 2 , which corresponds to the union of the triangles above the line ξ 1 = ξ 2 − ǫ and bellow the line To compute the first integral in (2.16), we first note that the integrand is a product of continuous functions. Thus, we can perform the integrations over ξ 1 and ξ 2 by simply invoking the mean value theorem to obtain Here, c 1 ∈ (ǫ, ∆) and c 2 ∈ − ∆ 2 , ∆ 2 − ǫ . Before proceeding with the computation of P (2) a , we verify that the line z 1 = z 2 is indeed excluded from the integration domain. To reach this conclusion, we first note that we can choose the support of the F -functions to be of the form: (2.18) in which µ ∈ R act as a regularization parameter. Hence, so that the condition for disjoint supports corresponds to: Finally, since c 1 < ǫ, it is sufficient to take µ 1 , µ 2 < ǫ 2 to completely remove the line z 1 = z 2 from the integration domain.
We can now take the limit ∆ → 3ǫ to obtain: Finally, collecting the terms corresponding to all permutations over the arguments of the F -regularized fields, we obtain: (2.26)

Quantum Hamiltonian
From the expression (2.15) and the examples involving two and four regularized fields, it is clear that the •-product essentially "smears" the product of operators around a volume of side ∆ avoiding the singularities at coinciding points. Thus, by combining the field regularization with the operator product regularization, we can formulate a quantum theory with a well-defined algebra of operators and construct the corresponding Hilbert space by diagonalizing the quantum Hamiltonian. Following this prescription, we propose as the quantum Hamiltonian corresponding to the anisotropic Belavin model (2.5).
Next, we will construct the Hilbert space of the theory by considering the action of (2.27) on the N-particle states defined by: with an anti-symmetric wave-function: Here, |0 denotes the pseudovacuum, i.e., the state annihilated by the field operator The vacuum energy is, therefore, given by H|0 = 0. Before proceeding with the diagonalization of the 1-and 2-particle sectors, it is important to emphasize that the wave-function φ i 1 ···i N α 1 ···α N (x 1 , . . . , x N ) is in general discontinuous and, therefore, not defined on the diagonal x 1 = x 2 = · · · = x N . Hence, the integral in (2.29) should be understood in the sense of removing these points of measure zero. Thus, for example, for the case N = 2, the integral in (2.29) should be understood as follows: x 2ˆ∞ This leaves open the question of whether one needs to add a compensating "localized" state at x 1 = x 2 . This is indeed the case for the AAF model, as explained in [34]. Such localized state was shown to be necessary in order to diagonalize the Hamiltonian of the AAF model, and was due to the fact that the corresponding wave-function was discontinuous together with its derivatives. Here, however, there is no need to add such a localized term.

1-Particle Sector
To evaluate the action of the regularized Hamiltonian (2.27) on the 1-particle state we use the regularized algebra to normal order each term in (2.33).
We start with the kinetic term In order to reduce it to the form (2.16), we use the anti-symmetry of the derivatives of the F -functions together with their rapidly decreasing behavior to move the derivative from the F -functions to the wave-function through a series of integrations by parts: Thus, we can integrate over u and v following the steps outlined in the section 2.1 arriving at: There remains to remove the regularization parameters µ and ν. Since µ, ν < ǫ 2 , we can simply take the limits µ, ν → 0 using (2.11) to obtain: The mass term in (2.33) is already in the form (2.16), so we can simply use the formula (2.22) to get: Then, removing the regularization parameters µ, ν, we can perform the integration over y and ξ: Finally, noting that the interaction term in (2.33) is of the fourth order in the fields, we can easily conclude from the algebra (2.34) that it trivially vanishes. Hence, we are left we the following expression for the action of the hamiltonian (2.27) on the 1-particle state: The wave functions describing the 1-particle state (2.29) can be derived by solving the differential equation: following from (2.42). The solution corresponding to the positive mass shell, is parametrized by a real rapidity θ. Here, k denotes the wave vector and u(θ) is the Dirac spinor: (2.45) On the other hand, the negative mass shell can be trivially obtained from (2.45) after performing the transformation: θ → iπ − θ.

2-Particle Sector
Next, we consider the action of the regularized Hamiltonian (2.27) on the 2-particle state. Using the property (2.36) and integrating by parts to remove all derivatives acting on Ffunctions, we obtain: Here, correspond to the boundary terms arising from the integrations by parts, while the kinetic, mass and interaction terms are given respectively by The derivation of (2.46) involves integrating by parts with respect to x 1 . Since before diagonalizing the Hamiltonian, we cannot know the exact dependence of the wave function or its derivatives on their arguments, we must ensure that all possible singularities arising from the line x 1 = x 2 be taken into consideration accordingly. Therefore, all integrals with respect to x 1 and x 2 should be carefully evaluated on the line x 1 = x 2 . In the following, we will discuss in more details how this can be consistently achieved as we evaluate the contributions from (2.48) and (2.49). We start by considering the simpler boundary term B (2) ξ . In this case, we can trivially exchange the integration over u 1 and u 2 with the derivative with respect to ξ to obtain an expression of the form (2.16). So that, we can follow the steps outlined in the section 2.1 to conclude that this boundary term is proportional to: For fixed x 1 = z or y = z and sufficiently small µ, the support condition of the F -functions (2.18) implies that a / ∈ supp F µ (z, a) in the limit a → ∞. Hence, this boundary term must trivially vanish, i.e., The second boundary term, B x , requires a more careful consideration involving the explicit use of (2.32) because of the possible discontinuities of the wave-function and its derivatives. Taking that into account, we can proceed as before and exchange the derivative with respect to x 1 with the integration over u 1 and u 2 to obtain: After integrating the total derivative and removing the regularization parameters µ, ν < ǫ 2 , we obtain: We can now shift the integration over y, so that it is possible to use the continuity of the fields on their arguments to remove the cut off η. The resulting contribution from the boundary terms amount to: The kinetic (2.49) and mass (2.50) terms require a similar careful analysis regarding the integration over x 1 and x 2 . Following exactly the same steps as before, we obtain: for the kinetic term and for the mass term. Finally, the contribution of the interaction term (2.51) can be evaluated by following the prescription outlined in the section 2.1. In this case, the interplay between the •-product and F -regularization guarantees that the singular hyperplane x 1 = x 2 is never reached without any further consideration. More precisely, the support condition of the F -functions (2.18) together with the requirement that µ 1 , µ 2 , ν 3 , ν 4 < ǫ 2 naturally exclude this hyperplane from the domain of integration. Thus, where the permutation is taken over u 1 , u 2 , u 3 and u 4 which depend both on x and on the regularizing parameter ǫ as: To simplify the sum over all permutations appearing in (2.59), we shifted the integration variables so that all terms had wave-functions with symmetrical arguments. Then, we rescaled the regularizing parameter for each term individually and used the continuity of the fields on their arguments to partially evaluate the limit ǫ → 0, leading to: Here, we introduced the following shorthand notation for the relevant combination of wavefunctions with respect to the dependence on the regularizing parameter: It follows from the anti-symmetry of the wave function and the continuity of the fields that the afore defined quantity ∆ (±) ij αβ satisfies the following property: In the remainder of this section, we consider the necessary conditions for diagonalizing the 2-particle sector. From the vanishing of the second line in (2.62), we obtain the following equations describing the discontinuities of the 2-particle wave function: where for further convenience we introducedξ = ζ 2 4 . All the remaining combinations ∆ (±) ij αβ are continuous by virtue of (2.64). To solve the system of equations (2.65) -(2.68), we invoke the Bethe Ansatz to write the 2-particle wave functions as: Here k 1 and k 2 are wave vectors, A ij αβ are the unknown A-amplitudes and u(θ) is the Dirac spinor (2.45). Substituting this Ansatz in the system (2.65) -(2.68), we obtain the following system of equations: where θ ij = 1 2 θ i − θ j , for the B-amplitudes defined as: For the sake of clarity, we introduced the following combinations of the coupling constants: (2.79) Note that the λ µ constants will enter into expressions for physicial amplitudes, and in order to obtain finite expressions for the S-matrices one has to appropriately renormalize the "bare" coupling constants g 0 , . . . , g 3 . This will be considered in the subsequent section. Thus, the B-amplitudes are related as with the γ µ (θ), 0 ≤ µ ≤ 3 coefficients given by: The n ≥ 3 case will be considered in section 4 in relation to the S-matrix factorization property, implying the quantum integrability of the model. First, however, we briefly explain in the next section how to generalize the results from the previous analysis by introducing additional interaction terms, the form of which are stipulated by the corresponding interaction terms of the AAF model. Although this was our main motivation, one can in principle easily modify the analysis below by considering different, more general interaction terms.

Extended anisotropic interactions
As discussed in the introduction, our main motivation is to find an integrable massive fermionic model which contains no derivatives of the fields, and which may produce the AAF model in the low-energy limit, possibly after a suitable redefinition of the fields. This stems from the observation that the interaction term of the corresponding Hamiltonian has the form of a decomposition in 1 /m and goes up to the third order, resembling a low-energy expansion of a more fundamental theory. 6 Recall, that the Lagrangian of the AAF model has the form: We note here, that, as was shown in [32], the S-matrix factorization, implying the quantum integrability of the model, occurs only when the condition (g 2 ) 2 = g 3 is satisfied. As a first step in this direction, we extend in this section the Lagrangian (2.1) to include two new interaction terms, with the coupling constants g 4 and g 5 , that resemble that of the AAF model. Thus, we consider the following model: 7 (3.2) 6 We do not reproduce here the explicit lenghty expression for the Hamiltonian, which can be found in [14]. 7 We emphasize that although we use the same notation for the fermionic field in the Lagrangians (3.1) and (3.2), ψ in the former Lagrangian is a two-component spinor, while in the latter Lagrangian it has four independent components.
Note that, without considering anisotropic extensions of the Thirring model, it would have been impossible to write a term proportional to ε αβ . We stress that we do not try here to find a Lagrangian that would exactly reproduce the AAF Lagrangian (3.1), but merely to consider the consequences on the intregrability by adding some terms that resemble the interaction terms in (3.1). Namely, in this case, the two new terms in (3.2), anisotropic with respect to τ 3 matrix, were added to resemble those quartic in the fermionic field interaction terms of the AAF Lagrangian (3.1). These terms indeed satisfy our requirements of absence of the derivatives in the interaction terms. The natural question is then how such derivative terms may be obtained from (3.2), and reproduce those of (3.1). First, it is easy to see that such derivative terms will appear upon choosing two fields to be integrated out. Another possible way to relate the two Lagrangians is to consider a field transformation with the generic form: ψ −→ ψ + γ α ∂ α ψ + · · · . It is clear that amongst the resulting terms after such transformation one will find exactly a term corresponding to the quartic terms in (3.1). Using the so-called equivalence theorem (see [34] for details), one can show that the S-matrix is unchanged under such field transformations, and moreover, the quartic terms are essentially the ones determining the complete S-matrix for an integrable system. These general arguments will be considered in details elsewhere.
The analysis of the previous section concerning the 1-and 2-particle sectors of the anisotropic Belavin model can be easily generalized. The original tensors G αβ,α ′ β ′ and Λ jl,km (c.f. (2.6) and (2.7)) should now be split into two pairs to accommodate the new interaction terms, the first, corresponding to terms proportional to g 0 , . . . , g 4 , is (3.4) and the second, accounting for the terms proportional to g 5 , is Then, the system of equations describing the discontinuities of the wave-function φ(x 1 , x 2 ) (2.29) extends the original system (2.65)-(2.68) as follows: in terms of the same shorthand combination of wave-functions (2.63). Again, all the remaining combinations ∆ (±) ij αβ are continuous by virtue of (2.64). To solve this system of equations, we first substitute the same wave function for the 2-particle sector obtained from the Bethe Ansatz for the original Belavin model (2.69) and (2.70). The resulting system of equations for the B-amplitudes (2.75) reads: Here, we also use the same combinations of coupling constants introduced in (2.76) -(2.79) and similarly define λ 4 = 2ξ g 4 and λ 5 = 2ξ g 5 . Thus, the B-amplitudes are related as: where the coefficients of the mixing matrix are given by: Note that, differently from the original Belavin model, in which the mixing matrix was diagonal (c.f. (2.80)), in the extende model, there is some mixing between different amplitudes due to the new interaction terms. We conclude this section by giving the explicit form of the Lagrangian (3.2) written in terms of the fields convenient for the purposes of considering the low-energy limit: The Lagrangian (3.2) in terms of these fields becomes: where the interaction term has the form:

S-matrix factorization, Baxter's eight-vertex model and its exceptional solutiions
In this section, we analyze the problem of S-matrix factorization for n ≥ 3, and its relation to the (generalized) eight-vertex model [56] and in particular its exceptional solutions [58]. We fix our notations following [43,55,57]. Let I = (1, 2, . . . , N) with x 1 < x 2 < . . . < x N denote the fundamental ordering, and P = (P 1 , P 2 , . . . , P N ) with x P 1 < x P 2 < . . . , < x P N be an arbitrary ordering sector. Then the Bethe Ansatz for the wave-function in the P-sector, denoted by φ P α q 1 ···α q N (x), is of the form: where PQ stands for the product of permutations, k i = m sinh(θ i ), and u j (θ) = 1 √ 2 cosh(θ) e − θ /2 e θ /2 is the spinor for the jth particle. This form of the wave-function satisfies the anti-symmetry condition (2.30). The quantum integrability of the model, expressed as the factorization of the N-particle S-matrix in terms of two-particle S-matrices, follows from (4.1), provided the Yang-Baxter equation is satisfied. This can be easily seen by recasting the above expression into the following equivalent form [57]: where we have combined the spinor and isotropic indices into σ i = {m i , α i }, and (σ, x) stands for (σ 1 , x 1 ; . . . , σ N , x N ). 8 The A-coefficients relate to the two-particle S-matrix via the exchange relation: The S-matrix should satisfy the standard normalization and unitarity conditions. The factorization property of the N-particle scattering matrix readily follows, and the consistency for N ≥ 3 requires the Yang-Baxter equation (YBE): (4.5) to be satisfied. For all the cases considered in this paper the S-matrix elements depend only on the differences of rapidities: θ ij := (θ i −θ j ) /2. In this case the YBE equation has the form: For the original model the S-matrix has the form: Thus, it has the form of the R-matrix for the eight-vertex model [56]. The YBE was investigated in [43,55] for the original model and it was shown that it corresponds to Baxter's general solution in terms of the elliptic functions [56] only in the massless case: m = 0. In the next section, we analyze this point more carefully, and show that there exists an exceptional solution to YBE for the massive case. Then, in the subsequent section we extend this result for the model including the interaction terms with g 4 and g 5 coupling constants, in which case the S-matrix is not of the type as in (4.7), but has a more general (inhomogeneous) form. The classification of R-matrices with inhomogeneous structure, together with the exceptional solutions to YBE, was considered in [58,59]. We show that the integrability of the model with g 4 and g 5 coupling constants corresponds exactly to such exceptional solutions.

Baxter's general and exceptional solutions for homogeneous S-matrices
It was shown by Baxter [56] that the most general solution given in terms of elliptic functions for the R-matrix of the eight-vertex model, which has the same form as the Smatrix in (4.7), requires the following ratios to be constants: (4.8) 8 The equivalence between (4.1) and (4.2) easily follows from identification: In this case, for the original model, the coefficients (a, b, c, d) in (4.7) are given by the following formulas: with γ µ (θ), 0 ≤ µ ≤ 3 given by (2.81).
The key point is that Baxter's condition (4.8) and consequently the general solution in terms of elliptic functions are valid provided all six Yang-Baxter equations are independent [56,57]. This indeed happens, as was shown in [43,55], only in the massless (m = 0) case. As a consequence, this requires taking the limit θ → ∞, in order to keep the momentum p = m sinh(θ) finite. In such a limit, the formulas (4.9)-(4.12) simplify considerably together with the six independent Yang-Baxter equations, and the massless case can be reduced to Baxter's general solution for the eight-vertex model.
We emphasize that, in general, verifying only the conditions (4.8) is not enough. Namely, one should also verify that the six Yang-Baxter equations are independent. This is indeed the case for the coefficients (a, b, c, d) in (4.9)-(4.12). Even though, the parameters Γ and ∆ defined in (4.8) are constants for some relations between λ 0 , . . . , λ 3 , they do not provide a solution for the Yang-Baxter equations, since the YBE system fail to be independent. There are, however, some exceptional solutions to the YBE which correspond to two possibilities: either some of the (a, b, c, d) coefficients in (4.7) coincide, or some of the (a, b, c, d) coefficients in (4.7) are equal to zero. For a detailed analysis and a subsequent classification of such exceptional cases, we refer to the original paper [58]. For the massive case (m = 0), we find exactly such exceptional solutions, which we list below. First, we address the finiteness of the S-matrix and renormalize the coupling constants (see the comment after (2.79)), and recall the relation between the original coupling constants g 0 , . . . , g 3 and λ µ = 2ξg µ ; µ = 0, . . . , 3. In order to write the S-matrix in terms of finite quantities, one should start with the Lagrangian written in terms of the renormalized coupling constants: (4.13) Below we list the (real) exceptional solutions in terms of both λ µ andg µ , and give, for each solution, the corresponding interaction Lagrangian, which can be readily found from the explicit form of the general Lagrangian given in (3.24).
In this case, the S-matrix has the form: The resulting interaction Lagrangian has the form: In this case, the S-matrix has the form: (4.16) The resulting interaction Lagrangian has the form: (4.17) In this case, the S-matrix has a same form similar to that of the Solution 2: (4.18) The resulting interaction Lagrangian has the form: In this case, the S-matrix can be written in the form: 20) where P is the permutation matrix: It is interesting to note that this solution imposes the following bound on the coupling constants:g 1,2,3 ≥ 1 /2. There are two resulting interaction Lagrangians corresponding to (±) solutions above: (4.23) and (4.24) Note that the square root structure is somewhat similar to that of the AAF model (3.1), for which, as we commented in section 3, S-matrix factorization requires the condition (g 2 ) 2 = g 3 for the coupling constants in (3.1).
In this case, the S-matrix has the form: The resulting interaction Lagrangian has the form: In the next section we extend these considerations for the interaction terms with g 4 and g 5 coupling constants.

Exceptional YBE solutions for inhomogeneous S-matrices
The inclusion of the interaction terms proportional to the coupling contants g 4 and g 5 results in the following (inhomogeneous) S-matrix: where The coefficients (a i , b i , c i , d i ), i = 1, 2 for the S-matrix can be trivially obtained from the coefficients of the mixing matrix (α i , β i , γ i , δ i ), i = 1, 2 (3.16) -(3.21) by simply changing from the B-amplitudes back to the original A-amplitudes using (2.75). Thus, we have a 1 = a 2 and b 1 = b 2 . In such a case, there are in general twelve independent Yang-Baxter equations following from (4.6) (for a detailed exposition see [58]): Here θ stands for θ 12 and θ ′ stands for θ 23 (c.f. (4.6)). The parameter η is defined via the relation d 2 = ηd 1 , and is required to be a constant for consistency of the Yang-Baxter equations. A necessary condition for a general solution of the inhomogeneous YBE to exist is [58]: However, for the extended Belavin model, the b-coefficients satisfy: Hence, only for λ 5 = 0 the condition (4.35) holds. Fixing g 5 = 0, there exists only one general real solution with non-zero g 4 , corresponding tog 0 =g 1 =g 2 =g 3 = 0,g 4 = 0. Nonetheless, this solution is rather trivial, since it corresponds to a pair of decoupled Thirring models (3.24). Besides this real solution, we note, for completeness sake, that it is possible to find some quite non-trivial solutions involving complex coupling constants, the physical meaning of which is currently not clear to us. Nevertheless, since the general analysis conducted in [58] relied on non-vanishing and non-coinciding S-matrix elements, it is still possible to find some exceptional solutions to the YBE with g 4 = 0 and g 5 = 0 by violating any of these hypotheses. One such solution corresponds to: 37) or, in terms of the renormalized coupling constants: so that the coupling constantsg 4 andg 5 are left unconstrained. The resulting S-matrix is: Note that in the limitg 4 → 0, the above S-matrix reduces to that of solution 5 (4.25). The interaction Lagrangian corresponding to this case has the following form: (4.40)

Relating XXZ and six-vertex models to exceptional solutions
In this section, we consider in details some of the solutions previously derived in section 4.1 and match their S-matrices to known models. To do so, we will first factorize each Smatrix as so that the c-coefficients ofS (i) are normalized to one. Thereby, we will be able to identify the transformation of variables which relate the models described by the S-matrices of section 4.1 to the XXZ and six-vertex models. For the sake of completeness, before proceeding with this analysis, we recall the form of the XYZ Hamiltonian on a closed chain of N sites where the spin-1 /2 operators are represented by Pauli spin operators.

Solution 1
The S-matrix for the six-vertex model on a square lattice can be parametrized in terms of hyperbolic functions as [56]: with respect to the complex parameters u, λ, ρ. It can then be considered an entire function of u, while the remaining variables λ and ρ are treated as constants, the latter being just some normalization factor usually taken as unit. The coupling constants J x , J y , J z of the Hamiltonian (5.2) describing the corresponding XXZ model can be obtained from the following relations (see, for example, [57]): To match the S-matrix (4.14) describing solution 1: (g 1 = 0,g 2 = 0,g 3 = −g 0 ), we introduce, according to (5.1), h (1) (θ,g 0 ) = 4g 0 4g 0 cosh 2θ + i 1 − 4g 2 0 sinh 2θ . (5.5) In this case, the resulting S-matrix can be reduced to Baxter's six-vertex S matrix (5.3) provided we take the rapidity θ = u /2, fix ρ = 1 and identify the parameter λ as: with |g 0 | ≥ 1 /2. Therefore, according to (5.4), solution 1 corresponds to a XXZ model described by Thus, up to an overall factor corresponding to the function h (1) (θ,g 0 ) given by (5.5), the S-matrix (4.14) amounts to the XXZ model with interaction constants given by (5.8).

Solutions 2 and 3
The analysis corresponding to solutions 2 and 3 requires that we consider the more general S-matrix describing the 8-vertex model on a square lattice. We recall that it can be parametrized in terms of elliptic functions as [56][57][58]: Here, sn(u, k) is the Jacobi elliptic function of argument u and modulus k, η and γ are complex parameters. The coupling constants of the Hamiltonian describing the corresponding XYZ model can be obtained from the following relations [57]: To match the S-matrix of the XYZ model (5.9) with the S-matrix (4.16) describing solution 2: (g 1 = −g 0 ,g 2 = 0,g 3 = 0), we introduce, according to (5.1), the function: The corresponding S-matrix becomes: We can, therefore, relate the S-matrix of the XYZ model (5.9) and the above matrices S (2) andS (3) as follows. Setting the parameter γ = 0 in (5.9), identifying the variable u in (5.9) with the rapidity θ, i.e., making u = θ, and taking the limit of the modulus k → 1 one reproduces exactly the matrix in (5.12) starting from the S-matrix of the XYZ model, provided the following relation for the η parameter: is satisfied. To obtain the matrix in (5.14), one proceeds exactly as before but with γ = 2πi. It then follows from (5.10) that the J x , J y , J z interaction constants in the Hamiltonian for the XYZ model take the form: J x = 1 − 4g 2 0 and J y = J z = 1 + 4g 2 0 . (5.16) Thus, up to an overall factor, the function h (2) (θ,g 0 ) in (5.11) of rapidities and the coupling constantg 0 , the S-matrices, S (2) (4.16) and S (3) (4.18), correspond to the XXZ model, with the interaction constants given by (5.16) depending on the coupling constantg 0 of the corresponding fermionic models (4.17) and (4.19).

Solution 4
The S-matrix (4.20) corresponding to solution 4:g 1,2,3 = 1 3 −g 0 ± 3 + 4g 2 0 can be easily mapped to one of the exceptional solutions considered by [58]: in (5.17) to obtain the S-matrix (4.20). This exceptional solution to the YBE corresponding to the case b 1 = b 2 = 0 can be reduced to the usual S-matrix of the XXZ model.

Conclusion
Motivated by finding a simpler theory with interaction terms similar to that of the AAF model, we considered in this work the integrable properties of the Belavin model, and some of its anisotropic extensions. The AAF model is a quite complex fermionic model which is hard to investigate using the standard methods in the context of integrable system, due to its highly non-ultalocal nature and singular potentials. On the other hand, the form of the AAF action does resemble the typical 1 /m expansion arising in the low-energy limit of some more fundamental theory. Thus, our motivation was to find a massive two-dimensional fermionic integrable model which had enough fermionic degrees of freedom to perform such a lowenergy 1 /m expansion (see for example [61] for recent methods).
We have found in this paper, for the simplest su(2) invariant Belavin model, that the integrability in the massive case requires the investigation of exceptional solutions to the eight-vertex model. Furthermore, we showed that under some conditions on the coupling constants such exceptional solutions indeed exist. This is to contrast with the massless case, which had been extensively investigated, for which one can write Baxter's general solution for the eight-vertex model. It is easy to perform the low-energy expansion, at least to the lowest order, for one of our integrable solutions, e.g., for the interaction Lagrangian given in (4.23). The resulting two-component massive fermionic model indeed has a form similar to that of the AAF model with the characteristic 1 /m expansion (the explicit lengthy expression for the AAF Hamiltonian can be found in [14]). Although in this case we do not exactly reproduce the terms of the desired form, the resulting low-energy model does exhibits some similar features. Namely, as discussed in the introduction section, the interaction terms contain the derivatives of the fields, which would make the investigation of this model quite a difficult task, had we not had known that it is the low-energy expansion of a much simpler integrable model. There is of course still much work to do in order to reproduce exactly the terms of the AAF model, and this problem remains open.
There are several ways to extend our results. Firstly, one can add more general interaction terms and investigate the existence of exceptional solutions to the YBE corresponding to the inhomogeneous form of the S-matrix. Secondly, it is straightforward to generalize our results to su(n) case, and it would be interesting to find all non-trivial integrable models for the massive case, which would involve the exceptional solutions. Finally, we mention the following open problem, which we leave to a future publication. It was shown previously in [14,15] that the AAF and the free fermion models are non-ultralocal. Since there is so far no lattice version of the AAF model, it is an interesting problem to first relate the lattice formulation of the free fermion model [62][63][64][65] to the continuous limit, which results in a non-trivial Lax pair and a non-ultralocal integrable structure.