An exactly solvable supersymmetric spin chain of BC_N type

We construct a new exactly solvable supersymmetric spin chain related to the BC_N extended root system, which includes as a particular case the BC_N version of the Polychronakos-Frahm spin chain. We also introduce a supersymmetric spin dynamical model of Calogero type which yields the new chain in the large coupling limit. This connection is exploited to derive two different closed-form expressions for the chain's partition function by means of Polychronakos's freezing trick. We establish a boson-fermion duality relation for the new chain's spectrum, which is in fact valid for a large class of (not necessarily integrable) spin chains of BC_N type. The exact expressions for the partition function are also used to study the chain's spectrum as a whole, showing that the level density is normally distributed even for a moderately large number of particles. We also determine a simple analytic approximation to the distribution of normalized spacings between consecutive levels which fits the numerical data with remarkable accuracy. Our results provide further evidence that spin chains of Haldane-Shastry type are exceptional integrable models, in the sense that their spacings distribution is not Poissonian, as posited by the Berry-Tabor conjecture for"generic'' quantum integrable systems.


Introduction
In the last few years exactly solvable (or integrable) supersymmetric spin chains and their associated dynamical models have been the subject of extensive research in connection with different topics of current interest, such as the theory of strongly correlated systems [1,2] or the AdS-CFT correspondence [3]. Among these models, the supersymmetric versions of the celebrated Haldane-Shastry chain [4,5] and its rational counterpart proposed by Polychronakos [6] and Frahm [7] occupy a distinguished position, due to the rich mathematical structures at the heart of their highly solvable character. The Haldane-Shastry (HS) chain was introduced in an attempt to construct a simple one-dimensional model whose ground state coincided with Gutzwiller's variational wave function for the Hubbard model in the limit of large on-site interaction [8,9,10]. It can also be obtained in this limit from the Hubbard model with long-range hopping studied in Ref. [11], in the half-filling regime. The original HS chain describes N spin 1/2 particles equally spaced on a circle, with pairwise interactions inversely proportional to the square of the chord distance. Its Hamiltonian can be written as where J i ≡ 1 2 (σ x i , σ y i , σ z i ), σ α i is a Pauli matrix at site i, and the summation indices run from 1 to N (as always hereafter, unless otherwise stated). A natural generalization of this chain to su(m) spin [12,13] is obtained by taking J i = (J 1 i , . . . , J m 2 −1 i ), where {J α i } are the generators of the fundamental representation of su(m) at site i. In fact, with the usual normalization tr(J α k J γ k ) = 1 2 δ αγ we have where S ij is the operator permuting the i-th and j-th spins, so that the su(m) spin Hamiltonian (1) is a linear combination of spin permutation operators.
The spectrum of the spin 1/2 HS Hamiltonian was numerically analyzed in the original papers of Haldane and Shastry. The more general su(m) spin case was studied in a subsequent publication by Haldane et al. [14], who empirically found a complete description of the spectrum and explained its high degeneracy by an underlying Y(sl m ) Yangian symmetry. These results were rigorously established in Ref. [15] by constructing a transfer matrix using the Dunkl operators [16,17] of the (trigonometric) Sutherland dynamical model [18,19]. Although this approach yields an explicit formula for the energies in terms of the so-called motifs, the computation of their corresponding degeneracies becomes quite cumbersome when m > 2. This is probably the reason why the partition function of this chain was computed only very recently [20], using a procedure known as Polychronakos's freezing trick which does not rely on the explicit knowledge of the spectrum [6,21].
The key idea behind the freezing trick is exploiting the connection between the HS chain and the Sutherland spin dynamical model in the strong coupling limit. Indeed, in this limit the particles in the latter model "freeze" at the coordinates of the (unique) equilibrium of the scalar part of the potential, which are exactly the chain sites. Thus, in this limit the dynamical and the spin degrees of freedom decouple, so that the spectrum of the spin Sutherland model becomes approximately the sum of the spectra of the scalar Sutherland model and the HS chain. This observation yields an explicit formula for the partition function of the HS spin chain as the strong coupling limit of the quotient of the partition functions of the spin and the scalar Sutherland models, both of which can be computed in this limit. In fact, this method was first applied by Polychronakos to evaluate the partition function of the su(m) spin chain associated with the spin Calogero (rational) model [22,23], usually known in the literature as the Polychronakos-Frahm (PF) chain, whose Hamiltonian (in the antiferromagnetic case) is given by The sites ζ i of this chain are the coordinates of the equilibrium of the scalar part of the Calogero potential, which coincide with the zeros of the Hermite polynomial of degree N [24].
Both the HS and the PF chains admit a natural su(m|n) supersymmetric version, in which there are m bosonic and n fermionic degrees of freedom at each site [25,26,27,28,29]. In practice, the Hamiltonians of these chains are obtained by replacing the spin permutation operators S ij by suitable supersymmetric generalizations thereof (cf. Eq. (4) below). The partition functions of the su(m|n) PF and HS chains have also been computed in closed form using the freezing trick method outlined above [25,27]. These partition functions were then used to establish a remarkable boson-fermion duality relation for the spectrum of each of these chains [25,28].
The original HS and PF chains (and their supersymmetric versions) are associated with the A N −1 root system, since the interaction between any two spins depends only on the difference of their site coordinates. It is also of interest to consider generalizations of these chains to the BC N extended root system, in which the spins interact not only among themselves but also with their mirror images with respect to a reflecting end located at the origin. In the non-supersymmetric case, the HS chain of BC N type was first discussed in Ref. [30], and its partition function was evaluated by means of the freezing trick in Ref. [31]. Likewise, the partition function of the su(m) PF chain of BC N type introduced in [32] was recently computed in closed form by the authors [33]. To the best of our knowledge, however, the supersymmetric counterparts of both the HS and PF chains of BC N type have not been studied so far. One of the main purposes of this paper is precisely to fill this gap in the case of the PF chain. More precisely, we shall apply the freezing trick to compute the partition function of this chain in closed form, relating it to the corresponding partition function of the su(m|n) PF chain of A N −1 type. We shall also extend the boson-fermion duality mentioned above to a large class of supersymmetric spin chains of BC N type which includes the PF chain as a particular case.
The spectra of the non-supersymmetric spin chains of Haldane-Shastry type whose partition function has been explicitly computed share some statistical properties that we shall now discuss. In the first place, the density of energy levels follows the Gaussian law with remarkable accuracy even for a moderately large number of particles [20,31,33,34]. Secondly (except for the HS chain of BC N type, which has not yet been studied in this respect), the density of spacings p(s) between consecutive (normalized) levels follows a characteristic distribution essentially different from Poisson's law p(s) = e −s . This fact is somewhat surprising, in view of the integrability of these chains [6,32] and the long-standing conjecture of Berry and Tabor [35], according to which the spacings distribution of a "generic" quantum integrable system should be Poissonian. The validity of this conjecture has been verified for a number of important integrable systems, including the Heisenberg chain, the t-J model, the Hubbard model [36], and the chiral Potts model [37]. On the other hand, the non-Poissonian character of the spacings distribution has recently been observed for the su(m|n) HS chain of A N −1 type [27], whose density of spacings is qualitatively similar to that of the non-supersymmetric HS and PF chains. Another important aim of this paper is the study of the level density and the spacings distribution for the su(m|n) PF chain of BC N type. Our main conclusion in this respect is that this chain behaves in essentially the same way as other chains of Haldane-Shastry type previously studied. In other words, when the number of spins is sufficiently large the level density follows with great accuracy the Gaussian distribution, and the spacings distribution also deviates substantially from Poisson's law. As a matter of fact, we shall see that the spacings density is well approximated by the same "square root of a logarithm" distribution as the original HS and PF chains, and the su(m) PF chain of BC N type [33,34]. This lends further support to the fact that spin chains of Haldane-Shastry type are exceptional integrable systems from the point of view of the Berry-Tabor conjecture.
The paper is organized as follows. In Section 2 we recall the definition of the supersymmetric spin permutation operators [25], and introduce new supersymmetric spin reversal operators. With the help of these operators, we construct the su(m|n) PF chain of BC N type and explain how the freezing trick can be used to compute its partition function from that of its associated spin dynamical model. In Section 3 we compute the spectrum of the latter model, which is then used in Section 4 to evaluate the chain's partition function. Using the grand canonical partition function of the spin dynamical model, in Section 5 we obtain an alternative expression for the chain's partition function in terms of that of the su(m|n) PF chain of A N −1 type. This expression turns out to be particularly efficient for the computation of the chain's spectrum for a relatively large number of spins. Section 6 is devoted to generalizing to spin chains of BC N type the boson-fermion duality uncovered in Refs. [25,28]. In the last section we use the explicit formulas for the partition function obtained in Section 5 to study the asymptotic behavior of the level density and the spacings distribution when the number of particles is very large. The paper also includes two appendices, the first one being a collection of several useful q-number identities, and the second a detailed derivation of the exact formulas for the mean and standard deviation of the chain's energy used in Section 7.
where α ij (s) = π(s i )π(s j ) + π(s i ) + π(s j ) f ij (s) , is the number of fermions in between the i-th and j-th spins. In other words, the sign (−1) α ij (s) is −1 when either both s i and s j are fermionic, or s i and s j are spins of different type with an odd number of fermionic spins between them. An equivalent way of defining the action of S (m|n) ij is by requiring that S (m|n) i,i+1 introduce a factor of −1 precisely when s i and s i+1 are fermionic. This clearly implies the previous (apparently more general) rule, since an arbitrary permutation is a product of permutations of consecutive elements.
The Hamiltonian of the supersymmetric su(m|n) PF chain of BC N type is defined by 1 whereS and β > 0. The chain sites ξ i are the same as those of the ordinary PF chain of BC N type studied in Ref. [33], namely the coordinates of the (unique) equilibrium of the scalar potential in the open set C = {x ≡ (x 1 , . . . , x N ) | 0 < x 1 < · · · < x N }. It can be shown [32] that ξ i = √ 2y i , where y i is the i-th zero of the Laguerre polynomial L β−1 N . Note also that the convention for the sign of the SUSY spin permutation operators in the definition of H (m|n) ǫǫ ′ is the opposite to that of Ref. [27], but (as we shall see in a moment) our choice is more consistent with the use of the names "bosonic" and "fermionic" for the two species of spins. Since where S i is the operator reversing the i-th spin, the purely bosonic (resp. purely fermionic) spin chain H (m|0) ǫ (resp. H (0|n) ǫ ) coincides with the ordinary ferromagnetic (resp. antiferromagnetic) PF chain of BC N type [33]. 1 Note that the operatorsS (m|n) ij actually depend on ǫ and ǫ ′ , although for simplicity we have chosen not to make this dependence explicit in the notation.
The spin chain (5) is naturally related to the su(m|n) spin dynamical model where x ± ij = x i ± x j , b = βa and a > 1/2, and to its scalar counterpart whereH (m|n) ǫǫ ′ is obtained from the spin chain Hamiltonian (5) by replacing the sites ξ i by the coordinates x i . Since as a grows to infinity the particles tend to concentrate at the chain sites ξ i , so that the dynamical and internal degrees of freedom decouple. By Eq. (10), in this limit the energies of the dynamical spin model (8) are approximately given by where E sc i and E j are any two energies of the scalar Hamiltonian (9) and the spin chain (5). The latter formula cannot be used directly to deduce the spectrum of H relating the partition functions Z and H sc , respectively. The partition function Z sc of the scalar model (9) was computed in Ref. [33]. We shall see in the next section that the spectrum of H (m|n) ǫǫ ′ is quite simple, a fact that shall be exploited in Section 4 to explicitly compute its partition function. Using Eq. (12) we shall then obtain the partition function Z (m|n) ǫǫ ′ in closed form.

Spectrum of the spin dynamical model
The first step in the evaluation of the partition function of the spin chain (5) by means of the freezing trick relation (12) consists in determining the spectrum of the spin dynamical model (8). To this end, we introduce the auxiliary "exchange Hamiltonian" where K ij is the operator permuting the i-th and j-th spatial coordinates, K i reverses the sign of the i-th coordinate, andK ij ≡ K i K j K ij . The rationale for introducing this operator is the fact that its restriction to states symmetric under the action of the operators coincides with that of the Hamiltonian H (m|n) ǫǫ ′ , as we shall next discuss. To begin with, note that all three sets of operators } generate a realization of the Weyl group of BC N type, namely they satisfy the commutation relations (where i, j, k, l are distinct) and similarly for the other two sets. We shall denote by Λ In particular, Λ (m|0) ǫǫ (resp. Λ (0|n) ǫǫ ) projects onto states totally symmetric (resp. antisymmetric) under particle permutations and with parity ǫ with respect to simultaneous sign reversals of spatial coordinates and spins. From the definition of the total symmetrizer Λ Hence H which, as we shall see, is precisely the relation needed to deduce the spectrum of H Indeed, recall [33] that the operator H ′ is represented by an upper triangular matrix in the (non-orthonormal) basis with elements where where the eigenvalues are given by Due to the impenetrable nature of the singularities of the Hamiltonian (8), its Hilbert space is the set L 2 0 (C) ⊗ Σ (m|n) of spin functions square integrable on the open set C and vanishing sufficiently fast on the singular hyperplanes form a basis of this Hilbert space provided that the following conditions hold: The last two conditions, which differ from the analogous ones for the su(m) case [33], deserve further remark. Condition (ii.a) simply states that acting with appropriate generalized spin permutation operators S (m|n) ij we can first reorder the spins so that within each sector of s (spin components corresponding to a maximal sequence of equal components of k) bosons precede fermions. Furthermore (condition (ii.b)), the bosonic spin components can be arranged in non-increasing order, while for the fermionic ones we can enforce strictly decreasing order since in this case the state vanishes if s i = s j by antisymmetry. (This is, by the way, the justification for the names "bosonic" and "fermionic" used for the two species of spins.) Finally, condition (iii) stems from the fact that the spin functions (19) A direct consequence of this condition is that the number of bosonic spin values of s j compatible with rule (iii) is given by Similarly, the number n(k j ) of fermionic spin values of s j compatible with the third rule is given by the previous expression, with m replaced by n and ǫ by ǫ ′ .
The Hamiltonian H is represented by an upper triangular matrix in the basis (19), ordered according to the degree |k|. Indeed, taking into account that H ′ commutes with Λ Hence the eigenvalues of H (m|n) ǫǫ ′ are given by where k and s satisfy conditions (i)-(iii) above. Note that, although the value of E k,s is independent of s, the degeneracy of each level clearly depends on the spin through the latter conditions.
From Eqs. (7)-(9) it follows that the Hamiltonian H sc of the scalar Calogero model is simply H (1|0) + , so that its spectrum is also given by the RHS of Eq. (21), with all k i even by condition (iii). We can thus subtract the ground energy E 0 from the spectra of both H sc and H (m|n) ǫǫ ′ without altering the partition function of the chain (5). With this proviso, the energies of both models are proportional to the coupling constant a, which implies that Z sc (aT ) and Z (m|n) ǫǫ ′ (aT ) are independent of a. Hence, in the calculation of the partition function of the chain (5) we can take without loss of generality E 0 = 0 and a = 1, and write the freezing trick formula (12) simply as

Evaluation of the partition function
We are now ready to compute the partition function of the spin chain (5) using the freezing trick formula (22). As for the ordinary (purely bosonic or fermionic) chain, the evaluation of the partition function Z (m|n) ǫǫ ′ of the spin dynamical model (8) depends on the parity of the integers m and n, so that one has to consider four different cases. In the calculations that follow, we shall use the fact that the number of ways into which we can arrange the spin components s i in a sector of s corresponding to a certain component κ of k repeated ν times is given by Clearly, this number depends on κ only through its parity, cf. (20) and its fermionic analogue. Note also that when m(κ) = 0 the second binomial coefficient in (23) should be interpreted as δ iν .

Case 1: m, n even
Recall, to begin with, that after setting E 0 = 0 and a = 1 the partition function of the scalar model (9) reads [33] Z sc (T ) = By Eq. (21), the partition function of the Hamiltonian (8) can be written as where the spin degeneracy factor d k is the number of spin states |s satisfying conditions (ii) and (iii) for a given multi-index k. Since both m and n are even, by Eq. (20) and its fermionic analogue we have m(k j ) = m/2 and n(k j ) = n/2. If where ν = (ν 1 , . . . , ν r ) and Since r i=1 ν i = N, the multi-index ν is an element of the set P N of partitions of N with order taken into account. Inserting Eqs. (26)- (28) into (25) we The RHS of the previous equation can be evaluated as in [33,Eq. (24)], with the result where ℓ(ν) = r is the number of components of the multi-index ν and From Eqs. (22), (24) and (29) it follows that in this case the partition function of the su(m|n) PF chain of BC N type is given by where the positive integers N ′ j are defined by Note that, as for the purely bosonic/fermionic model, in this case the partition function does not depend on ǫ, ǫ ′ . We shall therefore drop from now on the subscripts ǫ, ǫ ′ from Z (m|n) ǫǫ ′ when both m and n are even. The first product in Eq. (30) is simply the partition function of the su(2|0) chain (5) (cf. Eq. (34) in Ref. [33]), whereas the remaining sum coincides with the partition function Z of the su( m 2 | n 2 ) PF chain of type A N −1 computed in [29]. We thus obtain the remarkable relation Recall [25] that the partition function of the su(k|l) PF chain of type A is explicitly given by where the notation (q) k is defined in Eq.
Case 2: m odd, n even In this case, by Eq. (20) and its fermionic counterpart we have n(k j ) = n/2, while m(k j ) depends on the parity of k j . Hence it is convenient to modify condition (i) by first grouping separately the components of k with the same parity and then ordering separately the even and odd components. In other words, we shall now write k = (k e , k o ), where (23) it follows that the spin degeneracy factor d k is now given by where Equation (25) for the partition function of the dynamical spin model (8) then becomes where the RHS can be evaluated as in [33,Eq. (29)]. We thus obtain Substituting Eqs. (24) and (36) into (22) we arrive at the following explicit expression for the partition function of the spin chain (5) for odd m and even n: Note that in this case the partition function depends on ǫ but not on ǫ ′ .
As is the case for any finite system with integer energies, the partition function Z should be a polynomial in q. We shall now prove this fact by simplifying Eq. (37) with the help of the identities in Appendix A. To this end, we define the two sets of integers in terms of which we can rewrite Eq. (37) as From the discussion in Appendix A (cf. Eqs. (A.6)-(A.7)), it follows that the binomial coefficient N Ns q 2 is an even polynomial of degree 2N s (N − N s ) in q. Thus the RHS of Eq. (38) is also a polynomial in q, as expected.

Case 3: m even, n odd
The evaluation of the partition function Z (m|n) ǫǫ ′ is performed as in the previous case, the only difference being that the spin degeneracy factor d s ǫ (ν) in Eq. (38) should be replaced by where now In particular, the partition function Z (m|n) ǫǫ ′ does not depend on ǫ in this case.

Case 4: m, n odd
In this case both m(k j ) and n(k j ) depend on the parity of k j , so that the computation of the partition function proceeds as in Case 2. The final result is still Eq. (38), with the spin degeneracy factor d s ǫ (ν) replaced by where Thus in this case the partition function depends on both ǫ and ǫ ′ .

Relation with the spin chain of A N −1 type
We have seen in the previous section that when the integers m and n are both even, the partition function Z (m|n) of the chain (5) is related in a simple way to that of the su( m 2 | n 2 ) PF chain of A N −1 type, cf. Eq. (31). The purpose of this section is to derive an analogous relation for odd values of m and/or n. To this end, we shall deduce an alternative expression for the partition function by means of the grand canonical partition function of the spin dynamical model (8), following an approach similar to that of Refs. [21,25] for the ordinary and supersymmetric PF chains of A N −1 type.
In order to compute the grand canonical partition function, it is convenient to replace conditions (i)-(iii) in Section 3 ordering the labels k, s of the basis states (19) by the following equivalent set of rules: (a) If i < j then π(s i ) π(s j ) (b) If i < j and π(s i ) = π(s j ), then s 1 By Eq. (21) (setting, as before, E 0 = 0 and a = 1) the grand canonical partition function Z (m|n) ǫǫ ′ of the spin dynamical model (8) is given by where µ is the chemical potential and the inner sum runs over all values of k ∈ N N 0 and s compatible with conditions (a)-(d). As in the previous section, the computation of Z (m|n) ǫǫ ′ depends on the parity of m and n. Since the case of even m and n has already been dealt with above, we shall start with the case of m odd and n even: where ν i ∈ N 0 and m+ñ i=0 ν i = N. With a slight abuse of notation, we have suppressed the second component s 2 i of each spin s i , since clearly s 2 i = 0 for the first ν 0 +· · ·+νm spins and s 2 i = 1 for the remaining ones. By condition (c), the multi-index k is of the form where and the nonnegative integers κ i j satisfy By Eq. (43), the grand canonical partition function of the chain (5) is given by The first factor in the last equality is clearly the grand canonical partition function Z (1|0) ǫ of the purely bosonic spinless dynamical model (8) (which coincides with the scalar model (9) for ǫ = 1). On the other hand, in Ref. [25] it was shown that Using again the results in Ref. [25], we can write the previous formula as where Z denotes the grand canonical partition function of the su( m−1 2 | n 2 ) spin Calogero model of type A. On the other hand, the grand canonical partition function Z(T, µ) can be expressed in terms of the N-particle canonical partition function Z N (T ) as  22), the partition function of the N-particle su(m|n) PF chain of BC N type is thus given by The partition function Z (1|0) +,M is simply given by where we have used (24). On the other hand, from conditions (i) and (iii) in Section 3 we have The partition function of the su( m−1 2 | n 2 ) spin chain of type A with N − M particles is related to the partition function of its corresponding spin dynamical model by [25] Inserting Eqs. (49)-(51) into (48) we finally obtain where we have used the identity (A.2). Alternatively, from the first equality in (52) and Eq. (31) we can also derive the relation

Case 3: m even, n odd
This case is very similar to the previous one, the main difference being that Eq. (47) now becomes which leads to the following expression for the chain's partition function: By conditions (i)-(iii) in Section 3, the partition function Z (0|1) +,M of the spinless fermionic model (8) with even parity is given by Similarly, Inserting (55) and (56) into Eq. (54) and proceeding as before, we obtain the following two expressions for the partition function of the chain (5) with even m and odd n:

Case 4: m, n odd
In this case both the bosonic and fermionic spin components can take the zero value, so that instead of Eq. (47) we now have . (58) In particular, setting m = n = 1 in the previous this formula we have and therefore Z .
(60) Proceeding again as in Case 2, we easily obtain As before, this expression can be simplified in two alternative ways, leading to the following remarkable identities for the partition function of the chain (5) with odd m and n: where p is given by (46) and p ′ = (1 − ǫ ′ )/2. We thus have where we have used Eq. (A.3) with x = 1 and q replaced by q 2 . On the other hand, if ǫ = ǫ ′ = −1 Eq. (63) becomes Next, if ǫ = −ǫ ′ = 1 we have 6 Boson-fermion duality In Ref. [28] it was uncovered a remarkable boson-fermion duality satisfied by the spectrum of any spin chain of A N −1 type with global su(m|n) symmetry.
We shall now extend these ideas to the PF chain of BC N type (5), although our results will in fact be valid for any chain of BC N type possessing global su(m|n) symmetry.
Let us start by defining the star operator * : Σ (m|n) → Σ (m|n) by and the exchange operator X (m|n) : Σ (m|n) → Σ (n|m) by In other words, s ′ i has the same value as s i but opposite type. Following Ref. [28], we next define the operator U (m|n) : Σ (m|n) → Σ (n|m) as the composition Since, clearly, the operators * • * and X (n|m) • X (m|n) are the identity in Σ (m|n) , U (m|n) is invertible and Indeed, if |s ∈ Σ (n|m) , we have Moreover, it is also easy to see that the star operator is self-adjoint and X (m|n) † = X (n|m) , so that U (m|n) † = * • X (n|m) = U (m|n) −1 , i.e., U (m|n) is unitary. A straightforward computation [28] shows that U (m|n) S (71) The above considerations are enough for proving duality in the A N case. However, in the BC N case we need also consider the behavior of the spin reversal operators S ǫǫ ′ i with respect to conjugation by U (m|n) . For clarity's sake, in the following discussion we shall use the more precise notation S and therefore With the help of Eqs. (71) and (72) one easily obtains where the last sum was evaluated in Ref. [33]. Since U (m|n) −1 H (n|m) ǫ ′ ǫ U (m|n) and H (n|m) ǫ ′ ǫ are isospectral, from the previous equation we obtain the remarkable duality relation Z As an example, if m = n = 2 this property can be readily verified using the explicit formula Eq. (33).
As mentioned above, a duality relation of the form (73) clearly holds for an arbitrary chain with global su(m|n) symmetry, with Hamiltonian where c ij ,c ij and c i are real constants. More precisely, from Eqs. (71) and (72) it now follows that

Level density and spacings distribution
Having derived several explicit formulas for the partition function of the spin chain (5), it is natural to inquire on the global properties of its spectrum. In fact, the partition function (and hence the spectrum) can be computed in a very efficient way for relatively large values of N by expanding in powers of q Eqs. In the first place, our calculations for a wide range of values of N, m and n show that the energy levels of the chain (5) are equidistant, as already observed by us in the purely bosonic or fermionic case [33]. More precisely, we have found that the distance δE between consecutive levels is equal to one, except for the su(1|1) chain with ǫ = ǫ ′ , for which δE = 2 by Eqs. (64) and (65).
Secondly, our numerical calculations evidence that for sufficiently large N (N 20) the normalized level density where E 1 < · · · < E L are the distinct energy levels and d i is the degeneracy of E i , can be approximated with great accuracy by the Gaussian law with parameters µ and σ given by the mean and standard deviation of the chain's spectrum. Since the energy levels are equidistant, this means that As an illustration, in Fig. 1 we have plotted both sides of the latter equation in the case m = n = 1 and ǫ ′ = −ǫ = 1 for N = 15 and N = 30. The approximate Gaussian character of the level density is thus a property shared by all spin chains of HS type studied so far, both of type A and B [20,27,31,33,34]. Since, by Eq. (76), the level density is asymptotically given by the Gaussian law (75), it is of interest to compute the mean energy and its variance as functions of N, m, n, ǫ and ǫ ′ , as has been done for other spin chains of HS type (including the purely fermionic version of the chain (5)). Indeed, in Appendix B we show that µ and σ 2 are given by where p m , p n ∈ {0, 1} are the parities of m and n, respectively. In particular, when N tends to infinity µ and σ 2 grow as N 2 and N 3 , as for the ordinary (non-supersymmetric) PF chains of type A and B (cf. Refs. [34,33]).
Another property of the chain's spectrum worth studying is the distribution of spacings between consecutive "unfolded" levels. In general [38], the unfolding of the levels E i of a spectrum is the mapping is the continuous part of the cumulative level density It can be easily shown that the unfolded spectrum {η i } L i=1 is uniformly distributed regardless of the initial level density, making it meaningful to compare different spectra. In our case, by the previous discussion we can take η(E) as the cumulative Gaussian density (75), namely One then defines the normalized spacings is the mean spacing of the unfolded energies, so that {s i } L−1 i=1 has unit mean.
As mentioned in the Introduction, the main motivation for studying the spacings density p(s) lies in the Berry-Tabor conjecture, which states that the distribution of spacings of a "generic" quantum integrable system should obey Poisson's law p(s) = e −s . By contrast, for a chaotic system the spacings distribution is expected to follow Wigner's surmise p(s) = (πs/2) exp(−πs 2 /4), as is approximately the case for the Gaussian ensembles in random matrix theory [38]. A detailed study of the spacings distribution has been performed for many spin chains of HS type, namely the PF chains of types A and B [34,33], the original (type A) HS chain [20,34], as well as the supersymmetric version of the latter chain [27]. The spacings distributions of all these chains, which turn out to be qualitatively very similar, differ essentially from both Wigner's and Poisson's distributions. More precisely, for the ordinary (non-supersymmetric) chains mentioned above we showed in our recent papers [33,34] that the cumulative spacings distribution P (s) ≡ s 0 p(s ′ ) ds ′ is approximately given by where s max is the maximum spacing. As a matter of fact, in Ref. [33] we proved that the previous approximation holds for any spectrum E min ≡ E 1 < · · · < E L ≡ E max , provided that the following conditions are satisfied: (i) The energies are equispaced, i.e., E i+1 − E i = δE for i = 1, . . . , L − 1.
(ii) The level density (normalized to unity) is approximately given by the Gaussian law (75).
(iv) E min and E max are approximately symmetric with respect to µ, namely Moreover, when these conditions are satisfied the maximum spacing can be estimated with great accuracy as It should also be noted that Eq. (80) is valid only for spacings s s 0 , where s 0 is the unique zero of the RHS of this equation, namely Since, by Eq. (81) and condition (iii), we have (for instance) it follows that s 0 ≪ s max . The validity of the above conditions for the chain (5) when N ≫ 1 was checked in Ref. [33] in the purely fermionic case, which automatically implies their fulfillment in the purely bosonic case on account of the duality relation (73). For this reason, we shall assume in the rest of this section that m, n 1. We shall see below that in this case conditions (i)-(iii) are also satisfied when N ≫ 1, while condition (iv) only holds for m = n. It turns out, however, that when one drops condition (iv) the approximate formula (80) still holds, albeit in a slightly smaller range.
Indeed, it was proved in Ref. [33] that conditions (i)-(iii) imply that the spacings s i are approximately related to the energies E i by Hence where are the roots of the equation s = s max e − (E−µ) 2 2σ 2 , cf. Fig. 2. Using condition (i) to estimate the RHS of the latter equation we easily obtain If Since (again by condition (i)) E max − E min = (L − 1)δE, substituting Eq. (85) into (88) and using (81) we immediately arrive at the approximation (80) for P (s). We have thus shown that conditions (i)-(iii) imply that Eq. (80) holds for s s min . Note, finally, that the minimum spacing s min satisfies the inequalities s 0 s min ≪ s max .
Indeed, the second inequality is an immediate consequence of Eq. (87) and condition (iii). As to the first one, by Eqs. (81), (82) and (87) we have where the equality holds if and only if E min and E max are symmetric about µ. Let us next verify that the spectrum of the chain (5) satisfies conditions (i)-(iii) when N ≫ 1 and m, n 1. We have already seen at the beginning of this section that conditions (i) and (ii) hold for sufficiently large N. In order to check the validity of condition (iii), we need to compute in closed form the chain's maximum and minimum energies. By equation (11), the minimum energy is given by where E sc min and E min are the minimum energies of the scalar and spin dynamical models (9) and (8), respectively. From the discussion in Section 3 it follows that E sc min = E 0 , so that (cf. Eq. (21)) E min is the minimum value of |k|, where k is any multi-index compatible with conditions (i)-(iii) in Section 3. If m > 1 we can obviously take k = 0 (and, for instance, s = (s, 0), . . . , (s, 0) , where s is any positive spin component), so that E min = 0 in this case. Similarly, when m = ǫ = 1 the minimum energy is again zero, since the pair k = 0 and s = (0, 0), . . . , (0, 0) satisfies the required conditions. However, if m = 1 and ǫ = −1, condition (iii) implies that k i must be odd when s 1 i = 0. Hence, the multi-index k yielding the minimum energy is of the form As to the maximum energy, the duality relation (73) and the previous equation imply that so that condition (iv) is not satisfied unless m = n, as claimed above.
We have verified that the approximation (80)-(81) is in excellent agreement with the numerical data obtained using our exact formulas for the partition function for many different values of m, n 1 and N 15. For instance, in Fig. 3 we plot the cumulative spacings distribution P (s) versus its analytic approximation (80) in the cases m = n = −ǫ = ǫ ′ = 1 and m = 2, n = ǫ = ǫ ′ = 1 for N = 15 and N = 30 spins. The respective root mean square errors (normalized to the mean) drop from 9.4 × 10 −3 and 4.6 × 10 −3 for N = 15, to 2.4 × 10 −3 and 1.1 × 10 −3 for N = 30. It is worth noting that the approximation (80)-(81) for the cumulative spacings distribution contains no adjustable parameters, since the maximum spacing s max is completely determined by Eqs. (78), (89) and (90). In fact, from the previous equations one easily obtains the following asymptotic formula for the maximum spacing in the genuinely supersymmetric case m, n 1: Thus the behavior of s max for large N is qualitatively the same as in the non-supersymmetric cases studied in Ref. [33]. A Some useful q-number identities In this appendix we shall collect several identities involving q-numbers that are used in the simplification of the partition function Z (m|n) ǫǫ ′ of the chain (5) when either m or n are odd. Given a real number q ∈ (0, 1) and a nonnegative integer k, we define the symbols (q) k and [k] q by with (q) 0 ≡ 1. The q-factorial of k is then defined as Finally, if k, l are nonnegative integers with k l, we define the q-binomial coefficient From the above definitions it immediately follows the useful equality An important identity satisfied by the q-binomial coefficients is the following instance of Newton's q-binomial formula [39, Eq. (8)]: From the previous equation with q replaced by q 2 and x by q we obtain the formula Replacing q by q 2 and setting x = q −1 in Eq. (A.3) we obtain the analogous relation The last three identities are used in Section 5 to compute the partition function Z Likewise, if δ q is the q-dilation operator, whose action on a smooth function f : R → R is given by (δ q f )(x) = f (qx), we have [39,Eq. (11)] where 1 denotes the constant function identically equal to 1. Equating the coefficient of x i in both sides of the previous equation we obtain the remarkable identity where τ i (x 1 , . . . , x k ) denotes the elementary symmetric polynomial of degree i in k variables, defined by In particular, from Eq. (A.6) it easily follows that the q-binomial coefficient k i q is a polynomial of degree i(k − i) in q.

B Computation of the mean and variance of the energy
In this appendix we shall compute the mean and variance of the energy of the chain (5) for arbitrary values of N, m, n, ǫ and ǫ ′ . In the purely fermionic case, these quantities were computed in Ref. [33] using the formulas for the traces (of products) of the operators S ij and S i given in Ref. [31]. In the genuinely supersymmetric case, one should evaluate the traces of products involving the more general operators S , which are the only ones appearing in the chain of type A, were computed in detail in Ref. [27]. The calculation of the traces involving the spin reversal operator S ǫǫ ′ i is totally analogous, so that we shall limit ourselves to listing the final results in Table B.1. In this table p m , p n ∈ {0, 1} are the parities of m and n, respectively, and we have omitted the traces of products which reduce to the identity matrix or to one of the entries.  In the first place, calling h ij = (ξ i − ξ j ) −2 ,h ij = (ξ i + ξ j ) −2 , h i = β ξ −2 i , and using the formulas in Table B.1 for the traces of the spin operators we easily obtain The sums in the previous expression were evaluated in Ref. [33], with the result i =j (h ij +h ij ) = 1 2 N(N − 1) , Inserting these values into (B.1) we immediately arrive at Eq. (77) for the mean energy µ.
The evaluation of the standard deviation of the energy which is considerably more involved, is also performed in two steps. We begin by expanding Eq. (B.2) and simplifying the resulting expression using Eq. (B.1) and the trace formulas in Table B.1. After a long but straightforward calculation we obtain where the symbol