Separation of variables bases for integrable $gl_{\mathcal{M}|\mathcal{N}}$ and Hubbard models

We construct quantum Separation of Variables (SoV) bases for both the fundamental inhomogeneous $gl_{\mathcal{M}|\mathcal{N}}$ supersymmetric integrable models and for the inhomogeneous Hubbard model both defined with quasi-periodic twisted boundary conditions given by twist matrices having simple spectrum. The SoV bases are obtained by using the integrable structure of these quantum models,i.e. the associated commuting transfer matrices, following the general scheme introduced in [1]; namely, they are given by set of states generated by the multiple actions of the transfer matrices on generic co-vectors. The existence of such SoV bases implies that the corresponding transfer matrices have non degenerate spectrum and that they are diagonalizable with simple spectrum if the twist matrices defining the quasi-periodic boundary conditions have that property. Moreover, in these SoV bases the resolution of the transfer matrix eigenvalue problem leads to the resolution of the full spectral problem, i.e. both eigenvalues and eigenvectors. Indeed, to any eigenvalue is associated the unique (up to a trivial overall normalization) eigenvector whose wave-function in the SoV bases is factorized into products of the corresponding transfer matrix eigenvalue computed on the spectrum of the separate variables. As an application, we characterize completely the transfer matrix spectrum in our SoV framework for the fundamental $gl_{1|2}$ supersymmetric integrable model associated to a special class of twist matrices. From these results we also prove the completeness of the Bethe Ansatz for that case. The complete solution of the spectral problem for fundamental inhomogeneous $gl_{\mathcal{M}|\mathcal{N}}$ supersymmetric integrable models and for the inhomogeneous Hubbard model under the general twisted boundary conditions will be addressed in a future publication.

Integrable quantum models define the natural background to look for exact non-perturbative results toward the complete solution of some 1+1 dimensional quantum field theories or some equivalent two dimensional systems in statistical mechanics. They have found natural applications in the exact description of several important phenomena in condensed matter and have provided exact results to be compared with experiments. A prominent example is the quantum Heisenberg spin chain [60] introduced as a model to study phase transitions and critical points of magnetic systems. First exact results for the Hamiltonian's spectrum (eigenvalues and eigenvectors) have been obtained by Bethe for the spin 1/2 XXX chain, thanks to his famous coordinate ansatz [61]. Then it has been extended to the anisotropic spin 1/2 XXZ chain in [62,63], while Baxter [64,65] has obtained first exact results for the Hamiltonian's spectrum of the fully anisotropic spin 1/2 XYZ chain. In statistical mechanics, these quantum models correspond to the six-vertex and eight-vertex models. The icetype (six-vertex) models [66], accounting for the residual entropy of water ice for crystal lattices with hydrogen bonds, has first been described in the Bethe Ansatz framework in [67]. Instead for the eight-vertex model [68,69] first exact results in the two-dimensional square lattice are due to Baxter [65,70,71]. The integrable structure 1 of these spin chains and statistical mechanics models has been revealed in the Baxter's papers [64,71,77]. There, the one parameter family of eight-vertex transfer matrices have been shown to be commutative and the XYZ Hamiltonian to be proportional to its logarithmic derivative, when computed at a particular value of its spectral parameter. The subsequent development of a systematic description of quantum integrable models has been achieved through the development of the quantum inverse scattering framework [10][11][12][13][14][15][16][17][18]. In particular, the work of Faddeev, Sklyanin and Takhtajan [11] has set the basis for the classification of the Yang-Baxter algebra representations and the natural framework for the discovering of quantum groups [78][79][80][81]. The paper [11] has also introduced the Algebraic Bethe Ansatz (ABA), an algebraic version of the original coordinate Bethe Ansatz.
In this context, the form factor expansion has proven to be a very powerful tool: on the one hand, to compute dynamical structure factors [95,96], quantities directly accessible experimentally through neutron scattering [97]; on the other hand, to have access to the asymptotic behavior of correlation functions of these XXZ chains in the thermodynamic limit and explicit contact with conformal field theory [89][90][91][98][99][100][101][102][103][104].
Integrable quantum models also led to non-perturbative results in the out-of-equilibrium physics context, ranging from the relaxation behavior of some classical stochastic processes to quantum transport. The XXZ quantum spin chains, under general integrable boundary conditions, appear for example both in the description of the asymmetric simple exclusion processes [105][106][107][108][109][110][111] and the description of transport properties of quantum spin systems [112,113].
The new experiments allowing ultra-cold atoms to be trapped in optical lattices have produced concrete realizations of quantum integrable lattices, like the Heisenberg spin chains but also more sophisticated models like the Hubbard model [114,115]. They provide a further natural context for direct comparison of exact theoretical predictions with experiments.
The Hubbard model is of fundamental importance in physics. It is a celebrated quantum model in condensed matter theory, defining a first generalization beyond the band approach for modeling the solid state physics. It manages to describe interacting electrons in narrow energy bands and allows to account for important physical phenomena of different physical systems. Relevant examples are the high temperature superconductivity, band magnetism and the metal-insulator transitions. We refer to the review book [9] for a more detailed description of its physical applications and of the known exact results and relevant literature. Here, let us recall that the Hubbard chain is integrable in the quantum inverse scattering framework and it has been first analyzed by Bethe Ansatz techniques in a famous paper by Lieb and Wu [116,117]. Coordinate Bethe Ansatz wave-functions for the Hubbard Hamiltonian eigenvectors have been obtained in [118] and subsequent papers. The quantum inverse scattering formulation has been achieved thanks to the Shastry's derivation of the R-matrix [119][120][121]. From this, the one parameter family of commuting transfer matrices, generating the Hubbard Hamiltonian by standard logarithmic derivative, can be introduced. The proof that this R-matrix satisfies the Yang-Baxter equation has been given in [122]. In [123][124][125] a Nested Algebraic Bethe Ansatz for the Hubbard model has been introduced while the quantum transfer matrix approach to study the thermodynamics of the Hubbard model has been considered in [126]. Interestingly, the Hubbard Hamiltonian is invariant under the direct sum of two Y (sl 2 ) Yangians, as derived in [127][128][129][130], while the structure of the fusion relations for the Hubbard model have been studied in [131]. It is also relevant to remark that under a strong coupling limit and some special choice of the remaining parameters [9], the Hubbard Hamiltonian leads to the Hamiltonian of the t-J gl 1|2 -supersymmetric model, another well known model for the description of the high-temperature superconductivity, see [132][133][134] and references therein.
While integrable quantum models naturally emerge in the description of 1+1 quantum or 2dimensional statistical mechanics phenomena, they are not really confined to this realm. For example, they play a fundamental role in deriving exact results also for four-dimensional quantum field theory like the planar N = 4 Supersymmetric Yang-Mills (SYM) gauge theory, see the review paper [135] and references therein. In this context, integrability tools have been used to derive exact results, like characterizations of the scaling dimensions of local operators for general values of the coupling constant. Notably, such results can be used also as a test of the AdS/CFT correspondence in this planar limit. Indeed, holding for arbitrary values of the coupling, they allow for a verification of the agreement both at weak and strong couplings with the perturbative results obtained respectively in gauge and string theory contexts. Integrability is also becoming relevant in the exact computation of observables. Interestingly, quantum integrable higher rank spin and super-spin chains have found applications in the computation of correlation functions in N = 4 SYM, see e.g. [136][137][138][139][140][141]. The same integrable Hubbard model enters in the description of the planar N = 4 SYM gauge theory [135,142]. Indeed, relevant examples are the connection between its dilatation generator and the Hubbard Hamiltonian derived in [143] and the equivalence between the bound state S-matrix for AdS 5 × S 5 superstring [144][145][146] and the Shastry's R-matrix of the Hubbard model shown in [147,148]. Moreover, this S-matrix enjoys the Yangian symmetry associated to the centrally extended su(2|2) superalgebra [144][145][146] and an Analytic Bethe Ansatz description of the spectrum has been introduced on this basis in [147,148].
The large spectrum of applications of these higher rank quantum integrable models and of the Hubbard model, clearly motivate our interest in their analysis by quantum separation of variables.
We have already achieved such program in our previous publications [1,[149][150][151] for the spectrum of the higher rank quantum integrable models associated to fundamental representations of the Y (gl n ) and U q (ĝl n ) Yang-Baxter algebra and of the Y (gl n ) reflection algebra. While in [152], we have described in detail how our approach works beyond fundamental representations 2 for Y (gl 2 ).
Here, we construct the quantum Separation of Variables (SoV) bases in the representation spaces of both the fundamental inhomogeneous gl M|N Yang-Baxter superalgebras and the inhomogeneous Hubbard model under general quasi-periodic twisted boundary conditions defined by twist matrices having simple spectrum. The SoV bases are constructed by using the known integrable structure of these quantum models, i.e. the associated commuting transfer matrices, following our general ideas introduced in [1]. The SoV bases are generated by the multiple actions of the transfer matrices on a generic co-vector of the Hilbert space. The fact that we are able to prove that such sets of co-vectors indeed form bases of the space of states implies important consequences on the spectrum of the transfer matrices. In fact, it follows that the transfer matrices have non degenerate (simple) spectrum or that they are diagonalizable with simple spectrum if the twist matrix respectively has simple spectrum or is diagonalizable with simple spectrum. Moreover, in our SoV bases the resolution of the transfer matrix eigenvalue problem is equivalent to the resolution of the full transfer matrix spectrum (eigenvalues and eigenvectors). Indeed, our SoV bases allows us to associate uniquely to any eigenvalue an eigenvector whose wave-function has the factorized form in terms of product of the transfer matrix eigenvalues on the spectrum of the separate variables.
It is worth pointing out that for these classes of higher rank quantum integrable models, fewer exact results are available when compared to those described for the best known examples of the XXZ spin 1/2 quantum integrable chains. Exact results are mainly confined to the spectral problem and only recently some breakthrough has been achieved toward the dynamics in the framework of the Nested Algebraic Bethe Ansatz (NABA) [155][156][157][158], for some higher rank spin and super-spin chains [159][160][161][162][163][164][165][166][167][168].
More in detail, in the supersymmetric case, the associated spectral problem has been analyzed by using the transfer matrix functional relations, generated by fusion [169][170][171] of irreducible representations in the auxiliary space of the representation. The Analytic Bethe Ansatz [171][172][173] developed in this functional framework has been applied to the spectral problem of these supersymmetric models. An important step in the systematic description and analysis of these functional equations has been done by rewriting them in Bazhanov and Reshetikhin's determinant form in [174], and in a Hirota bilinear difference equation form in [175][176][177]. These so-called T -systems appear both in classical and quantum integrability. An interesting account for their relevance and different application areas can be found in [178]. The validity of these fusion rules and of Analytic Bethe Ansatz description in the supersymmetric case have been derived in [179][180][181][182]. In [183,184] a method has been developed and applied to the supersymmetric case based on the use of Bäcklund transformations on the Hirota-type functional equations [185][186][187]. It allows a systematic classification of the different Nested Algebraic Bethe Ansatz equations and T Q-functional equations, which emerge naturally in the supersymmetric case, due to different possible choices of the systems of simple roots. It also allows the identification of QQ-functional equations of Hirota type for the Baxter's Q-functions. Nested Algebraic Bethe Ansatz [158] has been successfully used to get Bethe vectors representations for fundamental representations of Y (gl M|N ) and U q (ĝl M|N ), see also the recent result [188], while determinant formulae for Bethe eigenvector norms, scalar products and some computations of form factors have been made accessible in [189][190][191][192] for the Y (gl 1|2 ) and Y (gl 2|1 ) case.
The development of the quantum separation of variables method for these supersymmetric integrable quantum models has the advantage to define a natural framework where to prove the completeness of the spectrum description. A simple algebraic representation for the transfer matrix eigenvectors is gained and determinant formulae for scalar products should be accessible. Indeed, completeness has been shown to be a built-in feature of separation of variables verified for a large class of quantum integrable models [38][39][40][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58]. The SoV representation of eigenvectors by factorized wave-functions in terms of the eigenvalues of the transfer matrix, or of the Baxter's Q-operator [65, is indeed extremely simple and universal. It also brings to Algebraic Bethe Ansatz rewriting of non-Nested type for the eigenvectors, i.e. as the action on a SoV induced "reference vector" of a single monomial of SoV induced "B-operators" over the zeros of an associated Q-operator eigenvalue [1]. It is worth to point out that this represents a strong simplification w.r.t. the eigenvector representation in NABA approach, where the holding different representations [158] are equivalent to an "explicit representation" which is written in the form of a sum over partitions. This type of results in the SoV framework is even more important in the case of the Hubbard model. Indeed, there, algebraic approaches like NABA are mainly limited to the two particle case [123] and other exact results are accessible only via coordinate Bethe Ansatz 3 .
The paper is organized as follows. In section 2, we first shortly present the graded formalism for the superalgebras gl M|N and their fundamental representations, we sum up the main properties of the hierarchy of the fused transfer matrices and their reconstruction in terms of the fundamental 4 one. The SoV basis is then constructed in subsection 2.4 by using the integrable structure of these models. In subsection 2.5, we make some general statement about the closure and admissibility conditions to fix the transfer matrix spectrum for these quantum integrable models. In section 3, we specialize the discussion on the gl 1|2 model. We state our conjecture on the corresponding closure conditions and we present some first arguments in favor of it in subsection 3.1. Then, we treat in detail a special twisted case in subsection 3.2, for which we prove that the entire spectrum of the transfer matrix is characterized by our conjecture. Then, we give a reformulation of the spectrum in terms of the solutions to a quantum spectral curve equation. Moreover, for these representations, we show the completeness of the Nested Algebraic Bethe Ansatz, by proving that any eigenvalue can be rewritten in a NABA form using our Q-functions. In section 4, we derive an SoV basis for the Hubbard model with general integrable twist matrix having simple spectrum. Finally, in appendix A for the gl 1|2 model with general integrable twist matrices, we verify that the NABA form of eigenvalues satisfies the closure and admissibility conditions, implying its compatibility with our conjecture in the SoV framework. In appendix B, we present the proof that our conjecture indeed completely characterizes the transfer matrix spectrum for any integrable twist matrix having simple spectrum for the model defined on two sites while we verify this property by numerical computations for three sites. In appendix C, we recall how to derive the closure relation for the gl M|N case.

Separation of variables for integrable gl M|N fundamental models
Graded structures and Lie superalgebras are treated in great details in [214,215]. The quantum inverse scattering construction for graded models was introduced in [2,3], and summarized in many articles, see e.g. [9,216,217]. Details on Yangians structures for Lie superalgebras can be found in [218,219].

Graded formalism and integrable gl M|N fundamental models
(2.1) Vectors of V 0 are even, while vectors of V 1 are odd. Object that have a well defined parity, either even or odd, are called homogeneous. The parity map, defined on homogeneous objects, writes Maps between Z 2 -graded objects are called even if they preserve the parity of objects, or odd if they flip it. An associative superalgebra is a super vector space with an even multiplication map that is associative and the algebra has a unit element for the multiplication. For a superalgebra V we have V i V j ⊆ V i+j(mod 2) . A Lie superalgebra is a super vector space g = g 0 ⊕ g 1 equipped with an even linear map [ , ] : g ⊗ g → g that is graded antisymetric and satisfies the graded Jacobi identity.
The set of linear maps from V to itself is noted End V , and it is a Z 2 -graded vector space as well. It is an associative superalgebra with multiplication given by the composition. It is also a Lie superalgebra with the Lie super-bracket defined as the graded commutator between homogeneous objects for the multiplication which extends linearly to the whole space. As a Lie superalgebra, it is denoted gl(V ).
Tensor products The tensor product of two super vector spaces V and W is the tensor product of the underlying vectors spaces, with the Z 2 -grading structure given by This also defines the tensor product of associative superalgebras, being defined on the underlying vector space structure, but then we have to define an associative multiplication compatible with the grading. For A, B two associative superalgebras, the multiplication rule on A ⊗ B is given by for a 1 , a 2 ∈ A and b 1 , b 2 ∈ B homogeneous, and extends linearly to A ⊗ B. This rule of sign also appears in the action of A ⊗ B on V ⊗ W , where V is an A-module and W is a B-module. We have (2.6) for a ∈ A, b ∈ B, v ∈ V and w ∈ W . This rule extends naturally to N -fold tensor product, N > 2. For example, Lie superalgebra gl M|N Let V = C M|N be the complex vector superspace with even part of dimension M and odd part of dimension N . The general linear Lie algebra gl M|N = gl(C M|N ) is the Z 2 -graded vector space End C M|N with the Lie super-bracket defined by the graded commutator (2.3). We fix a homogeneous basis where v i is even for i ≤ M and odd for i ≥ M + 1 and we assign a parity to the index themselves for convenience:ī = 0 for i ≤ M andī = 1 for i ≥ M + 1.
The elementary operators e j i of gl M|N have parity p(e j i ) =ī +j(mod 2). They are defined by their action on the basis of V by They satisfy the following (2.10) Elements of gl M|N decompose on the elementary operators as where in the last term the sum over repeated indexes is omitted, as we will do in the following. Elements of gl M|N ⊗N writes Note that due to the sign rule (2.5), the coordinates A i 1 ...i N j 1 ...j N do not coincide with the components of the image of v j 1 ⊗ . . . ⊗ v j N in the tensored basis of V ⊗N : In the non-graded case, these two tensors would be identical. One may use the coordinates expression to check the parity of a given operator of gl M|N ⊗N . An (2.14) The supertrace is defined on the elementary operators as str e j i = (−1)jδ j i . Elements of gl M|N write as a block matrix Dual space Let us denote | i ≡ v i . The dual basis { j |} j=1,...,M+N is defined by The covectors are graded by p( j |) =j. For V ⊗N , the dual basis covectors have an additional sign in their definition such that it compensate for the permutation of vectors and covectors: Permutation operator The permutation operator has to take account of the grading when flipping vectors : Thus we have Remark the additional signs in the action as discussed (2.13). On an N-fold tensor product V 1 ⊗. . .⊗V N , with V i C M|N , the permutation operator P ab between spaces V a and V b writes where the number of identity operator is obvious by the context. The permutation operator is globally even. It satisfies the usual identities P 12 = P 21 , (2.24) P 12 P 13 = P 13 P 23 = P 23 P 12 , (2.25) P 13 P 24 = P 24 P 13 , P 12 P 34 = P 34 P 12 (2.26) which extend naturally to a N -fold tensor product.
The Y(gl M|N ) fundamental model The R matrix for the fundamental model of the Yangian Y(gl M|N ) writes It is of difference type and decomposes on elementary operators as ⊗I . . . ⊗ I, (2.29) using the notation (2.28). The R matrix is globally even and satisfies the Yang-Baxter equation Sometimes the Yang-Baxter equation is written in coordinates and is explicitely referred to as a "graded" version, such as in [3]: If K is an homogeneous matrix of gl M|N , we have This is a scalar version of the Yang-Baxter equation (2.30), where we put a trivial representation on the third space. For a spin chain of length N, we denote the Hilbert space by H = V 1 ⊗ . . . ⊗ V N and the auxiliary space by V 0 , all the V j superspaces being isomorphic to C M|N . Taking an even twist where ξ 1 , . . . , ξ N are the inhomogeneities of the chain. The monodromy is globally even as a product and tensor product of even operators. In coordinates, using the notation R ik jl of (2.28), it writes where all the signs from the multiplication of the operators actually vanish because of the eveness of R. This expression allow to show that the monodromy elements M i j (λ) ∈ End(H), as a matrix in the auxiliary space, are homogeneous of parity p(M i j (λ)) =ī +j. The Yang-Baxter scheme generalizes to the monodromy thanks to the global eveness of the R matrix, and we have The transfer matrix is obtained by taking the supertrace over the axuiliary space V 0 It is an even operator of End(H) as a sum of diagonal elements of the monodromy. Because the supertrace vanishes on the supercommutator, ones proves the commutation of the transfer matrices

The tower of fused transfer matrices
Tensor products of fundamental representations of gl M|N decompose in direct sum of irreducible subrepresentations (irreps). Young diagrams are used to carry out calculations with a mechanic proper to superalgebra, though very similar to the non graded case [183,[220][221][222]. Finite dimensional irreducible representations are labelled in a unique way by Kac-Dynkin labels, but the correspondence between Kac-Dynkin labels and Young diagrams is not one-to-one [223]. Admissible Young diagrams lie inside a fat hook domain pictured in figure 1 . Young diagrams can expand infinitely in both a and b directions, but the box (a ≥ N + 1, b ≥ M + 1) is forbidden, leading to the hook shape. Remark 2.1. For N = 0, the fat hook degenerates to a vertical strip, forcing a ≤ M. We recover the usual Young diagram indexation of gl(M) irreducible representations, though the diagrams are here displayed vertically, corresponding to the transpose of the usual gl(M) ones. This orientation is consistent with the one used in [183].
The tensoring procedure is called fusion in the context of the quantum inverse scattering method [169,170]. It is used to generate higher dimensional gl M|N -invariant R-matrices. The Yang-Baxter scheme is preserved, as degeneracy points of the fundamental R-matrix allow to construct the projectors P λ : (C M|N ) ⊗n → V λ , that extract the wanted subrepresentation λ, as a product of R-matrices. Fusing on the auxiliary space from M (K) 0 , we obtain new monodromy operators and thus new transfer matrices of End(H).
For a rectangular Young tableau corresponding to the point (a, b) ∈ H M|N , with b rows and a columns, the monodromy matrix reads is obtained by taking the supertrace over the a×b auxiliary spaces V The shifts in (2.41) are given by filling the fat hook as in figure 2. We then read the rectangular Young diagram column by column, top to bottom from left to right, and tensor the shifted monodromy (2.35) corresponding to the current box to the left of the previous ones. Figure 2 The domain H M|N is filled with multiples of the deformation parameter η. Starting from 0 in box (1, 1), we add η when moving down and −η when going right.
As said earlier, the projectors P λ can be constructed as a product of R-matrices, like in the non graded case [169,183,187]. For rectangular diagrams (a, b), we have 5 where i, j run on the boxes of the diagram of figure 2 column by column, top to bottom from left to right, and s i , s j are the shift contained in the boxes. All these transfer matrices commutes with each other, as consequence of the Yang-Baxter equation (2.30) being true for any irreps taken in the spaces 1, 2, and 3 Let us comment that through the very nice coderivative formalism [224,225], in a different but equivalent manner, these fused transfer matrices and the following fusion properties can be also derived.
In particular, we make use of this formalism in appendix C to verify (2.50). Among many others, the fused transfer matrices satisfy the following important properties [179][180][181][182][183]: Polynomial structure The generic fused transfer matrix T (λ) is polynomial in λ of degree abN, with (ab − 1)N central zeros given by and therefore factorizes as T Fusion equations There is a bilinear relation between transfer matrices associated to adjacent rectangular diagrams where the following boundary conditions are imposed for consistency These relations come from the Jacobi identity applied on the determinant form of the transfer matrices given by the Bazhanov-Reshetikhin formula (2.53), (2.54). All the fused transfer matrices outside the Inner-boundary condition As the correspondence between Young diagrams and irreps is not bijective, there exist non-trivial relations linking transfer matrices coming from distinct Young diagrams. This is especially the case for rectangular diagrams saturating one of the branch of the fat hookH M|N [183,184,223]. We shall call the first of these relation the inner-boundary condition, which writes and coincides with the central element called the quantum Berezinian, for N = M, and it plays a role similar to the quantum determinant in the non graded case [217,218]. As anticipated, we verify this relation in appendix C.

Reconstruction of fused transfer matrix in terms of the fundamental one
Here we want to recall that all these fused transfer matrices are completely determined in terms of the transfer matrix T (λ) by: then our statement is proven once we prove it for the transfer matrix of type T (λ). Let us use in the following the simpler notations and similarly The fused transfer matrices T (a) (λ) are of degree N. We have the following properties for these matrices: Lemma 2.1. The following asymptotics holds: where, in agreement with (2.43), the projectors admit the following iterative representations in terms of the R-matrix:

60)
In the inhomogeneities the following fusion relations holds: for any positive integer n.
Proof. The asymptotics are an easy corollary of the definition of the fused transfer matrices of type T Let us now prove the fusion relations in the inhomogeneities, for n = 1 the identity: are obtained by the fusion equations (2.47) just remarking that being: Then we can proceed by induction to prove the identity, let us assume that it holds for n ≥ 1 and let us prove it for n + 1, the relevant fusion identities reads: which being: read: which lead to our identities for n + 1 once we use the induction hypothesis for n The known central zeros and asymptotic imply that the interpolation formula in the N special points defined by the fusion equations to write T and (λ) by the fusion equations and the following interpolation formulae:

SoV covector basis for gl M|N Yang-Baxter superalgebra
In the next section we construct a separation of variables basis for the integrable quantum model associated to the fundamental representations of the gl M|N -graded Yang-Baxter algebra. Here, the construction follows the general ideas presented in the Proposition 2.4 of [1] as the main steps in the proof rely on the reduction of the fundamental R-matrix to the permutation operator and on the centrality of the R-matrix asymptotics which hold true also in the current gl M|N -graded case.
Let K be a (M+N )×(M+N ) square matrix solution of the gl M|N -graded Yang-Baxter equation of the block form (2.34), then we use the following notation where K J is the Jordan form of the matrix K Then, the following similarity relation holds for the fundamental transfer matrices: i.e. they are isospectral. We can now state our main result on the form of the SoV basis: then for almost any choice of the covector S| and of the inhomogeneities under the condition (2.77), the following set of covectors: forms a covector basis of H. In particular, we can take the state S| of the following form, according to (2.18): where:

88)
as soon as we take: (λ) transfer matrix spectrum is simple.
Proof. As in the non graded case, the identity: holds true as a direct consequence of the definition of the transfer matrix T (K) 1 (λ) and the properties From this point we can essentially follow the proof of Proposition 2.4 of [1]. Indeed, the condition that the set (2.86) form a covector basis of H is equivalent to where we have defined: (λ) a polynomial in the inhomogeneities {ξ i } and in the parameters of the twist matrix K, the same statement holds true for the determinant on the l.h.s. of (2.92), which is moreover a polynomial in the coefficients S|e j of the covector S|.
Then it follows that the condition (2.92) holds true for almost any value of these parameters if one can show it under some special limit. Using this argument, the form of the transfer matrix in the inhomogeneities (2.90) and the central asymptotics of the gl M|N -graded R-matrix one can show that a sufficient criterion is that the following determinant is non-zero: (2.95) Because K is even and thanks to the expression (2.87) of S |, we have which is clearly nonzero under the condition that the twist K has simple spectrum and when (2.89) holds. The simplicity of the transfer matrix spectrum is then a trivial consequence of the fact that the set of covectors (2.86) is proven to be a basis. Indeed, it implies that given a generic eigenvalue t(λ) of T where |t and t| are the unique eigenvector and eigencovector associated to it.
Proof. Following the proof of Proposition 2.5 of [1] the non-ortogonality condition (2.100) can be derived. Such condition together with the simplicity of the spectrum implies that we cannot have non-trivial Jordan blocks in the transfer matrix spectrum so that it must be diagonalizable and with simple spectrum.

On closure relations and SoV spectrum characterization
In the previous two subsections, we have shown how the transfer matrix T (λ) associated to general inhomogeneous representations of the gl M|N -graded Yang-Baxter algebra allows to reconstruct all the fused transfer matrices (mainly by using the known fusion relations (2.61) and (2.62)). Moreover, we have shown that T (K) 1 (λ) allows to characterize an SoV basis, which also implies its spectrum simplicity or diagonalizability and spectrum simplicity if the twist matrix K is, respectively, with simple spectrum or diagonalizable with simple spectrum.
This analysis shows that the full integrable structure of the gl M|N -graded Yang-Baxter algebra can be recasted in its fundamental transfer matrix as well as the construction of quantum separation of variables. There is however still one missing information, which is a functional equation, or a discrete system of equations, allowing the complete characterization of the transfer matrix spectrum. As already mentioned, the fusion relations (2.47) alone only give the characterization of higher transfer matrices in terms of the first one. Some further algebra and representation dependent rules are required in order to complete them and extract a closure relation on the transfer matrix.
In the case of the quantum integrable models associated to the fundamental representations of the gl M and U q ( gl M ) Yang-Baxter and reflection algebras, such a closure relation comes from the quantum determinant [1,44,45,[149][150][151]. Indeed, P − 1...M is a rank 1 projector in these cases, implying that the corresponding transfer matrix T (K) (M) (λ) becomes a computable central element of the Yang-Baxter algebra, namely the quantum determinant. Then, substituting the quantum determinant in the fusion equation (2.61) for n = M−1 and using the same interpolation formulae for the higher fused transfer matrix eigenvalues, we produce a discrete system of polynomial equations with N equations in N unknowns which was proven [1,44,45,[149][150][151] to completely characterize the transfer matrix spectrum in quantum separation of variables. In the case of non-fundamental representations 6 of the same algebras the closure relation comes instead with the appearing of the first central zeros in the fused transfer matrices of type T (K) n (λ). In [152], this analysis has been developed in detail in the case of M = 2. There, it has been shown that imposing the central zeros of the fused transfer matrix T (K) 2s+1 (λ), for a spin s ≥ 1 representation, a discrete system of polynomial equations with N equations in N unknowns is derived for the transfer matrix eigenvalues. The set of its solutions completely characterizes the transfer matrix spectrum in quantum separation of variables. In the nonfundamental and cyclic representations of the U q ( gl M ) Yang-Baxter algebra for q a root of unit such closure relation comes from the so-called truncation identities. For M = 2, it has been shown in [39] how these identities emerge and are proven in the framework of the quantum separation of variables and how they are used to completely characterize the transfer matrix spectrum. While in [43] and [58,59] these results have been extended, respectively, to the most general cyclic representations of the U q ( gl 2 ) Yang-Baxter algebra and reflection algebra.
In the case of integrable quantum lattice models associated to the fundamental representations of the gl M|N -graded Yang-Baxter algebra, the natural candidate for the closure relation is the innerboundary condition (2.50). Indeed, once we impose it on the eigenvalues t (λ). Moreover, the inner-boundary condition (2.50) involved the quantum Berezinian as a central element hence playing a role similar to the quantum determinant in the bosonic case. More precisely, we can introduce the following polynomials: and from them recursively the following higher polynomials and N +m (λ|{x a }) = 0, ∀λ ∈ C and n, m ≥ 1. (2.107) In the next section we conjecture that the above system of functional equations completely characterizes the transfer matrix spectrum in the case of the gl 1|2 -graded Yang-Baxter algebra. We prove this characterization for some special class of twist matrices while we only give some first motivations of it for general representations and in appendix B we verify it for quantum chains with two and three sites. Let us also mention that this conjecture can be checked explicitly for the simple gl 1|1 case. It would be interesting in this respect to elucidate the relation of our method with the one developed recently in [188].

General statements and conjectured closure relation for general integrable twist
We use this subsection to clarify and justify the following Conjecture Taken the general twisted gl 1|2 -graded Yang-Baxter algebra, then t 1 (λ|{x a }) is an eigenvalue 8 of the transfer matrix T the inner-boundary condition (2.50) reads: where: Moreover, we have the following expressions 9 for the asymptotics of the fused transfer matrices Then imposing it for the corresponding eigenvalues we get which, once we express the eigenvalues t 2 (λ + η) = t 2 (λ)t 2 (λ + η) − t 1 (λ + η)t 3 (λ), (3.5) takes the following closure relation form: Now, using the interpolation formulae (2.101) and (2.102) for the transfer matrix eigenvalues, we get that the closure relation is indeed a functional equation whose unknowns coincide with the x i≤N ≡ t in the gl 1|2 -graded case under consideration. According to our Conjecture the transfer matrix spectrum coincides with the set of solutions to the functional equations (3.6) and (3.7) in the unknowns t (K) 1 (ξ i≤N ). This spectrum characterization for general twist matrix will be proven by direct action of the transfer matrices on our SoV basis in our next publication. Here we present some arguments in favour of it while in the next subsections 3.2 we prove it for a special choice of the twist matrix.
Let us consider the case of a twist matrix K invertible with simple spectrum 10 , then our SoV basis can be written as follows: (ξ n ) are invertible and so we can write our original vector S| as it follows: (3.9) If the t 1 (ξ i≤N ) solve the closure relation (3.6) and the null out-boundary conditions (3.7), to prove that a vector |t characterized by is indeed a transfer matrix eigenvector the main point is to be able to reduce the covectors containing a fourth order power of T (K) 1 (ξ n ) in those of the SoV basis, of maximal order three. Moreover, this reduction must came from relations which are satisfied identically by both the fused transfer matrices and the functions t r (λ|{x a }) defined in (2.102), in terms of the given solution t 1 (ξ i≤N ). Indeed, this is exactly what it is done by the closure relation for the transfer matrix: and by the corresponding one (3.6) for the eigenvalues. As one can easily remark that (3.11) and (3.6) are both of fourth order, respectively, in the T (K) 1 (ξ n ) and t 1 (ξ n ) on the right hand side while they are both of third order on the left hand side.
In appendix A, we will verify that Nested Algebraic and Analytic Bethe Ansatz are indeed compatible with these requirements, i.e. the functional ansatz for the eigenvalues t 1 (λ) indeed satisfies the closure relation (3.6) and the null out-boundary conditions (3.7). There, we moreover argue the completeness of the Bethe Ansatz which is compatible with our Conjecture.
It is also worth to mention that we have verified our Conjecture on small lattices, up to three sites. More in detail, we have solved the discrete system of N equations in the N unknowns t

SoV spectrum characterization for non-invertible and simple spectrum twist matrix
Let us study here the spectral problem for the transfer matrices associated to the fundamental representations of the gl 1|2 -graded Yang-Baxter algebra in the following class of non-invertible but having simple spectrumK twist matrices:K with K 2 any invertible, diagonalizable and simple 2 × 2 matrix, i.e. it holds: (λ) coincides with the following set of polynomials where S T (K) is the set of solutions to the following system of N cubic equations: (λ) is diagonalizable and with simple spectrum. For any t 1 (λ) ∈ Σ T (K) , the associated and unique eigenvector |t (up-to normalization) has the following wavefunctions in our SoV covector basis: (3.16) Proof. The main identity to be pointed out here is the following one: (λ) ≡ 0, (3.17) due to the closure relation (3.11) being k 1 = 0. So that the fusion equations (2.62) for n = 1 and n = 2 read: Now it is easy to verify that the system of equations (3.15) just coincides with the above fusion conditions once imposed to functions which have the analytic properties (polynomial form and asymptotics) of eigenvalues. So that it is clear that any eigenvalue has to satisfy them and one is left with the proof of the reverse statement. This proof can be done just showing that the state |t of the form (3.16) is indeed an eigenvector of the transfer matrix, i.e. that it holds: Finally, let us point out that Theorem 2.1 implies also that the transfer matrix T  which coincides with the closure relation (3.6) for the k 1 = 0 case under consideration. Then it is easy to verify that the null out-boundary conditions (3.7) are also verified. Note for example that the condition (3.21) together with the interpolation formula (2.102) for k 1 = 0 implies: so that the interpolation formula (2.104) implies: i.e. the null out-boundary conditions (3.7) for n = m = 0 and similarly for the others. system of polynomial equations admits 3 N solutions if the N polynomials, defining the system, have no common components 12 . So, under the condition of no common components, there are indeed exactly 3 N distinct solutions to the above system and each one is uniquely associated to a transfer matrix eigenvalue. The proof of the condition of no common components can be done following exactly the same steps presented in appendix B of [150]. (λ) is identically zero in these representations associated to noninvertible simple spectrum twist matrixK means, in particular, that it is central so that the algebra shows some strong resemblance to the twisted representations of the gl 3 Yang-Baxter algebra. In fact, taking the gl 3 -representation associated to the twist matrix K = −K and η = −η, then the SoV characterization of the spectrum implies that the transfer matrix T (λ|η), associated to the gl 1|2 -representation, is isospectral to T (−K) (2) (λ| − η), associated to the gl 3 -representation.
It is worth remarking that the same type of duality indeed holds between the gl 1|N -graded and the gl N +1 non-graded Yang-Baxter algebra when associated to the non-invertible but simple (N + 1) × (N + 1) twist matrixK with first eigenvalue zero. More in detail, we have the isospectrality of the transfer matrices T (m) (λ| − η), associated to the gl N +1 -representation, for any 1 ≤ m ≤ N . This in particular implies that we can characterize completely as well the spectrum of the transfer matrices of the gl 1|N -graded Yang-Baxter algebra for this special class of twist matrices just using the results of [149]. Then the results of the next two subsections can be as well generalized to these special classes of gl 1|N -graded Yang-Baxter algebras.

The quantum spectral curve equation
The transfer matrix spectrum in our SoV basis is equivalent to the quantum spectral curve functional reformulation as stated in the next theorem. such that t 1 (λ), t 2 (λ|{t 1 (ξ a≤N )}) and ϕ t (λ) are solutions of the following quantum spectral curve functional equation: where we have defined andᾱ is a nonzero solution of the characteristic equation: i.e.ᾱ = −k 2 orᾱ = −k 3 is a nonzero eigenvalue of the twist matrix −K. Moreover, up to a normalization, the common transfer matrix eigenvector |t admits the following separate representation: Proof. From the above Remark 3.3, we have that this theorem is a direct consequence of the Theorem 5.2 of [1]. To make the comparison easier, one has just to take the quantum spectral curve characterization of Theorem 4.1 of [149] and use it in the case n = 3, K → −K, η → −η to get the quantum spectral curve associated to the non-invertible simple spectrum twistK. Here, t 1 (λ) is the eigenvalue associated to the gl 1|2 -transfer matrix T (λ| − η). Then, just removing the common zeros, in the three nonzero terms of the equation and making the common shift λ → λ − 2η, we get our quantum spectral curve equation (3.26). Let us recall that the main elements in the proof of the theorem rely on the fact that the quantum spectral curve is equivalent to the following 3N conditions: once the asymptotics are fixed as stated above in this theorem. It is then easy to observe that the compatibility of this system of equations is equivalent to impose that t 1 (λ) and t 2 (λ|{t 1 (ξ a≤N )}) satisfy the fusion equations:

Completeness of Bethe Ansatz solutions by SoV
As detailed in the introduction, Nested Algebraic and Analytic Bethe Ansatz have been used to study the spectrum of the model associated to the fundamental representation of the gl M|N Yang-Baxter superalgebra. For the fundamental representation of the gl 1|2 Yang-Baxter superalgebra associated to a simple and diagonalizable twist matrix K, let us recall here the form of the Bethe Ansatz equations [3,[132][133][134]: where (3.38) and the Bethe Ansatz form of the transfer matrix eigenvalue defined by It is worth to observe that being it follows that the function t 1 (λ|{λ j≤L }, {µ h≤M }) has only apparent simple poles in the λ j≤L and µ h≤M . The regularity of t 1 (λ|{λ j≤L }, {µ h≤M }) for λ = λ j≤L is implied by the Bethe equation (3.36) while the regularity of t 1 (λ|{λ j≤L }, {µ h≤M }) for λ = µ j≤M is implied by the Bethe equation (3.37). Hence 13 t 1 (λ|{λ j≤L }, {µ h≤M }) is a polynomial of degree N with the correct asymptotic for a transfer matrix eigenvalue, i.e. it holds: So that the above ansatz is indeed consistent with the analytic properties enjoyed by the transfer matrix eigenvalues. Now, we show that the Bethe ansatz solutions are complete as a corollary of the completeness of the derived quantum spectral curve in the SoV framework, for the class of representations considered in this section. More precisely it holds the next Moreover, for any t 1 (λ) ∈ Σ T (K) the associated Bethe Ansatz solution {{λ j≤L }, {µ h≤M }} is unique and satisfies the admissibility conditions: Proof. This corollary directly follows from our previous theorem. The proof is done pointing out the consequences of the special form of the fusion equations for these representations associated to these non-invertible twists. In particular, from the fusion equations (3.33) and (3.34), which have to be satisfied by all the transfer matrix eigenvalues, we derive the following equation on the second transfer matrix eigenvalues only: Being by definition t (K) 2 (λ) a degree 2N polynomial in λ, zero in the points ξ a −η for any a ∈ {1, ..., N}, it follows that a solution to (3.49) can be obtained iff for any a ∈ {1, ..., N} there exists a unique h a ∈ {−1, 0} such that: (3.50) So we have that the system (3.49) has exactly 2 N distinct solutions associated to the 2 N distinct 2,h (λ). So, for any fixed h ∈ {0, −1} N , we have that the eigenvalues of the transfer matrix have the following form: 1,h (λ) is a degree N − m h polynomial in λ, and the function ϕ t,h (λ) associated by the quantum spectral curve to the eigenvalue t (K) 1,h (λ) has the form: Then, simplifying common prefactors, the quantum spectral curve rewrite as it follows: where we have defined:d (λ − ξ π h (a) ). (3.56) So, once we choseᾱ = −k 2 , we get the following representation of the transfer matrix eigenvalue: It is now trivial to verify that this coincides with the Bethe ansatz form (3.45) for k 1 = 0, once we fix: Clearly by definition Q
We can define the following monodromy matrix: where H = N n=1 V An , V An ∼ = C 4 . Then the transfer matrix: defines a one parameter family of commuting operators.

Our SoV covector basis
The general Proposition 2.4 and 2.5 of [1] for the construction of the SoV covector basis and the diagonalizability and simplicity of the transfer matrix spectrum can be adapted to the inhomogeneous Hubbard model. Let us denote with K J (a, α, β, γ) the diagonal form of the matrix K(a, α, β, γ) and W K the invertible matrix defining the change of basis to it: clearly W K is the identity for a = 1, then the following theorem holds: forms a covector basis of H, for almost any choice of S|. In particular, we can take the state S| of the following tensor product form: simply asking x y z w = 0.
Proof. We have just to remark that also in this case the following identity holds: T (K) (ξ n ) = R An,A n−1 (ξ n |ξ n−1 ) · · · R An,A 1 (ξ n |ξ 1 )K An R An,A N (ξ n |ξ N ) · · · R An,A n+1 (ξ n |ξ n+1 ). (4.21) Let us now point out that e h(λ,η) is an algebraic function of order two in η and e λ . Then the determinant of the matrix whose rows are the elements of these covectors in the elementary basis is also an algebraic function of {e ξm } m∈{1,...,N} and η. So that showing that this determinant is nonzero for a specific value of η one can prove that it is nonzero for almost any value of η and of the others parameters, i.e. the inhomogeneities satisfying (2.77) and the parameters α, β, γ of the twist matrix, for a = 1, 2, 3. We can study for example the case η = 0. In this case h(λ, η) has the following two different determinations: h(λ, η = 0) = 0, iπ/2 modiπ. (4.22) Note that in both the cases, we have that it holds: tanh(h(λ, 0) + h(µ, 0)) = 0 (4.23) so that the Shastry's R-matrix reduces to the tensor product of two XX R-matrix, i.e. it holds: In turn this implies that: so that we can repeat the same type of proof of the general Proposition 2.4 of [1] to show that for a covector S|, of the above tensor product form, the determinant of the full matrix factorizes in the product of the determinants of 4 × 4 matrices which are nonzero due to the simplicity of the spectrum of the matrix K. This already implies the w-simplicity of the transfer matrix T (K) (λ) then in the case η = 0 we can prove the non-ortogonality condition: for any transfer matrix eigenvector by the same argument developed in general Proposition 2.5 of [1], which implies the diagonalizability and simplicity of the transfer matrix spectrum for η = 0 and so for almost any value of η and of the others parameters of the representation.
Let us briefly comment about the consequences of the existence of such a basis. The first important point to stress is that whenever we have an eigenvalue for the transfer matrix, we can write the corresponding eigenvector in the above basis. It means that if we compute a set of solutions to the Nested Bethe Ansatz equations we can immediately write the transfer matrix eigenvalue and hence the corresponding eigenvector; in particular it will be a true eigenvector as soon as it is non zero. This could be of great use in practice when dealing with finite chains with a number of sites greater than the values accessible by direct diagonalization. In particular, scalar products and form factors could become accessible, at least numerically, from this procedure. For using the above basis on a more fundamental, analytical level, one needs to obtain the complete set of fusion relations that lead to the full closure relations enabling to compute the action of the transfer matrix in the SoV basis (see the discussion on this point given in [1,149]). This should lead to the full characterization of the spectrum. These fusion relations being rather involved for the Hubbard model, due in particular to the intricate dependence on spectral parameters [131], we will come back to this question in a future publication. Let us nevertheless anticipate that the results obtained for the gl M|N case will be of direct importance when dealing with the Hubbard model.

A Compatibility of SoV and Bethe Ansatz framework
In this appendix, we verify how the results obtained in the Nested Algebraic and Analytic Bethe Ansatz framework for the gl 1|2 Yang-Baxter superalgebra are compatible with the conjectured spectrum characterization in the SoV basis. This analysis is done in the fundamental representations of the gl 1|2 Yang-Baxter superalgebra associated to generic diagonalizable and simple spectrum twist matrices.

A.1 Compatibility conditions for higher transfer matrix eigenvalues
Here we use the Bethe Ansatz form (3.39) of the transfer matrix eigenvalues together with the Bethe ansatz equations (3.36) and (3.37) to describe the eigenvalues of the higher transfer matrices in order to verify that they satisfy both the null out-boundary (3.7) and the inner-boundary (3.4) conditions. Under these hypothesis, we get the following lemma: Lemma A.1. The eigenvalues of the higher transfer matrices admits the following representation in terms of the Λ i (λ) functions: Proof. We have already observed that due to the Corollary 2.1 the eigenvalues of the higher transfer matrices admits the interpolation formulae (2.102) in terms of the transfer matrix eigenvalue t 1 (λ). Equivalently, given t 1 (λ) an eigenvalue of the transfer matrix then those of the higher transfer matrices are of the form wheret n (λ) are degree N polynomials in λ, fixed uniquely by the recursive equations: t n+1 (ξ a ) = t 1 (ξ a )t n (ξ a + η), (A. 6) and the known asymptotics: So to prove the above Bethe Ansatz form for the higher transfer matrix eigenvalues we have just to verify these conditions. Concerning the asymptotic behavior, from the r.h.s. of formula (A.1) and (A.3) we get: so they are satisfied. So we are left with the proof of the fusion properties. It is easy to remark that by the definition of the Λ i (λ) it follows that the t n (λ) indeed factorizes the coefficients n−1 r=1 d(λ + rη), for n ≥ 2. Let us now show that is indeed a degree N polynomials in λ. This is the case iff the residues of this expression in the zeroes of Q 1 (λ) are vanishing, namely iff the following identities hold: and this is the case thanks to the Bethe equation (3.36) in λ j being: .
Similarly, we have that which is a degree N polynomials in λ due to the identity (3.36). So to show that the t n (λ) satisfy the characterization of the higher eigenvalues, we have just to verify that their values in the inhomogeneities agrees with (A.5) and (A.6). Indeed, it holds: where in the last line we have used the Bethe Ansatz form of t 1 (λ) and similarly: Here we have explicitly rewritten the eigenvalues form in Bethe Ansatz approach for the higher transfer matrix T (K) n (λ), by using the fusion we can easily derive those of the others. Now we are interested in showing that these expressions for the higher eigenvalues indeed imply the null outboundary (3.7) and the inner-boundary (3.4) conditions. Indeed, we have the following lemma:  Proof. By using the result of the previous lemma it is easy to show the following null conditions are satisfied: t (2) 3+n (λ) = 0, ∀n ≥ 0, (A. 19) indeed, the condition (A.3) implies: so that, in particular, it holds: t 3+n (λ) = (t 2+n (λ)/t 1+n (λ + η))(t 2+n (λ + η)), ∀n ≥ 1, (A. 21) or equivalently: which by the fusion equations implies the above null conditions. Similarly, we can derive all the others null out-boundary conditions (3.7).
Let us now show the inner-boundary condition, from the formula (A.1) we can write: and so which is equivalent to our closure relation: (k 1 t 1 (λ + η) + k 3 k 2 d(λ))t 3 (λ) = k 1 t 2 (λ)t 2 (λ + η), (A. 25) and taking into account the fusion equation: 3 (λ), (A. 26) we are led to the required identity: It should be noted that all these relations can be proven in a pure algebraic way using the general constructions of T and Q operators and the various relations they satisfy, as given in [225,[227][228][229]. Hence the computations presented here, although quite instructive could be considered merely as consistency checks.

A.2 On the relation between SoV and Nested Algebraic Bethe Ansatz
Let us consider the fundamental representation of the gl 1|2 Yang-Baxter superalgebra associated to generic values of the inhomogeneities {ξ a≤N }, satisfying the condition (2.77), and of the eigenvalues k 1 , k 2 and k 3 of a simple and diagonalizable twist matrix K. Note that in this case the transfer matrix is similar to the transfer matrix associated to a diagonal twist with entries the eigenvalues k 1 , k 2 and k 3 of K to which Nested Algebraic Bethe Ansatz (NABA) directly applies. Therefore, the following discussion on the connection between the SoV description and the NABA can be directly addressed in this diagonal case.
Let us recall that in the NABA framework, given a solution {{λ j≤L }, {µ h≤M }} of the Bethe Ansatz equations (3.36) and (3.37) satisfying the pair-wise distinct conditions (3.46), then the associated Bethe Ansatz vector |t {λ j≤L },{µ h≤M } is proven to satisfy the identity with t 1 (λ|{λ j≤L }, {µ h≤M }) defined in (3.39), so that it is a transfer matrix eigenvector as soon as it is proven to be nonzero. Then, such a Bethe Ansatz vector has in our SoV basis the following characterization:  {λ j≤L },{µ h≤M } is nonzero. Relying on some already existing results in the literature, we want to present a reasoning that allows to argue that the completeness of the Bethe Ansatz in the SoV framework, for the special representations of subsection 3.2, indeed implies the completeness for the NABA spectrum description. The reasoning goes as follows. In [1,149] we have shown in general that the SoV characterization of the transfer matrix eigenvectors allows for an Algebraic Bethe Ansatz rewriting on a well defined reference state, see for example section 5 of [149]. Adapting to the current fundamental gl 3 -representation associated the twist matrix K = −K and η = −η the analysis of [230], it can be argued 14 that the B-operator defined in our SoV basis indeed coincides with the one defined by Sklyanin [24]. Then by adapting the results presented in [168], one can deduce that these SoV eigenvectors rewritten in an Algebraic Bethe Ansatz form, in terms of the Sklyanin B-operator, in turn coincides (up to nonzero normalization) with Nested Algebraic Bethe Ansatz vectors, associated to the same Bethe Ansatz solutions. If implemented with all details this reasoning shows the completeness of the Nested Algebraic Bethe Ansatz as a consequence of the completeness of the SoV characterization derived in subsection 3.2.2 for the fundamental representations associated to non-invertible but simple spectrum K twist matrices, with eigenvalues satisfying (3.13).
It is also worth to comment that once the NABA completeness is derived for these special representations, it can be derived for general gl 1|2 -representations by adapting to them the proof given in [226] for the gl 2 fundamental representations associated to general diagonalizable twist. Indeed, one of the main ideas of the proof in [226] is that for a special value of the twist parameter, one can characterize the set of isolated Bethe Ansatz solutions that produce nonzero Bethe vectors, and which is proven to be complete. Then, the results on the completeness of Bethe Ansatz solutions by the SoV approach, derived in subsection 3.2.2, and the above argument on the NABA completeness for these gl 1|2 -representations associated to the twist matricesK can be as well the starting point for the proof of completeness by deformation w.r.t. the twist parameters like in [226]. Finally, let us add that relations with [188] would be interesting to explore.

B Verification of the Conjecture for the general twists up to 3 sites
Here, we make a verification of our conjecture on the form of the closure relations for the general twisted representation of the gl M|N Yang-Baxter superalgebra, in the case M = 1 and N = 2 for small chain representations, i.e. for a chain having up to N = 3 sites. The verification is done in the following way, we impose the closure relation (3.6) in N pairwise different values 15 of λ to the polynomials (2.101) and (2.102) for n = 1, 2. This determines a system of N polynomial equations of order 4 in the N unknown which are the values of the polynomial (2.101) in the inhomogeneities. We solve this system of equations by Mathematica and we select the solutions which generate polynomials (2.102) which satisfy the null out-boundary conditions 3.7. Our analysis shows that it is enough to impose 3.7 for n = m = 0 to select the correct solutions which generate exactly the N 3 different eigenvalues of the diagonalizable and simple spectrum transfer matrix T (K) (λ), obtained by diagonalizing it exactly with Mathematica. For N = 1, 2 the results of both the approaches are analytic and we present them here for the interesting N = 2 case. While for N = 3 we have verified our statements for different values of the parameters, i.e. the inhomogeneity parameters and the three eigenvalues of the twist matrix.

C Derivation of the inner-boundary condition
One may use the coderivative formalism introduced in [224] and developed in [225] to derive the innerboundary condition. The coderivative formalism allows to construct the transfer matrices associated to a given irreducible representation on the auxiliary space by acting on the associated character evaluated at the twist matrix. For a rectangular Young tableau (a, b), we have in our notation Let us take g = diag(x 1 , . . . , x M , y 1 , . . . , y N ) a diagonal twist. For k ≥ 1, the characters of the rectangular representations (a, b) which saturate an arm of the fat hook write [184]  (λ) matrices, we obtain (2.50).