“Generalized” algebraic Bethe ansatz, Gaudin-type models and Zp-graded classical r-matrices

We consider quantum integrable systems associated with reductive Lie algebra gl(n) and Cartaninvariant non-skew-symmetric classical r-matrices. We show that under certain restrictions on the form of classical r-matrices “nested” or “hierarchical” Bethe ansatz usually based on a chain of subalgebras gl(n) ⊃ gl(n − 1) ⊃ ... ⊃ gl(1) is generalized onto the other chains or “hierarchies” of subalgebras. We show that among the r-matrices satisfying such the restrictions there are “twisted” or Zp-graded nonskew-symmetric classical r-matrices. We consider in detail example of the generalized Gaudin models with and without external magnetic field associated with Zp-graded non-skew-symmetric classical r-matrices and find the spectrum of the corresponding Gaudin-type hamiltonians using nested Bethe ansatz scheme and a chain of subalgebras gl(n) ⊃ gl(n − n1) ⊃ gl(n − n1 − n2) ⊃ gl(n − (n1 + ... + np−1)), where n1 + n2 + ... + np = n. © 2016 The Author. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3.

The main example of integrable models with long-range interaction are the famous Gaudin spin chains [11] associated with simple (reductive) Lie algebras g and skew-symmetric g ⊗ g-valued classical r-matrices with spectral parameters. In the papers [12,13] we have proposed a generalization of classical and quantum Gaudin models with [13] and without [12] external magnetic field associated with arbitrary non-skew-symmetric g ⊗ g-valued non-dynamical classical r-matrices with spectral parameters that satisfy the so-called generalized or "permuted" classical Yang-Baxter equation [29]. Moreover, we have shown that these models are applied in order to construct new integrable fermion models of reduced BCS-type [14][15][16] and integrable "s"-type and "p + ip"-type proton-neutron models of nuclear physics [17,18,20]. This makes a study of the generalized Gaudin models associated with non-skew-symmetric classical r-matrices important both from the mathematical and physical points of view.
The present paper is devoted to study of quantum integrable models based on classical nonskew-symmetric r-matrices and reductive Lie algebras gl(n). The important problem of the theory of integrable systems is not only the construction of quantum integrable models themselves (i.e. not only construction of the mutually commuting hamiltonian and integrals of motion) but also the development of the methods of their exact solution. For the case of quantum integrable models based on higher rank Lie algebras the main method of the exact solution is the so-called "nested" or "hierarchical" algebraic Bethe ansatz. It was invented in [21] for the quantum group case and repeated in [22] for the Lie-algebraic case in the particular case of skew-symmetric rational and trigonometric r-matrices.
In our previous paper [23] we have considered quantum integrable systems associated with the Lie algebra gl(n) and Cartan-invariant non-skew-symmetric classical r-matrices for which exists the standard procedure of the nested Bethe ansatz associated with the chain of embeddings gl(n) ⊃ gl(n − 1) ⊃ gl(n − 2) ⊃ ... ⊃ gl (1). It is necessary to observe that, contrary to Cartan-invariant skew-symmetric classical r-matrices, that are exhausted by two r-matrices, namely, rational and trigonometric ones, there are many Cartan-invariant non-skew-symmetric classical r-matrices. Moreover, contrary to Cartan-invariant skew-symmetric classical r-matrices [22], majority of Cartan-invariant non-skew-symmetric r-matrices do not support standard nested Bethe ansatz scheme associated with the chain of embeddings gl(n) ⊃ gl(n − 1) ⊃ gl(n − 2) ⊃ ... ⊃ gl (1). In fact, only three classes of non-skew-symmetric Cartan-invariant classical r-matrices support the standard nested Bethe ansatz scheme [23]. That is why, in order to make it compatible with the structure of other r-matrices, it is necessary to generalize the scheme of the nested Bethe ansatz onto other chains of embedded subalgebras.
In the present paper we generalize the "nesting" procedure of the algebraic Bethe ansatz onto the wide class of the chains or "hierarchies" of the embedded subalgebras of gl(n). The standard idea of the nested Bethe ansatz in the gl(n) case is to reduce the problem of the diagonalization of the generating function of the commutative integrals defined for the gl(n)-valued Lax matrix to the same problem for the gl(n − 1)-valued Lax matrix and then apply this method recursively [21,22]. In the present paper we propose to reduce the problem of the diagonalization of the generating function of the commutative integrals of second order for gl(n)-valued Lax matrix to the same problem for gl(n − n 1 ) ⊕ gl(n 1 )-valued Lax matrix, where n 1 ∈ 1, n − 1. It occurred to be possible under certain requirements imposed onto the considered classical r-matrix. In the cases n 1 = 1 or n 1 = n − 1 we recover our previous results [23]. Other cases are new. In such a way, proceeding recursively, i.e. applying the similar steps to the Lie algebras gl(n − n 1 ) and gl(n 1 ) we obtain a lot of possible "nestings" or "hierarchies of embeddings" for the nested Bethe ansatz. We realize this scheme to the final end using the chain of subalgebras gl(n) ⊃ gl(n − n 1 ) ⊃ gl(n − n 1 − n 2 ) ⊃ gl(n − (n 1 + ... + n p−1 )), where n 1 + n 2 + ... + n p = n, complemented by the sub-chains gl(n k ) ⊃ gl(n k − 1) ⊃ gl(n k − 2) ⊃ ... ⊃ gl(1), k ∈ 1, p. We obtain the formula for the spectrum of the quadratic generating function of the quantum integrals for the integrable models, whose representation space possess the highest weight vector and whose r-matrix supports the above chain of the embedded subalgebras. The spectrum is given in terms of "rapidities" -solutions of the corresponding Bethe-type equations. In particular, we obtain the answer for the spectrum of the Gaudin-type hamiltonians (their representation spaces do possess the highest weight vectors) and the corresponding Bethe-type equations.
All the obtained answers are quite general and are given in terms of the components of the r-matrices that satisfy the above mentioned restrictions. Nevertheless, we also consider a concrete example of the r-matrix and Gaudin-type models with and without external magnetic field that are solved by our method. The r-matrices we concentrate on are the so-called "twisted" or Z p -graded classical r-matrices [32]. The corresponding generalized Gaudin hamiltonians have been studied in detail in [24]. The considered Z p -gradings are labeled by the decomposition of n into the sum of p integers: n = n 1 + n 2 + ... + n p . The corresponding Z p -graded gl(n) ⊗ gl(n)-valued classical r-matrices after the appropriate equivalence transformation satisfy the conditions sufficient for the applicability of the nested Bethe ansatz based on the chain gl(n) ⊃ gl(n − n 1 ) ⊃ gl(n − n 1 − n 2 ) ⊃ gl(n − (n 1 + ... + n p−1 )) and complemented by subchains gl(n k ) ⊃ gl(n k − 1) ⊃ gl(n k − 2) ⊃ ... ⊃ gl(1), k ∈ 1, p. Using this fact along with the developed in this paper general theory we find the spectrum and Bethe equations for the corresponding generalized Gaudin hamiltonians with and without external magnetic field. The case p = 2 is considered in some detail.
Observe, that example of Z p -graded r-matrices is not the only example of the r-matrices that fit in the proposed scheme -some more general r-matrices from [24] also fit into it, but Z p -graded example is the simplest one for which the nesting scheme based on the chain gl(n) ⊃ gl(n − n 1 ) ⊃ gl(n − n 1 − n 2 ) ⊃ gl(n − (n 1 + ... + n p−1 )) is applicable.
At the end of the Introduction let us emphasize, that the main motivation for the proposed generalization of the nested Bethe ansatz scheme consists in the fact that the standard nested Bethe ansatz based on the chain gl(n) ⊃ gl(n − 1) ⊃ ... ⊃ gl(1) does not work for many physically interesting models, i.e. is not compatible with the structure of the corresponding r-matrices. Among such the models are Gaudin-type models in magnetic field associated with gl(n) and Z 2 -graded classical r-matrices, where Z 2 -gradings are defined by decompositions n = n 1 + n 2 , such that n 1 > 1 and n 2 > 1. In particular, in the case when n = 4, n 1 = n 2 = 2 such Gaudin-type model produce integrable "p + ip" proton-neutron model [18], whose exact solution is not provided by standard nested Bethe ansatz scheme based on the chain gl(4) ⊃ gl(3) ⊃ gl(2) ⊃ gl (1).
We would like also to outline that general non-skew-symmetric r-matrices do not satisfy usual classical Yang-Baxter equation, lie out of the Belavin-Drienfeld classification [28] and are not, generally speaking, connected with the quantum groups or related structures. That is why it is methodologically important to develop a theory of the corresponding quantum integrable systems independently of much better elaborated quantum-group formalism [33,34].
The structure of the present paper is the following: in the second section we remind the construction of the quantum integrable systems based on the classical r-matrices. In the third (main) section we develop a theory of the nested Bethe ansatz based on non-standard chains of subalgebras gl(n) ⊃ gl(n − n 1 ) ⊃ gl(n − n 1 − n 2 ) ⊃ gl(n − (n 1 + ... + n p−1 )). In the fourth section we consider example of Z p -graded classical r-matrix and the corresponding Gaudin models with and without external magnetic field.

Definitions and notations
Let g = gl(n) be the Lie algebra of the general linear group over the field of complex numbers. Let X ij , i, j = 1, n be a standard basis in gl(n) with the commutation relations: (1) with values in the tensor square of the algebra g = gl(n) is called a classical r-matrix if it satisfies the following generalized classical Yang-Baxter equation [29][30][31][32]: where 1 ⊗ X kl , etc. and r ij,kl (u, v) are matrix elements of the r-matrix r (u, v).

Remark 1.
In the case of skew-symmetric r-matrices, i.e. when r 12 (u 1 , u 2 ) = −r 21 (u 2 , u 1 ), where r 21 (u 2 , u 1 ) = P 12 r 12 (u 1 , u 2 )P 12 and P 12 interchanges the first and second spaces in tensor product, the generalized classical Yang-Baxter equation reduces to the usual classical Yang-Baxter equation [35,28,36]: Observe, that gauge transformations and re-parametrizations are equivalences in the spaces of solutions of the equation (3), while "rescaling" is not, because it does not preserve skew-symmetry property of the r-matrix.
In the present paper we are interested only in the r-matrices that by re-parametrization and rescaling may be brought to the form possessing the following decomposition: where = n i,j =1 X ij ⊗ X ji is the tensor Casimir and r 0 (u 1 , u 2 ) is a regular on the "diagonal" For the subsequent we will also need the following three definitions: Definition 2. We will call classical r-matrix to be Observe, that on the level of Lie groups this definition means exactly the G 0 -invariance: where g ∈ G 0 and G 0 is a Lie group of the algebra g 0 . Definition 3. We will call classical gl(n)-valued r-matrix to be diagonal in the root basis if it has the following form: Diagonal in the root basis classical r-matrices are automatically Cartan-invariant. Inverse is not true: not all Cartan-invariant r-matrices are diagonal in the root basis.

Remark 2.
Observe that in the case of the diagonal r-matrices the generalized classical Yang-Baxter equation (2), is re-written in the component form as follows: where i, j, l ∈ 1, n and all three indices cannot coincide simultaneously. c ij (u)X ij of one complex variable is called "generalized shift element" if it satisfies the following equation: In the subsequent we will be interested only in the diagonal shift elements, i.e. shift elements of the following form: Remark 3. Observe, that for skew-symmetric classical r-matrices any element of its symmetry algebra is a shift element. For non-skew-symmetric r-matrices it is not true.

Algebra of Lax operators and classical r-matrices
Using the classical r-matrix r(u 1 , u 2 ) it is possible to define the following "tensor" Lie bracket: where ij (u)X ij , r 21 (u 2 , u 1 ) = P 12 r 12 (u 1 , u 2 )P 12 .
The tensor bracket (9) between the quantum Lax matrices L (u 1 ) and L (u 2 ) is a symbolical expression of the Poisson brackets between their matrix elements. In the case of the diagonal r-matrices (6) which are the main object of interest in the present article the corresponding commutation relations (9) have the following simple form: The explicit form of the operators L ij (u) as functions of the auxiliary spectral parameter u is not arbitrary. It agrees with the structure of the r-matrix r(u 1 , u 2 ) and depends on the concrete quantum system under consideration. In the present paper we will consider the simplest but the most important examples of such the systems.
Meantime, let us clarify the role of the shift elements in the algebra of Lax operators. The following proposition holds true: be the Lax operator satisfying the commutation relations (9) and c(u) be the shift element satisfying equation (8). Then the operator also satisfies the commutation relations (9).
(Proof of this proposition follows from the explicit form of commutation relations (9), definition (8) and the fact that c ij (u) are c-numbers.)

Quantum integrals
In this subsection we will explain the connection of classical non-skew-symmetric gl(n) ⊗ gl(n)-valued r-matrices with quantum integrability. It will be shown that, just like in the case of classical r-matrix Lie-Poisson brackets [29,31], the Lie bracket (9) leads to an algebra of mutually commuting quantum integrals.
Let us consider the following functions in generators of the Lax algebra: The following important theorem holds true [26]: be the Lax operator satisfying the commutation relations (9) with the diagonal in the root basis classical r-matrix. Assume that in some open region U × U ⊂ C 2 the function r(u, v) possesses the decomposition (4). Then the operator-valued functions τ n (u) and t n (u) are a generators of a commutative algebra, i.e.: (The proof of the theorem involves the Leibnitz rule and Jacobi identity for the commutator, the gl(n)-invariance of the corresponding quadratic form and some consequences of the generalized classical Yang-Baxter equations.)

Remark 4.
It is necessary to emphasize that, generally speaking, operator does not belong to a center of the algebra of Lax operators. Indeed, even in the simplest case of the diagonal r-matrices from the commutation relations (10) we obtain: This expression is not zero if r kk (u, v) = r ll (u, v).
From Theorem 2.1 follows the next corollary:

Remark 5.
Observe that the operators Ĥ ν l are exact analogs of the generalized Gaudin hamiltonians [12] and coincide with them for the special choice of the Lax operator. From the equality (12) it is easy to deduce that the operators Ĉ ν l belong to the center of the algebra of Lax operators.
To summarize: in this section we have constructed an algebra of quantum commutative operators that coincide with a quantization of linear and quadratic subalgebra of the algebra of Lie-Poisson commuting integrals of a classical integrable system admitting Lax representation. The problem of quantization of the other "higher" integrals is complicated. It is solved only in the partial case of the classical rational r-matrices (i.e. r-matrices for which r 0 (u, v) ≡ 0) in [37]. Fortunately, for the physically applications the most important are quadratic integrals. Moreover, as we will show below, their diagonalization does not depend on the higher integrals. That is why we do not consider the problem of the quantization of the higher integrals here.

Example: generalized Gaudin systems
Let us now consider the most important for the applications (and the simplest at the same time) examples of quantum integrable systems associated with classical r-matrices.
Let Ŝ (m) ij , i, j = 1, n, m = 1, N be linear operators in some Hilbert space that span Lie algebra isomorphic to gl(n) ⊕N with the commutation relations: Let us fix N distinct points of the complex plain ν m , m = 1, 2, ..., N . It is possible to introduce the following quantum Lax operator [12]: Using generalized classical Yang-Baxter equation it is possible to show that it satisfies a linear r-matrix algebra (9). This quantum Lax operator is a Lax operator of the generalized gl(n)-valued Gaudin spin chains. For the applications important are also the shifted Lax operators [13]: This quantum Lax operator is a Lax operator of the generalized gl(n)-valued Gaudin spin chains in an external magnetic field and a generalized shift element of an external nonhomogeneous magnetic field. Hereafter we will consider only diagonal r-matrices and diagonal shift elements.
As it was noted above, residues of generating function t n (u) produces only trivial integrals ii . Residues of the second order generating function τ (u) produces nontrivial integrals Ĥ ν l : (15) in the case of the "unshifted" Lax operators and ii (16) in the case of "shifted" Lax operators. The hamiltonians (15) are the generalized Gaudin hamiltonians corresponding to the diagonal gl(n) ⊗ gl(n)-valued r-matrix and hamiltonians (16) are the generalized Gaudin hamiltonians in external magnetic field corresponding to the same diagonal gl(n) ⊗ gl(n)-valued classical r-matrix.

Hierarchical Bethe ansatz
The standard idea of the nested or hierarchical algebraic Bethe ansatz [21,22] in the gl(n) case is to reduce the problem of the diagonalization of the generating function of the commutative integrals of gl(n) to the same problem for the Lie algebra gl(n − 1) and then apply this method recursively, using the hierarchy of subalgebras gl(n) ⊃ gl(n − 1) ⊃ ... ⊃ gl (1).
In this section we will show how to generalize the nested Bethe ansatz onto the wide class of hierarchies of subalgebras starting from the reduction of the problem of diagonalization of the generating function of the commutative integrals on gl(n) to the same problem on gl(n 1 ) ⊕ gl(n − n 1 ). This procedure allows to solve by this method quantum integrable models associated with the r-matrices not compatible with the usual sequence of embeddings gl(n) ⊃ gl(n − 1) ⊃ ... ⊃ gl(1).

Reduction to subalgebra
Let V be a space of an irreducible representation of the algebra of Lax operators. Let us assume that there exist a vacuum vector ∈ V such that: (17) and the whole space V is generated by the action of L kl (u), k < l on the vector . In order to diagonalize the generating function τ(u) and construct the corresponding set of eigen-vectors we will reduce this problem to the same problem on the subalgebra gl(n 1 ) ⊕ gl(n − n 1 ). We will do this in several steps.

Reduction of the problem to subalgebra on "vacuum" subspace
Let us at first consider a subspace V 0 ⊂ V consisting of the vectors v such that: where n 1 ∈ 1, n is some fixed integer.
Using the commutation relations in the algebra of the Lax operators it is easy to see that this subspace is invariant with respect to the action of the subalgebra of Lax operators taking values in subalgebra gl(n 1 ) ⊕ gl(n − n 1 ).
Let us return to the generating function τ n (u). We have: Let us take into account that τ n 1 (u) are the generating function of commutative integrals of a gl(n 1 )-valued and gl(n − n 1 )-valued subalgebras of the algebra of Lax operators and: . (This equality is obtained by taking the limit u 2 → u 1 = u in the corresponding commutation relations of the Lax algebra.) Let us assume that the following conditions on the r-matrix are satisfied: In this case we will have: In the result we obtain: where v ∈ V 0 . That is the problem of diagonalization of τ n (u) on this subspace is reduced to the problem of the simultaneous diagonalization of the linear and quadratic generating functions of the subalgebras gl(n 1 ) and gl(n − n 1 ). So we have performed the needed reduction on the space V 0 . Using this result, in the next sections we will perform it on the whole space V .

Preparation for the reduction of the problem to the subalgebra
In order to construct the needed Bethe vectors and perform the required reduction to subalgebra we need to introduce some notations and perform many auxiliary calculations.
By the direct calculations, using the commutation relations in the Lax algebra (9) and definition of the function τ n (u) one can prove the following formula: which will be used by us in the subsequent calculations. Let us introduce the following operator-valued matrix: T n (u) =τ n (u)1 n .
We will also use the following operator-valued matrices (parts of the Lax matrix): ij is a standard basis element of gl(n) acting in some auxiliary linear space C n carrying additional index 0. Hereafter it will be convenient for us to label auxiliary linear subspace by the indices 0 and 0 (instead of the indices 1 and 2 that were used for this purpose before).
Observe that the matrix Â (u) belonging to gl(n 1 ) may act (as a matrix) on B (u) only from the left, while the matrix D (u) belonging to gl(n − n 1 ) may act (as a matrix) on the matrix B (u) only from the right. On the other hand in order to find the spectrum we will need to write the components of Â (u) to the right of operators B (u). That is why we need to introduce the following symbolical notation for the multiplication of two operator-valued matrices Ŷ (u) = Now we can re-write relation (20) in the matrix form. The following proposition holds: Proposition 3.1. Let the r-matrix r(u, v) satisfy the conditions (18) then the commutation relations (20) are written as follows: where , 0 is a pairing in a second tensor space labeled by the index 0 and: Proof. In order to prove this proposition it is necessary to re-write the formula (20) in the following form: Now, using the very definitions of the matrices Â (u) and D (u), multiplying the both sides of (22) by X (0) kl , summing the both sides over k ∈ 1, n 1 , l ∈ n 1 + 1, n and taking into account that by our assumption on the form of the r-matrix: r ji (u, v) = r −,n 1 (u, v), ∀i ∈ 1, n 1 , j ∈ n 1 + 1, n we obtain from (22) that: where we have introduced the following notation: Let us now calculate this expression. Using the commutation relation of the Lax algebra (10) and re-grouping the summands we obtain: Now, using the fact that by our assumption on the form of the r-matrix: r ji (u, v) = r −,n 1 (u, v), r ij (u, v) = r +,n 1 (u, v) ∀i ∈ 1, n 1 , j ∈ n 1 + 1, n we obtain that the first two sums in the righthand-side of (24) is equal to the expression 2B(v)ρ LR (v, u). That is why it is left to transform the last two sums of the expression (24). For this purpose it is necessary to use the following identities: (v, u), which are the differential consequences of the generalized classical Yang-Baxter equation written in the component form (7).
Using this and the condition r ji (u, v) = r −,n 1 (u, v), r ij (u, v) = r +,n 1 (u, v) ∀i ∈ 1, n 1 , j ∈ n 1 + 1, n we finally obtain that the last sums in the left-hand-side of the expression (24) are equal to the expression r −,n 1

(u, v)B(u)ρ L,R (v). Proposition is proven. 2
Let us re-write the obtained in the above proposition formula in the different form, which will be important for the subsequent. For this purpose we need to introduce some new operators: It is easy to see that the operators D (1) (u) and Â (1) (u) are the gl(n − n 1 ) and gl(n 1 )-valued Lax operators, satisfying the Lax algebra (9) were the r-matrix is the corresponding gl(n − n 1 ) ⊗ gl(n − n 1 ) and gl(n 1 ) ⊗ gl(n 1 )-valued sub r-matrices of the initial gl(n) ⊗ gl(n)-valued r-matrix r(u, v). Let, furthermore, The following corollary of Proposition 3.1 holds true: Corollary 3.1. The equality (21) are written in the new notations as follows: Proof. To prove this corollary it is enough to observe that the expressions for ρ LR (v, u) and ρ LR 0 (v) are re-written in the following form: Then calculating explicitly 1 2 tr n−n 1 (D (1) (u)) 2 and 1 2 tr n 1 (Â (1) (u)) 2 and their residues and comparing them with the right-hand-side of (21) one obtains the equality (25). 2 For the subsequent we will need to prove similar proposition and corollary for the tensor products of many B i (v i ). For this purpose we will need to "lift" all the above matrices with the non-commuting entries to the tensor products of the M copies of gl(n) and we will use the following "tensorial" notations: B k (v k ) ≡ 1 n ⊗ · · · ⊗B(v k ) ⊗ · · · ⊗ 1 n , Â k (v k ) ≡ 1 n ⊗ · · · ⊗ A(v k ) ⊗ · · · ⊗ 1 n , D k (v k ) ≡ 1 n ⊗ · · · ⊗D(v k ) ⊗ · · · ⊗ 1 n , where Â (v k ), B (v k ), D (v k ) stand in the k-th component of the tensor product and 1 n is the unit matrix in the space C n . We will also "lift" the matrix τ n (u) to this tensor product in the following way: T n,M (u) ≡τ n (u)1 n ⊗ · · · ⊗ 1 n ⊗ · · · ⊗ 1 n .
The most important technical step we need to accomplish is to calculate the commutator The following important lemma holds: Lemma 3.1. The following commutation relations are valid: where "check" over the operator B i (v i ) means that it is absent in the product, τ n−n 1 (u) = 1 2 tr n−n 1 (D (1) (u)) 2 where Â (1) (u) and D (1) (u) are defined as follows:D (1) Idea of the proof. The lemma is proven in the direct way, using the Leibnitz rule for the commutator, Proposition 3.1, Corollary 3.1 and the classical Yang-Baxter equation written in the component form. We will not expose the detailed proof here, due to its long and technical character. 2

Reduction of the problem to subalgebra on the Bethe vectors
Now we are ready to reduce a problem of the diagonalization of the generating function of quadratic commutative quantum integrals on gl(n) to the same problem on the subalgebra gl(n 1 ) ⊕ gl(n − n 1 ).
Using the action formula (19) and Lemma 3.1 we obtain the following theorem: Theorem 3.1. Let the r-matrix r(u, v) be such that the conditions (18) are satisfied. Let, moreover, the r-matrix satisfy the following conditions: (u, v), ∀a, b ∈ 1, n 1 and r kk (u, v) = r ll (u, v), ∀k, l ∈ n 1 + 1, n (28) or (u, v), ∀k, l ∈ n 1 + 1, n and (r 0 −,n 1 (u, u) + r 0 +,n 1 (u, u)) = 0. (29) Let v i 1 (Here e i are basis vectors in the space C n -vector-columns with unit on the place i and zeros elsewhere and e * i are dual basis vectors in the space C n -vector-rows with unit on the place i and zeros elsewhere.) Then the eigen-values n (u) and (1) n−n 1 (u) of the operators τ n (u) and τ (1) n 1 (u), τ (1) n−n 1 (u) are connected as follows: where c n 1 (u) and c n−n 1 (u) is a spectrum of the linear Casimir operators t n 1 (u) and t n−n 1 (u) of Lax subalgebras with the values in subalgebras gl(n 1 ) and gl(n − n 1 ) in the representations of these Lax algebras with the highest weights ( 11 (u), ...., n 1 n 1 (u)) and ( n 1 +1n 1 +1 (u), ...., nn (u)) correspondingly.
Proof. In order to prove the theorem, let us calculate the action T n,M (u) · V(v 1 , ..., v M ), where the Bethe vector V(v 1 , ..., v M ) is defined as follows: We will have: Using the action formula (19) and the Lemma 3.1 we obtain: Now, in order to prove this theorem it is left to identify the action of the operators τ n (u) with the action of the operator T n,M (u) in the corresponding spaces. Let us consider the spaces V ⊗ (C n * ) ⊗M ⊗ (C n ) ⊗M . There is a canonical projection of this space onto V : V ⊗ (C n * ) ⊗M ⊗ (C n ) ⊗M → V made with the help of the pairing , :(C n * ) ⊗M ⊗ (C n ) ⊗M → C which is obtained as a prolongation of the pairing C n * × C n → C onto the tensor product of M copies of these spaces. It is easy to show that eigen-vectors of Taking into account all the above and the fact that under the condition (28) the operators t n 1 (u) and t n−n 1 (u) are the Casimir operators of Lax subalgebras with the values in subalgebras gl(n 1 ) and gl(n − n 1 ) and are proportional to unit operators, or, due to the condition (29) the operators ∂ utn 1 (u) and ∂ utn−n 1 (u) are the Casimir operators of Lax subalgebras with the values in subalgebras gl(n 1 ) and gl(n − n 1 ) and are proportional to unit operators, we obtain the statement of the theorem as a direct consequence of the formula (32).
The theorem is proven. 2

Diagonalization of τ (u)
In order to diagonalize quadratic generating function of the commutative quantum integrals it is necessary to apply in a recursive manner Theorem 3.1. For this purpose it is necessary to fix the chain of embedded subalgebras which is used in the recursion. We will use the following chain of the subalgebras: gl(n) ⊃ gl(n − n 1 ) ⊃ gl(n − n 1 − n 2 ) ⊃ gl(n − (n 1 + ... + n p−1 )), where n 1 + n 2 + ... + n p = n, complemented by the following sub-chains: gl(n k ) ⊃ gl(n k − 1) ⊃ ... ⊃ gl(1), k ∈ 1, p. The choice of the chain of the embedded subalgebras for the recursion in the nested Bethe ansatz will impose certain conditions on the corresponding classical r-matrices.
We will obtain the following formula for the spectrum of the generating function: Here n−(n 1 +..+n k ) (u). By the very definition these are eigen-values of the linear Casimir functions of the Lie subalgebras gl(n k ) and gl(n − (n 1 + .. + n k )) correspondingly. They have the following form: It is necessary only to find (k−1) ii (u), keeping in mind the procedure described in Lemma 3.1 that had been applied to the Lax operator k − 1 times. We will have: In order to see this it is necessary to take into account that, by the very procedure described in Lemma 3.1 and by the k − 1 steps based on the chain of subalgebras gl(n) ⊃ gl(n − n 1 ) ⊃ gl(n − n 1 − n 2 ) ⊃ ... ⊃ gl(n − (n 1 + ... + n k−1 )), we will have: c kk (ν m ) are the components of the shift element -external magnetic field and the "ra- g j be Z p = Z/pZ grading of g, i.e.: [g i , g j ] ⊂ g i+j where j denotes the class of equivalence of the elements j ∈ Z modp Z. It is known, that Z p -grading of g may be defined with the help of some automorphism σ of the order p, such that σ (g i ) = e 2πik/p g i and g 0 is the algebra of σ -invariants: σ (g 0 ) = g 0 . It is also known [32] (see also [24]), that using this data it is possible to define the following classical r-matrix: where X α,j is a basis in the space g j , X α,−j is a dual basis in the space g −j and (j) 12 = dim g j α=1 g αβ j X α,j ⊗ X β,−j is a projection operator onto the subspace g j , g αβ j ≡ (X α,−j , X β,j ). In particular, (0) is the tensor Casimir of the subalgebra g 0 .

Z p -graded classical r-matrices and generalized Gaudin models
For Z p -graded r-matrix (46) corresponding to g = gl(n), the specified Z p -grading and the gauge (47) the Gaudin-type Lax operator is written as follows: The Gaudin-type Lax operator with an external "diagonal" magnetic field has the form: The corresponding Gaudin-type hamiltonians (15) and Gaudin-type hamiltonians in an external magnetic field (16) are written as follows: In the next subsection we will consider the diagonalization of these hamiltonians by means of the nested Bethe ansatz, according to new "nesting" scheme described in the previous sections. But at first let us consider in some detail the simplest non-banal (i.e. different from the standard rational r-matrix case) example of the described in this subsection Lax operators and hamiltonians. .
Observe also, that the case p = 2 can be treated in the framework of our scheme directly, without the application of the gauge transform to the r-matrix (49). The Bethe equations obtained using this direct treatment are equivalent to those described above (see [27] where integrable spin, boson and spin-boson models corresponding to the case p = 2 have been considered). They are also equivalent to the Bethe equations obtained in the case p = 2 using the analytical Bethe ansatz [19].

Conclusion and discussion
In the present paper we have generalized nested Bethe ansatz onto a wide class of chains or "hierarchies" of embedded subalgebras of Lie algebra gl(n) and certain sub-classes of the Cartan-invariant non-skew-symmetric classical r-matrices. We have shown that among such the r-matrices there are "twisted" or Z p -graded non-skew-symmetric classical r-matrices with spectral parameters. We have considered an example of the corresponding generalized Gaudin models with and without external magnetic field and found the spectrum of their hamiltonians using nested Bethe ansatz.
The important on-going problem is a solution (at least a numerical one) of the obtained nested Bethe equations, which is necessary for concrete applications of the corresponding integrable models in quantum physics.
Finally we would like to suggest that some of the results of the present paper on the generalized nested algebraical Bethe ansatz for the case of linear Lax algebras can be prolonged onto the quadratic "quantum-group" cases. In particular, we suppose that for the case of the Reflection Equation Algebras corresponding to the considered "Z 2 -graded" classical r-matrices [19], there should be a generalization of the nested algebraic Bethe ansatz scheme onto a chain of subalgebras gl(n) ⊃ gl(n 1 ) or gl(n) ⊃ gl(n 2 ), n 1 + n 2 = n, n 1 > 1, n 2 > 1, complimented by sub-chains gl(n k ) ⊃ gl(n k − 1) ⊃ gl(n k − 2) ⊃ ... ⊃ gl(1), k ∈ 1, 2.