Non-skew-symmetric classical r -matrices and integrable “ p x + ip y ” proton–neutron BCS models

We study integrable cases of the pairing BCS Hamiltonians containing several types of fermions and possessing non-uniform coupling constants. We prove that there exist three classes of such the integrable models associated with “ Z 2 -graded” non-skew-symmetric classical r -matrices with spectral parameters and Lie algebras gl ( 2 m) , sp ( 2 m) and so ( 2 m) , respectively. The proposed models are higher rank generalizations of the so-called “ p x + ip y ” one-type fermion ( m = 1) BCS model. In the partial case of two types of fermions ( m = 2) the obtained models may be interpreted as N = Z , “ p x + ip y ” proton–neutron integrable models. In particular, in the case of sp ( 4 ) we obtain the “ p x + ip y ”-analogue of the famous integrable proton–neutron model of Richardson. We ﬁnd the spectrum of the constructed Hamiltonians in terms of solutions of Bethe-type equations.


Introduction
Proton-neutron pairing Hamiltonian plays an important role in the description of nuclear structure for proton-rich nuclei with N Z. In contrast to the proton-proton and neutronneutron pairing, the proton-neutron pairing may exist in two different forms, namely in the forms of the isoscalar (T = 0) and isovector (T = 1) pairing. A pairing formalism, which in-cludes T = 0 and T = 1 p-n correlations was proposed in [1]. Such the pairing models were studied in [2,3] using the mean-field approximation.
In our previous paper [5] we have constructed three types of integrable fermion Hamiltonians containing isoscalar (T = 0) and isovector (T = 1) pairing and possessing uniform coupling constants. The integrability of the models has permitted us to solve them exactly using the algebraic Bethe ansatz.
In the present paper we investigate integrable cases of more general pairing Hamiltonians containing isoscalar (T = 0) and isovector (T = 1) pairing and possessing non-uniform couplings. The Hamiltonian of the model under the consideration has the following form: where c † l,p,σ , c l,p,σ and c † l,n,σ , c l,n,σ are fermion creation-annihilation operators corresponding to protons and neutrons, σ ∈ ± is a sign of the angular momentum projection, l,p , l,n are the proton and neutron single particle energies, N is the number of protons and neutrons (their numbers are taken to be equal) and G T =0 pn,kl , G T =1 pn,kl are the interaction strengths of proton-neutron isoscalar and isovector pairings, G T =1 pp,kl and G T =1 nn,kl are interaction strengths of proton-proton and neutron-neutron isovector pairings.
Like in the case of uniform coupling constants [5], it turned out to be methodologically more convenient to consider a more general Hamiltonian, namely, the Hamiltonian containing instead of two types of fermions (i.e. instead of protons and neutrons only) an arbitrary number m of fermion types. In more details in the present paper we study the following Hamiltonian: the Hamiltonian (2) admits, similar to the case of the uniform coupling constants [5], three integrable cases. The form of the non-uniform coupling constants G T =1 ij,kl , G T =0 ij,kl in these cases is the following: l . These three cases are connected with the Lie algebras gl(2m), sp(2m) and so(2m), respectively. To the best of our knowledge for m > 1 all these integrable cases are new.
Let us briefly describe our method of obtaining of the above integrable cases of the Hamiltonian (2).
The usual method of obtaining of integrable fermion pairing models is based on applications of Gaudin spin models [7] in an external magnetic field. The idea to associate a "one type of fermions" BCS model with usual Gaudin spin models based on a standard skew-symmetric rational r-matrix and the Lie algebra sl (2) is attributed to the paper [8]. The idea to associate a "many type of fermions" BCS model with rational Gaudin spin models based on a skew-symmetric rational r-matrix and higher rank Lie algebras ascends to the paper [11]. It was used in many papers (see [9,10,12,13,5]).
Unfortunately, all the integrable fermion models that may be obtained using standard classical rational r-matrices possess uniform coupling constants. The integrable fermion models associated with the skew-symmetric trigonometric r-matrices are characterized by the Hamiltonians containing not only pairing interaction terms, but also other terms having no evident physical interpretation (see e.g. [14]). At last, there is no integrable Gaudin spin models in an external magnetic field associated with a skew-symmetric elliptic r-matrix and no corresponding integrable fermion model. Moreover there exists no other classical skew-symmetric r-matrices except for the rational, trigonometric and elliptic ones due to the Belavin-Drinfeld theorem [4].
Hence, in order to generalize the above approach and obtain more general integrable BCStype models, in particular, the models containing non-uniform coupling constants, it is necessary to generalize the very integrable Gaudin spin models in an external magnetic field. It was done in a series of our previous papers [16][17][18][19]23]. The main our idea is a possibility to associate quantum integrable Gaudin-type spin chains not only with skew-symmetric but also with non-skew-symmetric classical r-matrices. The non-skew-symmetric classical r-matrices satisfy a "generalized" classical Yang-Baxter equation instead of the usual one. This permits us to overcome the "bottle-neck" connected with Belavin-Drinfeld Theorem [4] and to obtain new integrable Gaudin-type spin models.
In our papers [20][21][22] we have applied the previously discovered generalized Gaudin spin models based on non-skew-symmetric r-matrices and Lie algebras sl (2) (and gl (2)) and have constructed integrable "one type of fermions" BCS models. We have shown in [21,22] that among the proposed models there is one especially interesting model with the Hamiltonian possessing only kinetic and pairing interaction terms with non-uniform coupling constants. This model was soon re-discovered by other method in [25] (see also [26]) and called p x + ip y integrable model.
In the present paper we generalize our idea of associating of integrable fermion models with non-skew-symmetric classical r-matrices onto the cases of higher rank Lie algebras g. For this purpose we describe the generalized Gaudin spin models in an external magnetic field [17][18][19] and study their symmetries. The symmetries of these models are very important -not all possible generalized Gaudin spin models in an external magnetic field are suitable for a construction of the generalized BCS models, but only those that are invariant with respect to the Cartan subalgebra h ⊂ g. This guaranties the existence of r = dim h particle number operators commuting with a Hamiltonian and permits one to interpret these models as possessing r different types of fermions.
We show that the "geometric" symmetry algebra of the generalized Gaudin spin models in an external magnetic field coincides with the intersection of the geometric symmetry algebra of the corresponding r-matrix and geometric symmetry algebra of "an external magnetic field". From this it follows that for a construction of BCS-type models only h-invariant classical r-matrices and h-valued external magnetic fields are useful. This substantially restricts the admissible class of the r-matrices. Nevertheless it is still large enough (see e.g. [16,17]).
In the present paper we concentrate on one class of such the r-matrices, namely, those leading to many type of fermions "p x + ip y " integrable models. We call such the r-matrices Z 2 -graded because they are connected with Z 2 -grading of the corresponding Lie algebra g, i.e. with the decompositions g = g0 + g1, where g0 is a closed subalgebra and g1 is a g0-module. For the construction of the BCS model one should require, moreover, that g0 = g K 0 and g1 = g K 1 + g K −1 , where K is an element of the Cartan subalgebra, g K 0 is its reductive centralizer in g and the subspaces g K ±1 are abelian subalgebras. It turned out that the integrable magnetic field for the corresponding r-matrices is proportional to the same element K of the Cartan subalgebra 1 which is used in the definition of the r-matrix itself. In this case the symmetry algebra both of the Z 2graded r-matrix and element K coincides with the subalgebra g K 0 . Having a fixed non-skew-symmetric classical r-matrix and "integrable magnetic field" one completely determines (at the Lie-algebraic level) the corresponding generalized Gaudin system and its integrals. We show that among the integrals of the model there is the following remarkably simple Hamiltonian: where is a system of roots of g and K is the system of roots of g K 0 , K = rank g i=1 k i H i ∈ h is the above element of the Cartan subalgebra, δ K is a sum of roots over the set ( / K ) + , and the operatorsŜ (l) i ,Ŝ (l) ±α , i ∈ 1, rank g, α ∈ represent the Cartan-Weyl basis of the lth copy of g in a representation of the direct sum Lie algebra g ⊕N and l are connected with the poles of the Lax operator of the generalized Gaudin system. The Hamiltonian (3) is our generalized g-valued p x + ip y BCS Hamiltonian written in the spin form. Due to the said above the symmetry algebra of the Hamiltonian (3) coincides with the reductive subalgebra g K 0 of the "global" spin operators. Such the symmetry algebra always contains the global Cartan subalgebra.
In order to obtain from the integrable spin Hamiltonian BCS-type fermion Hamiltonian (3) it is necessary to chose the fermionization formulas for the Lie algebra g (more precisely, a direct sum Lie algebra g ⊕N ) that should automatically define the representations of each Lie algebra g in g ⊕N and its highest weight λ (l) , l ∈ 1, N. We will in a standard way consider the case λ (l) = λ, ∀l ∈ 1, N, i.e. use the same fermionization formulas for all spins in the chain. In order to obtain the integrable Hamiltonian (2) we restrict ourselves to the consideration of the Lie algebras gl(2m) and their subalgebras so(2m) and sp(2m) and for the representations with the following highest weights: λ = (λ 1 , . . . , λ 2m ) with λ 1 = · · · = λ m = 1, λ m+1 = · · · = λ 2m = 0 for gl(2m) and λ = (λ 1 , λ 2 , . . . , λ m ) with λ 1 = · · · = λ m = 1 for the sp(2m) and so(2m). Using the fermionization formulas existing for these representation [5] we finally obtain the integrable cases of the l respectively, where the interaction strength g is connected with the coefficients of the "external magnetic field" -the chosen diagonal matrix K. Let us mention that the constructed Hamiltonians possess the geometric symmetry algebra, namely the algebra gl(m) ⊕ gl(m) in the first case, the algebra gl(m) in the second and the third cases. 2 The symmetry algebra coincides with the global spin algebra g K 0 written in the given fermion representation.
We diagonalize the constructed Hamiltonians using the general answers for the spectrum of the Z 2 -graded Gaudin Hamiltonians in terms of solutions of the Bethe-type equations in any representation of the Lie algebra g ⊕N [15]. It turned out that for all the three considered cases the spectrum of the integrable Hamiltonian (2) is given by the same formula: are the mth set of solutions of the Bethe equations, that depend on the Lie algebra, element K and highest weights λ (l) . For the cases at hand the Bethe equations have a sufficiently simple form and we write them explicitly in each case. Due to the existence of the non-abelian symmetry algebra the constructed energy levels are degenerated: the multiplicity of the degeneration of the level is equal to the dimension of the representation of the symmetry algebra g K 0 and depends on the given eigen-vector.
In the end of the introduction let us explain the usage of the name "p x +ip y " for the introduced integrable fermionic models. Using the fact that the considered subalgebras g K ±1 are abelian one may re-scale the spin operators as follows: , possessing the so-called "p x + ip y wave symmetry" [25].
The structure of the present paper is the following: in the second section we remind the necessary facts about semisimple Lie algebras and their reductive subalgebras, in the third section we describe the needed fermionization formulas for the Lie algebras gl(n), sp(2m) and so(2m). In the fourth section we consider quantum systems governed by general non-skew-symmetric classical r-matrices, in particular, the generalized Gaudin systems in an external magnetic field, and find their symmetry algebras. In the fifth section we concentrate on Z 2 -graded r-matrices and the corresponding Gaudin Hamiltonians. In the six section we construct integrable many types of fermions "p x + ip y " BCS Hamiltonians (both in the spin and fermion form). Finally in the seventh section we find the spectrum of the constructed Hamiltonians in terms of solutions of Bethe-type equations.

Simple Lie algebras and their reductive subalgebras
The results of the present paper are obtained using Lie-algebraic technique. That is why in this section we will remind the necessary facts from the theory of simple Lie algebras and their subalgebras [27,28].

General case
Let g be a complex simple Lie algebra or reductive Lie algebra gl(n). Let X a , a = 1, dim g be a basis in g with the commutation relations: where C c ab are the structure constants of Lie algebra. Let ( , ) be a bilinear symmetric invariant form on g (Killing-Cartan form), g ab are the components of the nondegenerate invariant form: g ab = (X a , X b ). Hereafter all raising and lowering of vector and tensor indices will be made with the help of the metric tensor with the components g ab and its inverse with the components g ab .
In the present paper we will use mainly the root basis in g. In more details, let be the set of roots of g, + be the set of its positive roots, g α be the root space corresponding to the root α, X α be a basis element of this root space and H i be a basis vector of Cartan subalgebra h. The basis in g consists of the elements {X α , α ∈ ; H i , i ∈ 1, rank g}. The corresponding commutation relations have the following form: where The "orthonormality" relations read in this basis as follows: Remark 1. Note, that hereafter the Latin indices a, b, c, d will index the basic elements of g independent of a division of them on root spaces and Cartan elements, indices i, j , etc., will index a basis in Cartan subalgebra, Greek indices α, β, γ , δ will denote the roots belonging to . We will need also a description of the reductive subalgebras of the algebra g. Let K ⊂ be a closed, symmetric subset of the set of all roots. Then there exist an element K = rank g i=1 k i H i such that α(K) = 0 for α ∈ K . The reductive subalgebras of g correspond to the elements K and the subset of roots K . We will denote such the subalgebras by g K 0 . The basis of these subalgebras consists of the elements {X α , α ∈ K ; H i , i ∈ 1, rank g} (rank g = dim h).
Each reductive subalgebra g K 0 provides its own "triangular" decomposition of g: ⊂ g K ±1 and the subalgebras g K ±1 are nilpotent subalgebras generated by the following basic elements: {X α , α ∈ ( / K ) ± }.
The reductive subalgebra g K 0 consists of its semisimple part -the ideal [g K 0 , g K 0 ] and the center -the abelian subalgebra c(g K 0 ) ⊂ h. Let us also note that using the definition of g K 0 it is easy to show that the element K = rank g i=1 k i H i belongs to c(g K 0 ).

Case of classical matrix Lie algebras
Let us briefly illustrate the above theory by examples of classical matrix Lie algebras, their root systems and reductive subalgebras. We will consider the examples of the elements K such that the subalgebras g K ±1 in the corresponding triangular decomposition are abelian. Only such the cases will be used for the construction of BCS-type systems.

Case of g = gl(n)
The algebra gl(n) consists of all n by n matrices X. The natural basis of gl(n) consists of the elements X ij , i, j ∈ 1, n which are n × n matrices with the matrix elements (X ij ) ab = δ ia δ jb . The Cartan subalgebra coincides with the algebra of diagonal matrices, and H i ≡ X ii , i ∈ 1, n is its orthonormal with respect to the form (X, Y ) = tr XY basis. The set of all roots coincides with the linear forms α ij = α i − α j , where α i (H j ) = δ ij , and the corresponding properly normalized element of the root space g α ij is X α ij ≡ X ij (i = j ). An arbitrary element of the Cartan subalgebra has the form K = n i=1 k i X ii and α ij (K) = k i − k j . The set of simple roots {α i , i ∈ 1, n − 1} are the roots α ii+1 , i ∈ 1, n − 1.

Case of g = sp(2m)
This Lie algebra consists of the matrices X of a dimension n = 2m such that Xw = −wX, where w is a matrix of the symplectic form in the linear space of dimension 2m: w = 0 1 m −1 m 0 . The corresponding matrix X is written explicitly as follows: X = a 11 a 12 a 21 a 22 , where the submatrices a ij act in the space of the dimension m and satisfy the following conditions: The Cartan subalgebra coincides with the set of the diagonal matrices, namely: This basis in the Cartan subalgebra is orthonormal with respect to the form (X, Y ) = 1 2 tr XY . The set of positive roots coincides with the linear forms The corresponding normalized basic elements of the root spaces g ±α ± i,j , g ±2α i are given by the formulas: where X ij is a basic element of gl(2m) defined as in the previous example. The set of the simple . The set of roots K has the form: are abelian and have the form:

Case of g = so(2m)
This Lie algebra consists of the matrices X of dimension 2m such that Xs = −sX, where s is a matrix of the symmetric bilinear form in the linear space of dimension 2m: s = 0 1 m 1 m 0 . The corresponding matrix X is written explicitly as follows: X = a 11 a 12 a 21 a 22 , where the sub-matrices a ij act in the space of the dimension m and satisfy the following conditions: The Cartan subalgebra coincides with the set of the diagonal matrices, namely: This basis in the Cartan subalgebra is orthonormal with respect to the form (X, Y ) = 1 2 tr XY . The set of positive roots coincides with the linear forms The corresponding normalized basic elements of the root spaces g ±α ± i,j , g ±2α i are given by the formulas: where X ij is a basic element of gl(2m) defined as previously. The set of the simple roots . The set of roots K has the form: The nilpotent subalgebras so(2m) K ±1 are abelian and have the form: so(2m) The corresponding subset of roots are:

Fermionization
The important technical tool in the construction of integrable BCS models is the "fermionization" of integrable spin Hamiltonians based on Lie algebras. By the fermionization we will mean the realization of complex Lie algebras (or direct sums of their N -copies) via n (or nN ) fermion creation-annihilation operators (i.e. via complex Clifford algebra in the space of the dimension 2n (or 2nN )). The commutation relations of the fermions will be the standard anticommutation relations. In more details, let c k,i ,c † l,j , k, l ∈ 1, N, i, j ∈ 1, n be the fermionic creation-annihilation operators, i.e: We will use the standard representation of the fermion algebra with the vacuum vector |0 such that The representation space is spanned over the elements There are different fermionization formulas that depend on a chosen Lie algebra and a chosen class of its representations. In the next subsections we will construct the "fermionization" of certain classes of representations of classical matrix Lie algebras that we will use in the present paper.

Remark 2.
Observe that in the applications the index i will label different types of fermions and the index l will label fermions of the same type.

Case of g = gl(n)
In this subsection we describe a fermionic realization of the fundamental representations of g = gl(n).
The following proposition holds true [5]: N satisfy the anticommutation relations (10). Let us fix some integer p ∈ 1, n − 1. Then (i) the following operators: constitute a fermion realization of the direct sum Lie algebra gl(n) ⊕N , i.e.: sj .

Example 1.
In the simplest gl(2) ⊕N case the only meaningful choice is p = 1 and we obtain the following fermionization formulas [21]: After the restriction to the subalgebra sl(2) ⊕N and the substitutionŜ 3 we obtain the standard fermionization of the Lie algebra sl(2) used in the BCS models [8,14].

Case of g = sp(2m)
The Lie algebra sp(2m) is important for the applications in nuclear physics. In this subsection we describe its fermionization. For this purpose it is convenient to use the root basis of sp(2m) described in Section 2.2.2 rather than the matrix one.

Example 2.
In the simplest sp(2) ⊕N case the only non-trivial generators areŜ ±2α 1 which are written in terms of the fermionization formulas equivalent to those discussed in Example 1. This equivalence is explained by the well-known isomorphism: sp(2) sl (2).

Example 3.
In sp(4) ⊕N case, i.e. in when m = 2 we obtain the following fermionization formulas:Ŝ where l ∈ 1, N. This fermionization is equivalent to the fermionization constructed in [9] for g = so (5). It is explained by the isomorphism sp(4) so (5) existing in this case.

Case of g = so(2m)
The Lie algebra so(2m) may also be important for the applications. In this subsection we describe the fermionization of so(2m) based on its realization as a subalgebra in gl(2m). This fermionization will be similar to the one written for a symplectic Lie algebra.

Remark 3.
Observe that the operators (14) in this case, like in the case of g = gl(n), do not span the whole fermion Fock space (which is completely reducible so(2m) ⊕N -module) but only its irreducible submodule containing the vacuum vector.

General classical r-matrices and "shift elements"
In the construction of BCS models classical r-matrices are used. That is why we have to remind several main facts about them. We will need the following definition [29][30][31]: Definition 1. A function of two complex variables r(u 1 , u 2 ) with values in the tensor square of the algebra g is called a classical r-matrix if it satisfies the following "generalized Yang-Baxter equation": where (15) pass into the usual classical Yang-Baxter equation [4]: Hereafter we will not be interested in the global geometry associated with the solutions of (15) preferring to work in some local parametrization. In more details, we will assume that the parameters u and v are such that in some open region U ⊂ C 2 the r matrix r(u, v) possesses the decomposition: where Ω ∈ g ⊗ g is the tensor Casimir: Ω = dim g a,b=1 g ab X α ⊗ X b and r 0 (u, v) is a regular in U g ⊗ g-valued function, i.e. is decomposed into a Taylor power series in u and v.
Lat g 0 be a Lie subalgebra of g. We will use the following definition: if for all X ∈ g 0 the following identity holds true: Note, that on the level of Lie groups the property (18) means exactly the invariance: for g ∈ G 0 , where G 0 is a Lie group of g 0 .
We will need also the following definition [19]: a=1 c a (u)X a is called a "generalized shift element" if it solves the following equation: Remark 5. Observe that the generalized shift elements constitute a linear space, i.e. one can add them and multiply by complex numbers. Observe also that for the skew-symmetric g 0 -invariant r-matrices (and for some special non-skew-symmetric ones) any element X ∈ g 0 is a shift element. For general non-skew-symmetric r-matrices shift elements are not connected with the symmetry of the r-matrix.

Linear algebra of the Lax operators
Having general solution of Eq. (15) one can define the "Lie algebra of quantum Lax operators" by one convoluted "tensor" formula: Remark 6. The form of the Lax operator as a function of quantum dynamical variables depends on the physical model under the consideration. In this paper we will consider the most important form of such the Lax operators corresponding to generalized Gaudin systems in an external magnetic field.
The role of the generating function of the quantum integrals usually plays the function trL 2 (u). For the subsequent we will need to know its symmetry algebra, coinciding with the symmetry algebra of the integrals it generates. The following proposition holds true: Proposition 4.1. Let g be a semisimple Lie algebra and g 0 its subalgebra. LetM a , a ∈ 1, dim g 0 be some operators that act on the Lax operatorsL(u) by the adjoint representation of g 0 : Then (i) the operatorsM a constitute the Lie algebra isomorphic to g 0 : (ii) the Lie algebra M a centralizes the algebra of integrals generated by trL 2 (u) i.e.: M a , trL 2 (u) = 0.
Proof. Let us prove the item (i) of the proposition. For this purpose it is enough to utilize the relations (21) and Jacobi identity. Indeed, using the Jacobi identity we obtain: On the other hand, using two times the relations (21) The proof of the item (ii) is also straightforward. Indeed we have: due to the skew-symmetry of the structure tensor C cd a : C cd a = −C dc a . Proposition is proven. 2

Remark 7.
This proposition permits us to widen the algebra of integrals generated by trL 2 (u).

Generalized Gaudin Hamiltonians in external magnetic field
LetŜ (l) a , a ∈ 1, dim g, l ∈ 1, N be quantum operators that constitute a representation of the Lie algebra g ⊕N , i.e.: Let ν k , ν k = ν l , k, l ∈ 1, . . . , N be some fixed points on the complex plane belonging to the open region U in which the r-matrix r(u, v) possesses the decomposition (17). Let c(u) be a generalized shift element and ν k , k, ∈ 1, N be its regular points. In this case it is possible to show [17,19] that the following operator: satisfies the linear Lax algebra (20). Let us assume that an algebra Lie g is realized as a subalgebra of the full matrix algebra. Let us define the following set of quantum operators 3 : The following theorem holds true [17,19]: a , k ∈ 1, . . . , N a ∈ 1, dim g be as above. Then: (i) the operatorsĤ l , l ∈ 1, N constitute an abelian (commutative) algebra.
(ii) the operatorsĤ l have the following explicit form: Remark 8. The operators (24) are generalized Gaudin Hamiltonians in an external magnetic field [19]. They may be interpreted as energies of the spinning particle living at the site l, that interact with the spins living at other sites, with itself and with the external non-uniform magnetic field. The case c(u) = 0 corresponds to the generalized Gaudin Hamiltonians without external magnetic field, that were introduced in [17].

Remark 9.
Observe that additionally to the HamiltoniansĤ l there always exist the "trivial" integrals -Casimir operators, in particular, the quadratic Casimir operatorsĈ 2 l commuting with the Hamiltonians and having the form: Remark 10. Generally speaking, for higher rank Lie algebras the constructed quantum operators (24) do not constitute a complete family of "non-trivial" commutative integrals. In order to construct such a family one needs to construct the so-called "higher Gaudin Hamiltonians". This is a complicated mathematical problem solved only for the partial case g = gl(n) and standard rational r-matrices [32]. Fortunately for physical applications one needs to know only quadratic Gaudin-type Hamiltonians (24) and the "additional" integralsM a stemming from the "geometric" symmetry of the model.

Symmetries of the generalized Gaudin models in magnetic field
In this subsection we will describe the additional integralsM a for the case of the generalized Gaudin models in an external magnetic field.
The following proposition is true: Proof. The proposition is proven directly, using the invariance of the r-matrix. Indeed by direct calculation one can show that the g 0 -invariance of the r-matrix has the following expression in the component form: Using (25) we obtain: where we have used the equality dim g b=1 C d ab c b (u) = 0 which is equivalent to the condition [X a , c(u)] = 0, i.e. X a ∈ z(c(u)). Now, using the properties of the structure tensor C c ab , and invariant metric g de , it is straightforward to show that this condition is equivalent to the condition (21). The proposition is proven. 2 Remark 11. Observe that in the case of the generalized Gaudin models in an external magnetic field the symmetry of the model coincides with the intersection of the symmetry of the r-matrix and symmetry of the shift element c(u).
For obtaining of the meaningful BCS Hamiltonians starting from Gaudin-type models it will be very important that the symmetry algebra of the model contains the Cartan subalgebra h of g elements of which will play the role of the particle number operators. This is achieved, when the Cartan subalgebra belongs to a symmetry algebra of the r-matrix and the shift element c(u) takes its values in the Cartan subalgebra. In the next subsection we will consider one of the simplest examples of such the r-matrices, shift elements and Gaudin-type systems.

Example: Z 2 -graded r-matrices and Gaudin Hamiltonians
In this subsection we will consider a class of examples of classical r-matrices and generalized Gaudin models that will be basic throughout the rest of the article.
It is straightforward to show that the r-matrix (27) is g0-invariant. Due to the fact that the r-matrix (27) is not skew-symmetric, elements of X ∈ g0 are not shift elements. Nevertheless, for the reductive g0 it is possible to show [24] that shift elements in this case lie in the quotient space u −1 (g0/[g0, g0]) * .

Remark 12.
The r-matrix (26) is a particular example of σ -twisted classical r-matrices of the following form [18]: where r 12 (u − v) is a σ -invariant classical skew-symmetric r-matrix (see [18]). In more details, the r-matrix (26) is a σ -twisted rational r-matrix i.e. in this case: In order to proceed with the construction of BCS-type models we will hereafter assume that automorphism σ and the corresponding Z 2 -grading of g is such that where the subalgebra g K 0 is reductive, the subalgebras g K 0 are nilpotent and the decomposition g = g K −1 + g K 0 + g K 1 is triangular. In order for these definitions to be consistent one has to require moreover g K ±1 be abelian, i.e. if α, β ∈ ( / K ) + then α + β / ∈ + . The corresponding r-matrix (26) may be written in this case as follows: Due to the said above, in this case we may take the following shift element: c(u) = u −1 K.

Z 2 -graded Gaudin systems in a magnetic field
Let us consider the Lax operator corresponding to the r-matrix (26) and generalized Gaudin Hamiltonians in magnetic field. By virtue of the general formula (23) it has the form: The generalized Gaudin Hamiltonians in a magnetic field associated with a Z 2 -graded rmatrix are more complicated and less symmetric than usual rational ones. They have the following form: where K = dim g0 a=1 k a X a is not arbitrary element of g * 0 but belongs to the quotient space (g0/[g0, g0]) * .
In the case of Z 2 -gradings connected with triangular decompositions the Lax operator (28) is written as follows: and the Hamiltonians (29) have the following explicit form: where K = rank g i=1 k i H i coincides with the element of the Cartan subalgebra centralized by g K 0 , because in this case the quotient algebra coincides with a center of g K 0 and belongs to h. 4

Remark 13.
Observe that due to the fact that symmetry algebra of the element K coincides with the symmetry algebra of the r-matrix, introduction of the external magnetic field does not decrease the symmetry of the corresponding generalized Gaudin Hamiltonians in a magnetic field. Their symmetry algebra is the global spin algebra g K 0 . It always contains Cartan algebra as it subalgebra. For the case of the arbitrary Z 2 -graded r-matrix (26) situation is similarthe symmetry algebra of the corresponding generalized Gaudin Hamiltonians is a global spin algebra g0. But, in this more general case g0 may contain only a part of the Cartan subalgebra of g.

Generalized (p x + ip y )-type BCS models
In this section we write an explicit form of the higher rank generalized BCS Hamiltonians associated with non-skew-symmetric classical r-matrices. The form of the Hamiltonians will be different for different r-matrices, that is why we will concentrate on the main example of the present paper -Z 2 -graded r-matrices. They will produce the generalized BCS models of the p x + ip y type.

Higher rank (p x + ip y )-type BCS Hamiltonian in spin form
In order to obtain the BCS Hamiltonians of p x + ip y type it is necessary to consider their spin version. For this purpose let us consider the following linear combination of the Z 2 -graded Gaudin Hamiltonians in a magnetic field: whereĈ 2 l are Casimir operators of lth copy of g in g ⊕N or, more explicitly: where K = dim g0 a=1 k a X a belongs to the quotient algebra g0/[g0, g0]. The Hamiltonian (32) is the general integrable BCS-Richardsons Hamiltonian associated with Z 2 -grading of g and written in the spin form. In order to give it more applicable form we assume that the decomposition g = g0 + g1 is based on the triangular decomposition of g.
In this the Hamiltonian (32) acquires the following form: where K = rank g i=1 k i H i belongs to the center of the Lie algebra g K 0 . We will also need to represent the Hamiltonian (33) in the "normal-ordered form": where δ K ≡ α∈( / K ) + α. The Hamiltonian (34) is our general integrable p x + ip y BCS Hamiltonian written in the spin form.

Remark 14.
Observe also that the Hamiltonian (34) commutes with the global operatorsM i = N n=1Ŝ (k) i , i ∈ 1, rank g, representing the Cartan symmetry, which thus are also included in the constructed commutative family and are diagonalized simultaneously withĤ p x +ip y gBCS . Physically this will mean that the number of each type of fermions is a good quantum number.

Higher rank (p x + ip y )-type BCS Hamiltonian: case of gl(n)
Now we are ready to obtain the p x + ip y BCS Hamiltonian with many types of fermions. For this purpose we will use the fermionization formulas (11) and the specially chosen element K described in Section 2.2.1. Applying the general formula (34) for this choice of the algebra g and root subsystem ( / K ) + we obtain the following spin Hamiltonian: which yields (up to the constant -"vacuum energy") after the fermionization (11) the following BCS-type Hamiltonian: where we have introduced the following obvious notation: l ≡ ν −2 l .

Higher rank (p x + ip y )-type BCS Hamiltonian: case of gl(2m)
Let us restrict ourselves to the case n = 2m, p = m and assume that k 2m = −k m and the operators c l,i , c † l,i and c l,i+m , c † l,i+m correspond to the two "time-reversed" states of the fermion of number l and of type i: In such the notations we see that the Hamiltonian (36) (after division by −(k m + m)) acquires the form: where g ≡ 2(k m + m) −1 .
The first term in the Hamiltonian (37) is a kinetic term of fermions with the free energies l ≡ ν −2 l . The second describes the pairing interaction of fermions of m different types. The interaction takes place only between "time-reversed" fermions, its strength does not depend on a fermion type but does depend on the fermion number, i.e.: g k,i;l,j = g √ k √ l , where i, j ∈ 1, m, k, l ∈ 1, N.
Observe that from the general Lie-theoretical consideration exposed above it follows that there exists the "geometric" symmetry of the Hamiltonian (37) given by the Lie algebra gl(m) ⊕ gl(m) with the following generators: This Lie algebra contains the commutative subalgebra -the Cartan algebra of "number of particles" operatorsM ii ,M ss , i ∈ 1, m, s ∈ m + 1, 2m, i.e. number of particles of each type i ∈ 1, m and each "spin" δ ∈ (+, −) are good quantum numbers.

Example 5.
Let us consider the simplest one type of fermions case m = 1 (g = gl (2)). In this case we obtain the following Hamiltonian: It coincides with the Hamiltonian of the second integrable case of the reduced BCS model discovered in [21] and [25]. The "geometric" symmetry algebra of the Hamiltonian (39) is the Lie algebra gl(1) ⊕ gl(1) with the following generators: From the general consideration above the "geometric" symmetry of the Hamiltonian (39) is given by the Lie algebra gl(2) ⊕ gl (2) with the following generators: This Lie algebra contains the commutative subalgebra -the Cartan algebra of "number of particles" operatorsM ii , i ∈ 1, 4, i.e. number of protons and neutrons of each "spin" δ ∈ (+, −) are good quantum numbers.

Higher rank (p x + ip y )-type BCS Hamiltonian: case of sp(2m)
Let us obtain the (p x + ip y )-type BCS Hamiltonian with many types of fermions associated with sp(2m). For this purpose we will us use the fermionization formulas (13) and the specially chosen element K described in Section 2.2.2. The external magnetic field and the r-matrix are defined with the help of this element K and share the same symmetry with it.
Applying the general formula (34) for this choice of the algebra g, element K and root subsystem ( / K ) + , we obtain the following spin-BCS Hamiltonian: Using the fermionization formulas (13) and introducing there the notations N we obtain the following expression for the Hamiltonian (40): where l = ν −2 l , g = 2(k m + (m + 1)) −1 . In order to obtain the Hamiltonian (41) we have extracted from the Hamiltonian (40) the vacuum energy (k m + (m + 1))m N l=1 l , multiplied the Hamiltonian by −(k m + (m + 1)) −1 and used that in the considered case k 1 = · · · = k m .
The first term in the Hamiltonian (41) is a kinetic term. The second and third terms describe a pairing interaction of fermions of m different types. The interaction takes place only between "time-reversed" fermions and does not depend on their type but depends on a number of fermions of the fixed type. The interaction is less symmetric than in the considered in the previous example case of g = gl(2m). From the general consideration exposed above it follows that the symmetry algebra for this model coincides with the algebra gl(m) generated by the following operators: Let us consider the physically most interesting following example.
Example 7. Let m = 2 (g = sp (4)). In this case there are two types of fermions which may be interpreted as protons and neutrons. Introducing like in the previous example the following notations: c l,2,+ ≡ p l,+ , c † l,2,+ ≡ p † l,+ , we obtain the following integrable proton-neutron BCS Hamiltonian of the "p x + ip y " type: The Hamiltonian (42) is an exact p x + ip y analogue of proton-neutron Hamiltonian of Richardson [6]. From the general consideration above the "geometric" symmetry of the Hamiltonian (42) is given by the Lie algebra gl(2) with the following generators:

Higher rank (p x + ip y )-type BCS Hamiltonian: case of so(2m)
Let us obtain the (p x + ip y )-type BCS Hamiltonian with many types of fermions associated with so(2m). For this purpose we will us use the fermionization formulas (14) and the specially chosen element K described in Section 2.2.3.The external magnetic field and the r-matrix are defined with the help of this element K and share with it the same symmetry.
Applying the general formula (34) for this choice of the algebra g, shift element K and root subsystem ( / K ) + , we obtain the following spin-BCS Hamiltonian: Using the fermionization formulas (14), and introducing there the notations c † l,i+m ≡ c † l,i,+ , i ∈ 1, m, l ∈ 1, N we obtain the following expression for the Hamiltonian (43): where l = ν −2 l , g = 2(k m + (m − 1)) −1 and we have extracted from the Hamiltonian (43) the vacuum energy (k m + (m − 1))m N l=1 l , multiplied it by −(k m + (m − 1)) −1 , and used that in the considered case k 1 = · · · = k m .
The first term in the Hamiltonian (44) is a kinetic term of fermions. The second describes the pairing interaction of fermions of m different types. The interaction does not depend on the type of a fermion but depends on the number of fermions of the fixed type. The interaction is less symmetric than in the case of g = gl(2m). From the general consideration exposed above it follows that the symmetry algebra for this model is gl(m) and it is given by the same formulas as in the sp(2m) case.
From the above considerations it follows that the "geometric" symmetry of the Hamiltonian (45) is the Lie algebra gl (2) and it is given by the same formulas as in the considered example of sp(4).

Diagonalization
In this section we will describe the spectrum of the obtained in the previous section fermion BCS-type Hamiltonians based on non-skew-symmetric classical r-matrices in terms of solutions of the Bethe-type equations.
Let us remind that a representations of the algebra g ⊕N is called a representation with the highest weight Λ = (λ (1) , . . . , λ (N ) ), λ (l) = (λ (l) 1 , . . . , λ (l) rankg ), l ∈ 1, N if there exists such a vector Ω that: Our treatment of the problem of a diagonalization of the Gaudin Hamiltonians in an external magnetic field will be based on the following general theorem 5 : Theorem 7.1. Let g = gl(n), so(2m) or sp(2m) and the element K of the Cartan subalgebra be as the described in a Section 2. Then the spectrum of the Gaudin Hamiltonians in the external magnetic fieldĤ l (30) in the finite-dimensional representation of g ⊕N with the highest weigh (λ (1) , λ (2) , . . . , λ (N ) ) is characterized by rapidities v (i) k , k ∈ 1, M i , i ∈ 1, rank g 6 and has the following form: where is an eigenvalue ofĤ l on the highest-weight vector Ω, ρ K = α∈( K ) + α, ρ K = α∈( / K ) + α, and the rapidities v (i) k satisfy the following set of Bethe equations: andα i , i ∈ 1, rank g are simple roots.

Remark 16.
Taking into account our definition (31) of the HamiltonianĤ p x +ip y gBCS written in the spin form and the explicit form of the spectrum of the Casimir operator we obtain for its spectrum the following answer: Remark 17. Observe, that due to the special form of the elements K and ρ K , only one type of Eq. (49), namely that corresponds to the only simple root belonging to ( / K ) + , contains the summand proportional to (v (i) k ) −2 . 5 The proof of the theorem is beyond the scope of the present article. We re-direct reader to the paper [15] for details. 6 In the case of the reductive Lie algebra gl(n) one has to put everywhere in this theorem rank g → rank g − 1.

Case of gl(n)
In this section we will consider the diagonalization of the constructed BCS Hamiltonians in the case of gl(n). For simplicity we will restrict ourselves to the consideration of the physically most meaningful case n = 2m.
A space of an irreducible representation of the algebra gl(2m) ⊕N has the form: is the space of an irreducible representation of the lth copy of gl(2m) labeled by the highest weights λ (l) = (λ (l) 1 , . . . , λ (l) 2m ), l ∈ 1, N. In the representation space V there exists the highest weight vector Ω such that: We will consider the representation with the following highest weight: Observe, that in this case Ω coincides with the fermion vacuum |0 .
The following corollary of the Theorem 7.1 holds true: where "energies" E (k) i , i ∈ 1, M k , k ∈ 1, 2m − 1 satisfy the following Bethe-type equations: Proof. In order to prove the corollary it is necessary to specify the statement of the Theorem 7.1 for the Lie algebra gl(2m) and the chosen form of the highest weights λ (l) and element K.
Using the explicit form of simple roots in the case of the Lie algebra gl(2m), and identifying h * with h we may write Using the chosen form of the highest weights λ (l) and element K we obtain that (λ (l) The direct calculation gives: Using this we obtain that (ρ K ,α i ) = 2mδ im . The direct calculation also gives: Taking into account that in the case of gl(n) we have that (α i ,α j ) = 2δ ij − δ i−1,j − δ i+1,j , substituting all this into the Bethe-type equations (49), making the change of variables and taking into account that in this case g = 2(k m + m) −1 we come to the Bethe-type equations (53). Now it is left to obtain the spectrum of the BCS Hamiltonian (37). Using the formula (50) and the above calculations we obtain: .
The first summand in this term is the vacuum energy that we have to extract in accordance with the definition of the Hamiltonian (37). It is left to calculate the last term. We have: Using the mth set of the Bethe equations (53d) we obtain: .
The second summand in the right-hand side of this equality is equal to zero. Let us calculate its third and fourth summands. Making use of Eqs. (53c) and (53e) we obtain .
The first summands in the right-hand sides of these equalities are equal to zero. In order to calculate their second summands we proceed recursively using the other Bethe equations (53).
Finally we come to the Bethe equations (53a) and (53g) and obtain (because there is no E In the result we will have that Finally, taking into account that the fermion Hamiltonian (37) was defined by the division of the spin Hamiltonian (31) by −(k m + m) = − 2 g , we will obtain that the spectrum of (37) is given by the formula (52). The corollary is proven. 2 Remark 18. The spectrum of the additional integralsM ii , i ∈ 1, 2m has the same form as in the case of rational r-matrices [21] and is written as follows [15]: (54) Example 9. In the simplest m = 1 case the spectrum of the ordinary "one type of fermions" p x + ip y BCS Hamiltonian [21,22]:Ĥ p x +ip y gBCS (38) has the form: where energies E (1) i , i ∈ 1, M 1 satisfy the following Bethe-type equations: .
The spectrum m i , i ∈ 1, 2 of the additional integrals which are particle number operatorsM ii , i ∈ 1, 2 (see the formula (54)) depends only on N and M 1 and has the form: Example 10. In the important "proton-neutron" m = 2 case the spectrum of the generalized BCS HamiltonianĤ gBCS (39) has the form: where the energies E (k) i , i ∈ 1, M k , k ∈ 1, 3 satisfy the following Bethe-type equations: The spectrum m i , i ∈ 1, 4 of the additional integrals, i.e. of the particle number operatorsM ii (see the formula (54)) i ∈ 1, 4 depends only on N and M i and has the form:

Diagonalization: case of sp(2m)
In this subsection we will consider the diagonalization of the constructed BCS Hamiltonians in the case of g = sp(2m). A space of an irreducible representation of the algebra sp(2m) ⊕N has the form: is the space of an irreducible representation of the lth copy of sp(2m) labeled by the highest weights λ (l) = (λ (l) 1 , . . . , λ (l) m ), l ∈ 1, N. There exists the highest weight vector Ω such that: We will consider the representation with the following highest weight: 1 = · · · = λ (l) m = 1, ∀l ∈ 1, N. Observe that in this case Ω coincides with the fermion vacuum |0 .
where the energies E (k) i , i ∈ 1, M k , k ∈ 1, m satisfy the following Bethe-type equations: . . . . (57e) Proof. In order to prove the corollary it is necessary to specify the statement of the Theorem 7.1 for the Lie algebra sp(2m) and the chosen form of the highest weights λ (l) and element K.
Using the explicit form of simple roots in the case of Lie algebra sp(2m), and identifying h * with h we may write α i = H i andα i = H i − H i+1 , i ∈ 1, m − 1,α m = 2H m . Using the chosen form of the highest weights λ (l) and element K we obtain that (λ (l) ,α i ) = 2δ im , (K,α i ) = 2k m δ im .
Taking into account that in the case of sp(2m) we have that (α i ,α j ) = 2δ ij − δ i−1,j − δ i+1,j , i, j ∈ 1, m − 1, (α i ,α m ) = −2δ i+1,m , (α m ,α m ) = 4, substituting all this into the Bethe-type equations (49), making the change of variables and taking into account that in this case g ≡ 2(k m + m + 1) −1 we come to the Bethe equations (57). Now it is left to obtain the spectrum of the BCS Hamiltonian (41). Using the formula (50) and the calculations above we obtain: .
The first summand in this term is the vacuum energy that we have to extract in accordance with the definition of the Hamiltonian (41). It is left to calculate the last term. We have: .
The first summand in the right-hand sides of this equality is equal to zero. In order to calculate its second summand we proceed recursively using the other Bethe equations (57). Finally we come to the Bethe equations (57a) and obtain in the result (because there is no E In the result we will have that Finally, taking into account that the fermion Hamiltonian (41) was defined by the division of the spin Hamiltonian (31) by −(k m + m + 1) = − 2 g , we will obtain that the spectrum of (41) is given by the formula (56). The corollary is proven. 2 Remark 19. The spectrum of the additional integralsM i , i ∈ 1, m is written as follows [15]: where energies E (k) i , i ∈ 1, M k , k ∈ 1, 2 satisfy the following Bethe-type equations: .
The spectrum m i , i ∈ 1, 2 of the additional integrals -the number of particles operatorsM i , i ∈ 1, 2 depends only on N and M 1 , M 2 and has the form:

Diagonalization: case of so(2m)
In this subsection we will consider the diagonalization of the constructed BCS Hamiltonians in the case of g = so(2m). A space of an irreducible representation of the algebra so(2m) ⊕N has the form: V = ( N l=1 V λ (l) ), where V λ (l) is the space of an irreducible representation of the ith copy of so(2m) labeled by the highest weights λ (l) = (λ We will consider the representation with the following highest weight: λ (l) 1 = λ (l) 1 = · · · = λ (l) m = 1, ∀l ∈ 1, N. Observe, that in this case Ω coincides with the fermion vacuum |0 .
The following corollary of the Theorem 7.1 holds true: where the energies E (k) i , i ∈ 1, M k , k ∈ 1, m satisfy the following Bethe-type equations: . . .
and taking into account that in this case g ≡ 2(k m + m − 1) −1 we come to the Bethe equations (61). Now it is left to obtain the spectrum of the BCS Hamiltonian (44). Using the formula (50) and the calculations above we obtain: .
The first summand in this term is the vacuum energy that we have to extract in accordance with the definition of the Hamiltonian (44). It is left to calculate the last term. We have: In the result we will also have that

Conclusion and discussion
In the present paper we have constructed three classes of the integrable many type of fermions models associated with non-skew-symmetric classical r-matrices and the Lie algebras gl(2m), sp(2m) and so(2m). In the partial case of two types of fermions (m = 2) the obtained models are interpreted as N = Z proton-neutron integrable models with non-uniform coupling constants. We have diagonalized the constructed Hamiltonians by means of the Bethe ansatz. The important open problem is to study the solutions of the discovered Bethe-type equations in order to obtain final numerical answers for the spectrum of the constructed pairing Hamiltonians.