Classical r -matrices, “elliptic” BCS and Gaudin-type hamiltonians and spectral problem

In this article we make a brief review of the generalized Gaudin spin chains with [21] and without [18], [19] external magnetic ﬁeld in the simplest case of g = so( 3 ) and provide the r -matrix interpretation and generalization of the recent results of [38], [39], [40], on their exact solution. We present the classiﬁcation theorem for the r -matrices that are diagonal in the natural basis and for the “shift elements” playing the role of the interaction strengths and external magnetic ﬁelds in the corresponding Gaudin-type models. We apply these results to the description, classiﬁcation and solution of BCS-Richardson-type models. In particular, we introduce “elliptic” BCS-Richardson hamiltonian corresponding to non-skew-symmetric elliptic r -matrix and calculate its spectrum. We show that the “closed” [23], [30] and the “open” [36] BCS-Richardson hamiltonian of p + ip -type coincide with its trigonometric degenerations.


Introduction
The Gaudin models [1] are quantum integrable models describing spin chain of N spins with a long-range interaction among the spins in the chain that are characterized by quadratic spin hamiltonians in which the interaction strengths r αβ (ν k , ν l ) of the componentwise spin-spin interaction are the matrix elements of the classical r-matrix: r(ν k , ν l ) = 3 α,β=1 r αβ (ν k , ν l )X α ⊗ X β , satisfying the classical Yang-Baxter equation [3] and skew-symmetry property r αβ (ν k , ν l ) = −r αβ (ν l , ν k ). The parameters ν k , k ∈ 1, N may be interpreted as the coordinates of the spins of the chain, {X α , α ∈ 1, 3} is the basis of so (3).
A slight modification of the Gaudin model -the so-called Gaudin models in an external magnetic field are characterized by linear-quadratic hamiltonians [7] where the coefficients of the linear part of the spin hamiltonians are characterized by an element of so(3) c = 3 α=1 c α X α belonging to the symmetry algebra of the classical r-matrix.
The revival of the interest to the Gaudin models has begun after the paper [11], where it was shown that the famous Richardson model [2] of nuclear physics is a partial case of the Gaudin model in the external magnetic field corresponding to the standard rational r-matrix. Since that there have been a number of papers (see e.g. [13], [14], [15]) proposing to construct the Richardson-type models with the help of some other skew-symmetric, in particular, trigonometric solutions of the usual classical Yang-Baxter equation. The review of this approach may be found in [16]. Some applications of the algebra of Lax operators governed by skew-symmetric r-matrices were considered also in [17].
In our papers [18], [19], [21] we have constructed classical and quantum generalizations of Gaudin models with and without external magnetic field associated with general non-skewsymmetric classical r-matrices satisfying, instead of the usual classical Yang-Baxter equation, the so-called "permuted" or "generalized" classical Yang-Baxter equation [10]. The generalized classical Yang-Baxter equation possesses much more solutions than the usual classical Yang-Baxter equations does and includes all non-degenerated skew-symmetric solutions of the classical Yang-Baxter equations as a subset of its solutions. That is why the variety of the generalized Gaudin models with or without external magnetic field is much wider than the variety of the usual Gaudin models with or without external magnetic field. Observe that the role of the external magnetic field is played by the so-called "shift elements" [21], [25] that are not in general connected with the symmetries of the r-matrix. In the subsequent papers [22], [23], [24] we have proposed to construct Richardson-type models with the help of our generalized Gaudin models in an external magnetic field. In such a way we have constructed, in particular, the so-called p + ip BCS-Richardson model [23], [24]. We remark, that the p + ip BCS-Richardson model was independently constructed in [30] by other method.
The aim of the present article is threefold. The first aim is connected with the recent interest in the generalized Gaudin models (see e.g. [36][37][38], [39,40], [47], [51]). In the connection with this interest it becomes necessary to make more order in the subject, in particular to give a classification of the most physically important class of the non-skew-symmetric classical r-matrices, namely, diagonal in the standard basis classical r-matrices having the following form: r(u, v) = 3 α=1 r α (u, v)X α ⊗ X α . In this paper we show that there are only four classes of non-equivalent non-degenerated non-skew-symmetric so(3) ⊗ so(3)-valued classical r-matrices diagonal in the standard basis. They are: 3. Non-skew-symmetric elliptic r-matrix [18], 4. Skew-symmetric elliptic r-matrix of Sklyanin [3].
Observe that from this result follows that there are only two classes of the Cartan-invariant r-matrices, namely the "shifted" rational and "shifted" trigonometric r-matrices. In the present paper we also classify the shift elements having physical interpretations of the integrable magnetic fields in the corresponding generalized Gaudin models which are important in the BCS-Richardson theory.
The second aim of the present article is to obtain some new interesting generalizations of the BCS-Richardson model. By the virtue of the above classification there are only four possibilities for the BCS-Richardson model associated with the "diagonal" r-matrices. The first two of them -shifted rational and shifted trigonometric cases have been in details considered in our previous paper [22]. The fourth -skew-symmetric elliptic case is not useful for the BCS-type modelsthe corresponding r-matrix possesses no non-trivial shift elements and no non-trivial integrable magnetic fields. It is, hence, left to analyze the third -non-skew-symmetric elliptic case. The associated generalized Gaudin hamiltonians in the external magnetic field read as follows [18], [19], [21], [25], [26]: where Ŝ (l) α is α-th component of spin operator living in the site ν l , j α , α ∈ 1, 3 are the branching points of the elliptic curve and the representation of the spin algebra is arbitrary. The hamiltonians (1) mutually commute for all values of c α , j α and ν l (here ν l = ν k , ν l = −j α , k, l ∈ 1, N, α ∈ 1, 3): The corresponding "elliptic" BCS-Richardson hamiltonian is constructed from (1) by the fermionization of the spin operators Ŝ (l) α using the combination of the hamiltonians (1) in the limit j 3 → ∞. It has (up to the identity operator) the following simple form: where 1,l ≡ √ ν l + j 1 , 2,l ≡ √ ν l + j 2 , g = 1 c 3 and we have put for simplicity c 1 = c 2 = 0. Here c † i, , c j, , i, j ∈ 1, N, , ∈ +, − are standard fermion creation-annihilation operators.
The important problem is diagonalization of the generalized Gaudin and BCS hamiltonians, in particular, of the hamiltonians (1), (2). There are two standard methods used in order to find the spectrum of the Gaudin-type hamiltonians: algebraic Bethe ansatz and quantum separation of variables method [7]. The algebraic Bethe ansatz technique for the case of non-skew-symmetric classical r-matrices and Lie algebra so (3) was developed in our papers [20], [22]. Unfortunately its standard form works only for the Cartan-invariant classical r-matrices [20]. The standard form of the separation of variables method [7] works for a wider [29], but yet restricted class of sl(2) ⊗ sl(2) so(3) ⊗ so(3)-valued classical r-matrices that does not contain elliptic r-matrices. For the skew-symmetric elliptic r-matrix quantum separation of variables was performed in the paper [12]. The variable separation for the non-skew-symmetric elliptic r-matrix and the corresponding integrable models was performed in [27], [28] in the classical case. In the quantum case it is still missing.
A simple alternative approach to the diagonalization problem of the generalized Gaudin models has been proposed in recent papers [38], [39], [40]. Its idea goes back to the papers [43], [44], [45] and (in a certain sense) further to the paper [35]. The approach works in the simplest case when the values of all spin in the chain is equal to one-half and is based on the additional quadratic identities satisfied in this case by quantum Gaudin-type hamiltonians. So the third aim of the present article is to give the r-matrix interpretation and generalization of the above method. We show that this method works when the corresponding r-matrix satisfies the following conditions: where , 3 We show, that among the "diagonal" r-matrices the conditions (3) are satisfied by skewsymmetric rational r-matrix, two equivalent non-skew-symmetric trigonometric r-matrices and non-skew-symmetric elliptic r-matrix. We call the corresponding r-matrices non-skewsymmetric triad. It seems to be more important for the applications than the standard skewsymmetric triad of XXX, XXZ and XYZ Gaudin models. We show also that the method of "quadratic identities" can be prolonged onto the Cartan-invariant r-matrices not satisfying the conditions (3) if one includes into these identities also the number of particle operator M 3 .
We apply the results on the quadratic identities in order to find the spectra of the hamiltonian (2). The eigenvalues of the hamiltonian (2) are given by the formula: where h l are (up to the constants) eigenvalues of the "re-scaled" lim 3Ĥ l elliptic Gaudintype hamiltonians Ĥ l given by (1) in the representation of the spin algebra so(3) ⊕N with the highest wight λ = ( 1 2 , 1 2 , ...., 1 2 ). The eigenvalues h l satisfy the following system of quadratic equations: replacing the (yet unknown in this case) set of Bethe equations. Let us also remark that the discussed quadratic identities may be used in the proof of the completeness of the set of Bethe vectors (if known) (see [46] for the case of standard rational r-matrix).
The structure of this paper is the following: in the second section we remind the basic facts about the classical r-matrices with spectral parameters and perform the classification of the "diagonal" classical r-matrices and their shift elements. In the third section we review the generalized Gaudin models with and without external magnetic field and describe in details the generalized Gaudin models connected with the diagonal r-matrices. In the section four we obtain the quadratic identities on the generalized Gaudin hamiltonians and derive the conditions (3). In the fifth section we apply the obtained results to the BCS-Richardson models. At last in the sixth section we conclude and discuss the ongoing problems.

Generalities
Let g = so(3) be the Lie algebra of the three dimensional rotation group over C. Let X α , α ∈ 1, 3 be a standard basis in so(3) sl (2) with the commutation relations: We will use the following definition [8], [9], [10]: is called a classical r-matrix if it satisfies the generalized classical Yang-Baxter equation: where etc. and r αβ (u, v) are matrix elements of the r-matrix r(u, v).

Remark 2.
Note that integrable systems associated with equivalent r-matrices are equivalent.

Definition 2.
The r-matrix r(u, v) is called non-degenerated in an open region U if there is no such a subalgebra g ⊂ g and subspace g ⊂ g that r(u, v) ∈ g ⊗ g or r(u, v) ∈ g ⊗ g for all u, v ∈ U .
We will also use the following definition [26]: ∈ U is called a "special point" of the classical r-matrix if one of the following two conditions is satisfied: where g 0 is a subalgebra and g 0 is a subspace of g not coinciding with g itself.
Hereafter we will additionally assume that the function r(u 1 , u 2 ) is meromorphic in an open region U × U ⊂ C 2 and possesses the following decomposition: where r 0 (u 1 , u 2 ) is a holomorphic function with the values in so(3) ⊗ so (3), is a tensor Casimir:

Remark 3.
Observe that the classical r-matrices satisfying (8) are non-degenerate in the sense of the Definition 2. Moreover, for skew-symmetric r-matrices the regularity property (8) in an open region U is equivalent to the non-degeneracy in this region in the sense of the Definition 2 [5].
is called a symmetry algebra of the r-matrix if each X ∈ h satisfies the following equation: Definition 5. The r-matrix r(u, v) is called to be h-invariant if its symmetry algebra coincides with h. In particular, if h is a Cartan subalgebra, then the r-matrix is called to be Cartan-invariant.
We will need one more definition: c α (u)X α is called a generalized shift element (or simply "a shift element") if it satisfies the following equation: Remark 4. The non-trivial solutions of the equation (10) not always exist. Their existence is connected with the existence of the special points of the classical r-matrices r(u, v) [25], [26]. For skew-symmetric r-matrices any element of its symmetry algebra (if it exists) is a constant shift element. The shift elements (if exist) constitute a linear space.
Finally, the last definition we need in this section is the following one: (3)-valued r-matrix is so-called to be "diagonal" in the natural basis if it has the following form: Remark 5. Observe that standard skew-symmetric elliptic, trigonometric and rational r-matrices are diagonal in the natural basis of so(3).

Classification of the diagonal r-matrices
Without loss of generality we will hereafter assume that Cartan subalgebra of h ⊂ so(3) is spanned over the element X 3 . The diagonal r-matrices are then divided into two subclasses: Cartan-invariant "partially anisotropic" r-matrices and non-Cartan-invariant "completely anisotropic" r-matrices.
The following Theorem holds true: Up to the equivalence Cartan-invariant r-matrices possessing the decomposition (8) are divided into two subclasses: a) "Shifted rational" r-matrices: (ii) Up to equivalence there are only two Cartan non-invariant r-matrices with decomposition (8): a) Non-skew-symmetric elliptic r-matrix of the following form [18], [19], [27]: where u 2 α = u + j α , v 2 α = v + j α , α ∈ 1, 3 and parameters j α are arbitrary. b) Skew-symmetric elliptic r-matrix of the following form [3]: Idea of the Proof 1 . The proof of the item (i) of the Theorem repeats in its major part the lines of the proof of a Theorem 4.1 from [23], where a bit more general gl(2) case was considered. The idea of the proof is that the generalized Yang-Baxter equation in the component form is reduced to two non-equivalent quadratic equations on the functions r 3 (u, v) and r 1 (u, v) = r 2 (u, v), which, after the usage of the regularity condition (8) are reduced to one functional equation [23] for the details). It is then shown that the corresponding functional equation has two types of non-equivalent solutions, leading (after the equivalence transformation) to the r-matrices (12) and (13).
The statement of the item (ii) of the Theorem follows already from the results of [6] stating that there are only two completely anisotropic integrable models on so(3) * ⊕ so(3) * with the diagonal in the natural basis quadratic hamiltonians. They are Shottky-Frahm and Steklov models that correspond to the cases b) and a) of the item (ii).

Remark 6.
We call the r-matrix (14) to be elliptic because it possesses parametrization in terms of the elliptic Jacobi functions and j α are interpreted as the branching points of the elliptic curve: Nevertheless, the irrational parametrization used in this article is more convenient for the usage.
The above Theorem has the following important Corollary: Corollary 2.1. Up to the equivalence there exist only one so(3)-invariant classical r-matrix possessing the decomposition (8), namely the standard rational r-matrix, i.e. the r-matrix (12) with c 0 = 0.
From this Corollary follows, in particular, that the r-matrix implicitly present in the paper [34] is equivalent to the rational one and the corresponding Gaudin-type model is also equivalent to the standard rational Gaudin model.

Remark 7.
Note, that in the skew-symmetric case all non-degenerated classical r-matrices are classified [4]. They include three above diagonal r-matrices, as well as two non-diagonal rational and one non-diagonal trigonometric r-matrix [4].

Remark 8.
Observe that the families of the r-matrices (12), (13) have a skew-symmetric point c 0 = 0. As it is easy to show, after re-parametrization they coincide with the rational and trigonometric degeneration of skew-symmetric elliptic r-matrix of Sklyanin. Together with the elliptic r-matrix of Sklyanin they produce a famous triad of XXX, XXZ and XYZ Gaudin models.
Like the elliptic r-matrix of Sklyanin (15) the r-matrix (14) defines a triad of non-skewsymmetric classical r-matrices: the r-matrix (14) itself, its trigonometric and rational degenerations.
The following Proposition holds true: Proof. The Proposition is proven by the direct calculation. The item (i) is proven by the multiplication of the r-matrix (14) by √ v + j and by the change of the spectral parameters: To prove the item (ii) we at first observe that by the shift of the spectral parameter u → u − j , v → v − j we can make j 1 and j 2 = j 1 to be zero. Then multiplying the r-matrix by √ v + j 3 j 3 and changing the spectral parameters: Proposition is proven.

Remark 9.
Observe, that the rational r-matrix (12) with c 0 = 0 can be more simply obtained from the r-matrix (14) dividing it by √ j 1 j 2 j 3 and taking the limit j 1 , j 2 , j 3 → 0. The r-matrix equivalent to the shifted trigonometric r-matrix (12) with c 0 = ± 1 2 is obtained from the r-matrix (14) dividing it by √ j 3 , taking the limit j 3 → 0 and putting j 1 = j 2 = j . We will use this connection while considering the corresponding BCS hamiltonians.

Shift elements for the diagonal r-matrices
Let us now explicitly describe all possible shift elements for the diagonal classical r-matrices classified in the Theorem 2.1. They will be necessary for the construction of "integrable magnetic fields" in the corresponding Gaudin-type models.
Sketch of the Proof. The proof of the Proposition follows from the results of [25], [26]. In more details, the corresponding spaces of shift elements are identified with the dual spaces of the quotient of r-matrix algebras so (3) ν 0 where ν 0 is a special point of the classical r-matrix. The shifted rational r-matrix has only one special point u = ∞. Depending on the value of c 0 , the corresponding dual of the quotient space coincide either with the Cartan subalgebra or with so(3) itself. The shifted trigonometric r-matrix has two special points u = 0, ∞. Depending on the value of c 0 , the corresponding dual of the quotient space coincide either with the Cartan subalgebra or contains it, is three dimensional and its elements are described by the formula (17) [25]. Non-skew-symmetric elliptic r-matrix (14) has three special points u = −j α , α ∈ 1, 3 and each of the corresponding quotient spaces is one-dimensional [25]. In the result the corresponding generic shift element has the form (18). At last, elliptic r-matrix of Sklyanin has no special points, and, hence, no shift elements.

Remark 10.
Observe, that only the r-matrices of the described above non-skew-symmetric triad possess three-dimensional spaces of shift elements. They are special in this, and in some other senses too. In particular, they satisfy the additional equalities (3), which are useful in the diagonalization of the corresponding Gaudin-type hamiltonians.

Generalized Gaudin hamiltonians
Let Ŝ (k) a , α ∈ 1, 3, k ∈ 1, N be linear operators acting in some linear space that constitute a Lie algebra isomorphic to so(3) ⊕N with the commutation relations: The operator Ŝ (k) α is an α-component of the spin operator living in the k-site of the spin chain. Let ν k , ν k = ν l , k, l = 1, ..., N be some fixed points on the complex plane belonging to the open region U such that in U × U the r-matrix r(u, v) possesses the decomposition (8). Using this data, arbitrary non-skew-symmetric r-matrix satisfying (6) and possessing the decomposition (8) we have introduced in the papers [18], [19] the generalized Gaudin hamiltonians 2 : In the partial case of skew-symmetric r-matrices the second summand in the rhs of (20) disappears and hamiltonian (20) turns into the ordinary Gaudin hamiltonian. Like the ordinary Gaudin hamiltonians, the hamiltonians (20) mutually commute in arbitrary representation of spin algebra so(3) ⊕N (see [19] for the proof), defining in such a way quantum a integrable system.
The evident difference of the hamiltonians (20) from the standard Gaudin hamiltonians of [1] is in the presence in them of the additional second "self-interaction" term. It is interesting to remark, that there is an example when this second term disappears: Example 1. Let us consider the simplest possible case when the representation of the spin algebra so(3) ⊕N is given by the Pauli matrices (all representation with the spin equal to one-half): In this case due to the existence of the additional associative relations among the Pauli matrices: Nevertheless, the hamiltonians (21) do not coincide with the standard Gaudin hamiltonians: they have in general non-skew-symmetric interaction strengths r αβ (ν k , ν l ).
Remark 11. Some sporadic particular examples of the hamiltonians (20) that generalize standard Gaudin hamiltonians have appeared in the literature also previously (mainly) in the context of reflection equation algebras (see [31][32][33]). Nevertheless, until our papers [18], [19] their connection with the solutions of the equation (6) was neither observed nor understood.

Generalized Gaudin hamiltonians in a magnetic field
For the physical applications it is usually necessary to modify the hamiltonians (20) adding to it a linear term. How to do this in general without spoiling the integrability was explained in our paper [21]. It was shown that if the algebra-valued function c(u) =  (22) mutually commute in the arbitrary representation of the Lie algebra so(3) ⊕N . The hamiltonians (22) are the generalized Gaudin hamiltonians in the external magnetic field, where the role of external (non-homogeneous in general) magnetic field is played by generalized shift element c(u). Example 2. Let us consider the simplest possible case when the representation of the Lie algebra so(3) ⊕N is given by the Pauli matrices (all representation with the spin 1 2 ): S (l) α = σ (l) α , α ∈ 1, 3, l ∈ 1, N. In this case the hamiltonians (22) acquire (up to the identity operator) the Gaudin-like form: r αα 0 (ν l , ν l )I d. Their spectra also differ by this term.
In the case of diagonal r-matrices generalized Gaudin hamiltonians in magnetic field are simpler: This partial case of the hamiltonians (22) was considered in the papers [39], [40]. The coincidence of the hamiltonian (25) with that of [39], [40] is achieved by a substitution: r α (ν k , ν l ) = α kl , c α (ν l ) = B α l .

"Diagonal" Gaudin-type hamiltonians: case by case study
Having at hand the classification of the diagonal so(3) ⊗ so(3) valued r-matrices we can (up to the equivalence) explicitly describe all the corresponding Gaudin-type hamiltonians with and without external magnetic field, i.e. specify the hamiltonians (24) for all the above four classes of r-matrices.

Shifted rational case
The Gaudin Hamiltonians in magnetic field corresponding to the r-matrix (13) have the form: Due to their U(1) symmetry they commute with the "global spin operator" M 3 : With the help of this operator one can write the Hamiltonian (26) as follows: 3 , l ∈ 1, N. Here Ĥ rat l is a standard rational Gaudin hamiltonian without external magnetic field corresponding to the choice c 0 = 0, c 3 = 0 in the hamiltonian (26). By the other words, the appearance of c 0 in the r-matrix leads to the appearance of the "dynamical" external magnetic field in the corresponding rational Gaudin hamiltonian. The operator M 3 is connected with Ĥ l as follows: In the case c 0 = 0 we obtain standard rational Gaudin hamiltonians in the constant magnetic field:Ĥ It is Cartan-invariant with respect to the chosen Cartan subalgebra only in the case c 1 = c 2 = 0.

Shifted trigonometric case
The Gaudin Hamiltonians in magnetic field corresponding to the r-matrix (13) in the case c 0 = ± 1 2 have the form: Due to their U(1) symmetry they commute with the "global spin operator" (27). Similar to the rational case one can write the Hamiltonian (29) as follows: Here Ĥ trig l is a standard "skew-symmetric" trigonometric Gaudin hamiltonian without external magnetic field corresponding to c 0 = 0, c 3 = 0 in the formula (29) and M 3 is a global spin operator. By the virtue of (30) the relation (28) holds true also in this case.
In the case c 0 = − 1 2 the Gaudin-type hamiltonians in external magnetic field have the form: Such the Gaudin-type model with non-diagonal magnetic field has been studied in the papers [36], [38], [42]. For the case of the diagonal magnetic field they have been studied previously in [24].

Remark 13.
One can connect arbitrary hamiltonian of the "shifted trigonometric series" also with the hamiltonian (29). In general hamiltonians Ĥ l of the "shifted series" differ from each other by the "dynamical magnetic field term" proportional to M 3Ŝ (l) 3 . By other words the models associated with the r-matrices that are connected with each other by the tensorial shift c 0 X 3 ⊗ X 3 are very close, but non-equivalent. They have different (although similar) algebras of integrals, different (although similar) equations for the spectra, may have different spaces of shift elements (admissible magnetic field) etc.

Non-skew-symmetric elliptic case
Let us consider the most interesting non-skew-symmetric elliptic case. The direct calculation gives:Ĥ In the case when Ŝ (l)

Skew-symmetric elliptic case
Let us finally briefly consider the case of the elliptic skew-symmetric r-matrix. The standard elliptic "XYZ" Gaudin hamiltonians correspond to it: Unfortunately elliptic "XYZ" Gaudin hamiltonians (34) will be of no use for us because they do not have integrable magnetic fields and do not have several other good properties, in particular, they do not satisfy the additional quadratic identities that will be considered in the next section.

Additional identities among the Gaudin-type integrals
In this section we will describe additional identities existing for some of the Gaudin-type hamiltonians (20) when the values of all spin in the chain are equal to one-half. They are used in order to find the quadratic identities for the eigenvalues of the Gaudin-type hamiltonians (20). These equations replace the famous Bethe equations which are known not for all classical r-matrices.

The case of the generalized Gaudin models
The following Proposition holds true: Proposition 4.1. Let the classical r-matrix r(u, v) satisfies the following two conditions: where , 3 is a scalar product in the third tensor component and the bilinear operation [[ , ]] denotes simultaneous commutator in the both components of the tensor product, e.g.

The case of the generalized Gaudin models in magnetic field
Let us now consider the generalized Gaudin models in the external magnetic field and find the analog of the equation (36).
The following Proposition holds true: Let the r-matrix satisfies the conditions (35) and its shift element c(u) satisfy the following additional condition: where , 2 is a scalar product in the second multiplier of the tensor product. Then the hamiltonians (23) satisfy the following quadratic identities: (r αβ (ν k , ν l )) 2 I d, l ∈ 1, N.

Remark 15.
Observe the covariance of the conditions (37) with respect to the three equivalence transformations of the r-matrix.

The case of the diagonal r-matrices
Let us consider the most interesting case of the diagonal classical r-matrices and find out which of them satisfy the conditions (35), so that the equation (36) be applicable to the corresponding Gaudin-type hamiltonians.
The following Theorem holds true: (3) valued r-matrices satisfying the condition (35) are the anisotropic non-skew-symmetric r-matrix (14) and the r-matrices that are equivalent to its trigonometric and rational degenerations.
Proof. In order to prove the theorem it is necessary to re-write the conditions (36) for the case of the diagonal r-matrices. They are written in the component form as follows: We will prove the Theorem using our classification of the r-matrices and case by case calculation. Let us at first consider the r-matrix (14). The direct calculation show that i.e.
The direct check shows that for the r-matrix (12) the conditions (39) are satisfied only for The direct check shows also that conditions (39) are satisfied for the r-matrix (13) only for c 0 = ± 1 2 and in this last cases we have: Finally, using the properties of the elliptic Jacobi functions one shows that skew-symmetric elliptic r-matrix of Sklyanin satisfies neither the condition (39a) nor the condition (39b).
This proves the Theorem.
Let us now check whether the r-matrices satisfying the condition (35) satisfy the condition (37).
The following Proposition holds true: The non-skew-symmetric elliptic r-matrix (14) and the r-matrices that are equivalent to its trigonometric and rational degenerations satisfy the conditions (37).
Proof. The Proof is achieved by the direct calculation. Indeed, in the case of the diagonal r-matrices one has that the condition (37) reads as follows: Substituting in (42) the components of the r-matrix (14)  In the same direct way it is checked that the shifted rational and shifted trigonometric r-matrices in the points c 0 = 0 and c 0 = ± 1 2 correspondingly satisfy the conditions (42). Proposition is proven. Now let us present the equations (38) for the above "non-skew-symmetric triad" case by case. 1) For the case of the r-matrix (12) with c 0 = 0 the equations (38) reads as follows: 2) For the case the r-matrix (13) with c 0 = − 1 2 the equations (38) reads as follows: 3) For the case of the r-matrix (14) the equations (38) reads as follows: where

Additional identities for other models with diagonal r-matrices
Shifted rational r-matrices with c 0 = 0 and shifted trigonometric r-matrices with c 0 = ± 1 2 (including standard trigonometric r-matrix with c 0 = 0) do not satisfy the conditions (39). Nevertheless, for the corresponding Gaudin hamiltonians in the representations with all spins equal to one-half, it is sill possible to write the additional identities. Indeed, as it follows from the results of the previous section the generalized Gaudin hamiltonians Ĥ l whose r-matrices differ by the term proportional to X 3 ⊗ X 3 differ from each other by the "dynamical magnetic field term" proportional to M 3Ŝ (l) 3 . Hence, taking as a starting point the generalized Gaudin hamiltonians, corresponding to the rational r-matrix with c 0 = 0 and trigonometric r-matrix with c 0 = − 1 2 and external magnetic field with c 1 = c 2 = 0 we obtain from (43) and (44) correspondingly, the following equations: 1) for the case of the r-matrix (12): 2) For the case the r-matrix (13): where we have taken into account that M 3 is an integral of motion, commuting with Ĥ l independently of the values of c 0 and commuting with all Ŝ (l) 3 . It behaves as a scalar with respect to them and does not influence the proof of the Proposition 4.2, leading only to the appropriate shift of the constant c 3 in the equations (43) and (44).

Remark 16.
Observe that using the relation (28) one can exclude the integral M 3 from the equations (46) and (47)

Fermionization
Having obtained quantum integrable spin systems it is possible to derive, using them, integrable fermionic systems. For this purpose it is necessary to consider realization of the corresponding spin operators in terms of fermionic creation-annihilation operators. We will consider here only simplest "fermionization" of the Lie algebra so(3) ⊕N corresponding to its representation with all the spins in the chain equal to one half.

Elliptic case
The hamiltonian To obtain the most physical BCS-Richardson's hamiltonian it is necessary to chose the "right" combinations of the generalized Gaudin hamiltonians in the external magnetic field. We will consider the following combination of the hamiltonians (33): After the simple transformation we obtain: Now for the "elliptic BCS hamiltonian" we will take the following singular limit of Ĥ : where we have also re-scaled the "shift constant" c 3 → j 3 c 3 , c 2 → √ j 3 c 2 , c 1 → √ j 3 c 1 . In the result we obtain the following hamiltonian: The terms proportional to c 1 , c 2 in it are interpreted, similar to [36], [38] as interaction with the environment.
Using the relation between so(3) and sl(2) basis we obtain the following form of (53): (c 1 ν l + j 2 + ic 2 ν l + j 1 )σ Now, applying the fermionization formulae and putting 1,l ≡ ν l + j 1 , 2,l ≡ ν l + j 2 we obtain (up to the identity operator) the following integrable fermion hamiltonian: We call it "open" elliptic BCS-type hamiltonian. In the case c 1 = c 2 = 0 it acquires "closed" form:Ĥ The spectra The eigenvalues of the hamiltonian (53), by its definition, are given by the formula: where h l are the eigenvalues of the "re-scaled" elliptic Gaudin-type hamiltonians lim that, by the virtue of the results of the previous section satisfy the following system of quadratic equations: where we have taken in the equations (45) the limit j 3 → ∞, previously dividing the equation (45) by j 3 and re-scaling appropriately c α , α ∈ 1, 3 and Ĥ l .

Remark 17.
The equations (57) replace the Bethe equations unknown in this case.

Trigonometric degeneration
The hamiltonian In the trigonometric degeneration j 1 = j 2 = j the elliptic hamiltonian (54) has the form After the fermionization and substitution l = √ ν l + j it acquires "open" p + ip BCS form: To obtain standard "closed" p + ip BCS hamiltonian one has simply to put in (59) c 1 = 0, c 2 = 0: The spectra The eigenvalues of the hamiltonian (58) are given by the formula: where the eigenvalues h l of the "trigonometric" Gaudin-type hamiltonians Ĥ l satisfy the following system of quadratic equations: These equations coincide with the equations (44) after the change of variables ν k → ν −2 k − j .

Rational degeneration
The hamiltonian In order to obtain the rational BCS-Richardson hamiltonians one should consider another combination of the elliptic Gaudin-type hamiltonians and take the following limit: where we have also to re-scale the constants c α : c α → √ j α √ j 1 j 2 j 3 c α . In the result we obtain the following combination of the rational Gaudin hamiltonians: 3 α=1 ν l c α σ (l) α , l ∈ 1, N.
It can be put to the standard Richardson-type form only in the case c 1 = 0, c 2 = 0 when M 3 is an integral of motion. Putting c 1 = 0, c 2 = 0, extracting from Ĥ the integral M 2 3 and passing to sl(2) basis we obtain (up to the identity) the s-type BCS-Richardson hamiltonian written in the spin form: where l ≡ ν l . The integral M 3 is connected with Ĥ l as follows: The spectra The eigenvalues of the hamiltonian (64), by its definition, are given by the formula: where h l are the eigenvalues of the "re-scaled" elliptic Gaudin-type hamiltonians 3Ĥ l coinciding with the rational Gaudin hamiltonians in the constant magnetic field and satisfying the following system of quadratic equations:

Conclusion and discussion
In this article we have considered generalized Gaudin spin chains with [21] and without [18], [19] external magnetic field in the simplest case of g = so (3) and presented the classification theorem for the diagonal in the natural basis r-matrices and their "shift elements" playing the role of external magnetic fields. We have applied these results to the description, classification and solution of BCS-Richardson-type models. In particular, we have introduced "elliptic" BCS-Richardson hamiltonian corresponding to non-skew-symmetric elliptic r-matrix.
We have developed and generalized the technique of the spectrum calculation based on the additional quadratic identities among the Gaudin hamiltonians proposed (in the partial cases) in [39], [40] finding the general conditions on the corresponding r-matrix to be hold true for their existence. We have shown that for all the Gaudin models based on the diagonal r-matrices, except for the XYZ Gaudin model based on the elliptic r-matrix of Sklyanin, this technique is applicable. Using this technique we have calculated, in particular, the spectrum of our elliptic BCS-Richardson hamiltonian in terms of solutions of the set of quadratic equations for the spectrum of the corresponding Gaudin-type hamiltonians. It will be interesting to find out for what non-diagonal non-skew-symmetric r-matrices the discovered conditions are satisfied. For this purpose it will be necessary to classify such the r-matrices. This problem is open.
Another interesting ongoing problem is to generalize the technique of spectrum calculation based on the additional quadratic identities among the Gaudin hamiltonians onto the case of integrable spin-boson models of the Jaynes-Cummings-Dicke type. As we have shown in our papers [48], [49], [50] such the models exist for a wide class of the non-skew-symmetric classical r-matrices, in particular for all diagonal r-matrices (except for the elliptic r-matrix of Sklyanin). Observe that in the case of the standard rational r-matrix the existence of such the quadratic identities follows from the results of [35]. Their existence for the Jaynes-Cummings-Dicke type models associated with the skew-symmetric trigonometric r-matrices was shown in [41]. In order to obtain the quadratic identities among the Jaynes-Cummings-Dicke type hamiltonians based on the general r-matrix it will be necessary to impose yet another additional constraint on this r-matrix. We conjecture that such the additional quadratic identities exist for the case of the models based on the shifted rational and shifted trigonometric r-matrices but do not exist for the Jaynes-Cummings-Dicke type model based the non skew-symmetric elliptic r-matrix [49], [50], [28].