Can a spin chain relate combinatorics to number theory?

The Motzkin spin chain is a spin-$1$ frustration-free model introduced by Shor&Movassagh. The ground state is constructed by mapping random walks on the upper half of the square lattice to spin configurations. It has unusually large entanglement entropy [quantum fluctuations]. The ground state of the Motzkin chain can be analytically described by the Motzkin paths. There is no analytical description of the excited states. The model is not solvable. We simplify the model by removing one of the local equivalence moves of the Motzkin paths. The system becomes integrable [similar to the XXX spin chain]. We call it free Motzkin chain. From the point of view of quantum integrability, the model is special since its $R$-matrix does not have crossing unitarity. We solve the periodic free Motzkin chain by generalizing the functional Bethe Ansatz method. We construct a $T-Q$ relation with an additional parameter to formulate the energy spectrum. This new parameter is related to the roots of unity and can be described by the M\"obius function in number theory. We observe further patterns of number theory.


Introduction
Entanglement is a primary feature of quantum mechanics lacking in classical mechanics. It can be used for building revolutionary new technologies (e.g., quantum computer) and information processing. One way of measuring this quantity is by computing entanglement entropy (EE). In this regard quantum many-body systems described by local Hamiltonians gain importance. Restricting ourselves to 1D systems we know that gapped systems obey the area law [1,2]. The area law says that EE between two subsystems of a system is proportional to the area of the boundary between them. This is expected given the local nature of the interactions. It is widely believed that for "physically reasonable" quantum systems, the area law can be violated by at most a logarithm log(N ) factor [3], with N the size of the quantum system. This belief has been seriously challenged by both quantum information and condensed matter theorists in recent years. It is of interest to find local Hamiltonians that exhibit more quantum fluctuations with their EE scaling as either the logarithm or the volume of the system size [4,5].
More recently, a spin-1 local, frustration-free Hamiltonian called the Motzkin spin chain [5][6][7] (Its half-integer analog is called the Fredkin spin chain [8,9]) was shown to have unique ground states. The model is constructed based on random walks, called the Motzkin walks. Each Motzkin walk (path) can be mapped to a wave function of spins. The Hamiltonian of the Motzkin chain is the nearest neighbors with projectors as interactions. It has unique ground state that can be seen as a uniform superposition of Motzkin walks. As the volume of the system N → ∞, the random Motzkin walks converge to Brownian excursions [10]. The EE of the system is proportional to √ N . It is exponentially more than the logarithm 1 . * haoke72@163.com † vladimir.korepin@stonybrook.edu 1 Actually, there is an infinite class of integer spin chain models that are physically reasonable and exact calculation of the EE shows that they violate the area law to the leading order by √ N .
Such a spin chain can in principle be experimentally realized. Hydrodynamics of the model was studied in [11]. The large amount of entanglement may be used as a resource for quantum technologies and computation. The model has since then been deformed to show phase transitions [12,13], generalized to possess new features using symmetric inverse semigroups [14,15] and simplified to preserve translational invariance in [16]. The ground state of the Motzkin chain can be analytically expressed as a matrix product state of the Motzkin walks using techniques from enumerative combinatorics [17]. While there is no analytical description of the excited states, the model is not exactly solvable.
The full Hamiltonian of the uncolored 2 Motzkin chain was described by Bravyi et al. [5]. It is a spin chain with a 3-dimensional local Hilbert space given by a basis |u , |f , |d . Here u, f , and d are used to abbreviate "up", "flat", and "down" of the Motzkin walks respectively. This identification comes from that the states in this system can be mapped to paths (Motzkin walks) in the "x-y" plane as done in [6]. The Hamiltonian of the Motzkin spin chain with periodic boundary conditions is given by (the coupling constant g was equal to 1 in [6]): The local Hamiltonians H f ree and H int are given in terms of two-site projectors, and The term H int can be viewed as a "kinetic part" of the Motzkin model. By turning on g, the system loses its symmetries and starts to "interact". Only the ground state of the chain can be analytically described by the Motzkin paths. The model is not exactly solvable. Thus we will call the case with g = 0 the non-interacting Motzkin chain or the free Motzkin spin chain. While it is still strongly interacting, it is indeed integrable in the sense of quantum inverse scattering method [18,19] (or the algebraic Bethe Ansatz method (ABA)), as expanded on below. This fact allows us to solve the spectrum of the excited states completely.
In this paper we will analyze the case of the free Motzkin spin chain We observe that the resulting system is integrable. Its structure can be described (generated) by Yang-Baxter equation (Eq.2.10). Then we will calculate its spectrum for periodic boundary conditions 3 . The model can be seen as a generalized Temperley-Lieb system. This is understood from the local Hamiltonian, which is a representation of the generator of Temperley-Lieb algebra. Due to the lack of commutation relations among the monodromy-matrix operators(Eq.2.16), it is hard to construct the Bethe state and further steps of the standard ABA. We generalize Sklyanin's functional Bethe Ansatz [21,22] and some ideas of Off-diagonal Bethe Ansatz [23,24] to approach the spectrum problem of the model 4 . This paper is organized as follows. In section 2, we introduce the integrable Hamiltonian and its associated integrable structure. In section 3, we briefly review the results of XXX spin-1 2 chain case. In section 4, in contrast with the XXX spin-1 2 chain case, we derive certain operator product identities for the transfer matrix of the model by using some intrinsic properties of the R-matrix. The asymptotic behavior of the transfer matrix are also obtained. Section 5 is devoted to the construction of the T − Q relation and the corresponding Bethe Ansatz equations. In section 6, we summarize our results and give some discussions. Some detailed technical derivation is given in Appendixes.

The Free Motzkin spin chain as an integrable model
To readers familiar with other similar models but not with the Bethe Ansatz method, the central building block is the R-matrix. It essentially describes the scattering between two quasiparticles in the spin chain. The R-matrix satisfies the Yang-Baxter equation, which states that many-body scattering can be reduced to two-body scattering. This fact ensures the integrability of the periodic free Motzkin chain.
The R-matrix R 12 (λ) = P 12 [(λ + η)1 − λê 12 ] for the free Motzkin spin chain can be written as Here λ is the spectral parameter. η is the crossing parameter. P 12 is the permutation operator,ê 12 is the corresponding generator [25] of the Temperley-Lieb algebra, and 1 2ê 12 is a projection operator.
Unless otherwise specified, we use the standard notation: for any matrix A ∈ End(V), A j is an embedding operator in the tensor space V ⊗ V ⊗ · · · . It acts as A on the j-th space and as an identity on the other factor spaces. So R ij (λ) acts as an identity on the factor spaces except for the i-th and j-th ones. We denote the basis states by {|u , |f , |d } where u, f and d are used to abbreviate "up", "flat" and "down" respectively, The operatorsÛ j,j+1 andD j,j+1 are projectors to the states |u j , f j+1 −|f j , u j+1 , and |d j , f j+1 −|f j , d j+1 respectively, and Alternatively, they also can be expressed in terms of standard basis and Pauli matrices like representation, In a similar manner, we have the expression for operatorD j,j+1 , The R-matrix satisfies the quantum Yang-Baxter equation (QYBE) and possesses the following properties: Initial condition: The partial transpose of this R-matrix is degenerate det(R t1 12 (λ)) = 0. (2.14) Therefore it does not satisfy crossing unitarity [26]. It can be seen as a special solution of Yang-Baxter equation [27]. This fact also makes the standard ABA procedure for Temperley-Lieb algebra [28] does not work on this model. The monodromy-matrix T (λ) can be constructed as where the subscript 0 indicates the auxiliary space and the indices {1, · · · , N } denote the physical or quantum spaces. N is the number of sites and {θ j |j = 1, · · · , N } are the inhomogeneous parameters.
Then the monodromy-matrix can be written as Here the subscripts of B, C, D operators indicate their position in auxiliary space 0. The transfer matrix t(λ) is defined as the partial trace of monodromy matrix in the auxiliary space Let us further introduce the vacuum state |0 , The operators acting on the vacuum state give rise to The QYBE (Eq.2.10) implies that the monodromy-matrix T (λ) satisfies the Yang-Baxter relation Therefore, one can prove that the transfer matrices with different spectral parameters commute with each other So t(λ) serves as the generating functional of the conserved quantities of the system. This ensures the integrability of the periodic free Motzkin spin chain. The Hamiltonian of periodic free Motzkin spin chain can be recovered from the transfer matrix through a logarithmic trace identity, such that Here the local Hamiltonian is with the periodic boundary conditions H N,N +1 = H N,1 .

A brief review of XXX spin-1 2 chain
The final results show that the spectrum of the free Motzkin chain can be constructed based on the XXX spin-1 2 chain. Here we take the XXX spin-1 2 chain as an example to further explain the situation in the free Motzkin chain. Following the concept of functional Bethe Ansatz or Separation of Variables (SoV) method, we have the operator identities for transfer matrix of XXX spin-1 2 chain 5 , Next, we consider the asymptotic behaviors of fused transfer matrix Then let us assume the eigenstate of t ( 1 2 ) (λ) is |ψ , with its eigenvalue Λ ( 1 2 ) (λ). Acting the above operator identities (Eq.3.1) and the asymptotic relation (Eq.3.2) on |ψ , they should satisfy the corresponding functional relations 6 As we know, the T − Q relation of the XXX spin-1 2 chain is a polynomial function of degree N . This can be seen from the construction of the monodromy matrix. The above N + 1 functional relations are sufficient to determine such a T − Q relation with Bethe equations The function Q(λ) is parameterized by N Bethe roots {λ j |j = 1, · · · , N } as follows 7 : To ensure Λ ( 1 2 ) (λ) to be a polynomial, the residues at the poles λ j must vanish. This requires the N Bethe roots to satisfy the Bethe equations

Operator identity of the transfer matrix
For the current free Motzkin chain model, it is difficult to implement the standard ABA method. Let us go along with the concept of functional BA [21,22] and consider the products of the transfer matrices at special points. To do this, we follow the method in [23] and take advantage of the symmetry and properties (Eq.2.11), (Eq.2.12) and (Eq.2.13) of the R-matrix. After a long calculation, we find that the transfer matrix (Eq.2.17) satisfies the following operator identities Other than the identity operator (matrix) on the r.h.s of (Eq.3.1) for the spin-1 2 chain case, here the op-eratorÔ is a special operator. Thus we can not simply convert the operator identities to the functional 5 The superscript ( 1 2 ) refers to spin-1 2 , and the functions a(λ) and d(λ) are the same as that of free Motzkin case. 6 People mainly focused on functional relations in previous works, since the transition between the operator identities and functional relations is straightforward. Also, the calculation of the functional relations is relatively easy. However, in our current model, we should emphasize the operator identities. See the next section for details. 7 We remark that the form of function Q(λ) here is the same as that of the free Motzkin case in a later section.
relations like in the XXX spin-1 2 chain and all other cases. The final results show thatÔ is a global operator acting highly non-trivial on all quantum spaces. One can not derive the functional relation by projecting the transfer matrix to its subspace. The fusion techniques [29] do not work here.
So for the current model, the key point is to understand the structure of theÔ operator. The form of operatorÔ consists of the leading terms of operators A(λ), D 11 (λ) and D 22 (λ) as polynomials of variable λ. It can be expressed as followŝ The second step of the formula above results from the leading terms of operator D 11 (λ) is proportional to the identity operator (matrix).
Next, let us briefly analyze the asymptotic behaviors of operators A(λ) and D 22 (λ) as λ tends to infinity. The operators A(λ) and D 22 (λ) are centrally symmetric to each other with respect to the element E (3 N +1)/2,(3 N +1)/2 (the center element of the matrix acting on the full quantum spaces). This means each element E i,j (we adopt the standard basis) in operator A(λ) can be mapped to the element E 3 N +1−i,3 N +1−j in operator D 22 (λ) after a 180-degree rotation around the central element. Here 3 N is the dimension of N quantum spaces.
The operators A(λ) and D 22 (λ) have the following asymptotic behaviors as λ goes to infinity. The leading terms of A(λ) operator is a degenerate operator acting on column vectors in the up half quantum space x 1 , x 2 , . . . , x M , 0, . . . , 0 T , and M = 3 N +1 2 . This degenerate operator can be seen as a series of row transformations acting on this diagonal matrix Here, from the first diagonal element E 1,1 to the M -th diagonal element E M,M , M elements each have the value 1, and the other elements are all 0.
Considering the symmetry between A(λ) and D 22 (λ), the leading terms of the composite operator A(λ) + D 22 (λ) is a non degenerate operator. It can be expressed as a series of row transformations acting on the diagonal matrix diag(1, · · · , 1, 2, 1, · · · , 1) = This is just an identity matrix with its center (M -th) diagonal element plus 1. Moreover, the asymptotic behavior of the transfer matrix is The operator identities (Eq.4.1) of the transfer matrix together with the asymptotic behavior (Eq.4.5) completely determine the spectrum of the transfer matrix. The detailed analysis of the structure of operatorÔ (also operators A(λ) and D 22 (λ)) and its eigenvalues are given in Appendix A.

Functional relations and the T − Q relation
The commutativity of the transfer matrices (Eq.2.21) means that they have common eigenstates. Let us assume that |Ψ is an eigenstate of t(λ). |Ψ does not depend upon λ, and its corresponding eigenvalue is Λ(λ), t(λ)|Ψ = Λ(λ)|Ψ . (5.1) The operator identities (Eq.4.1) imply that the corresponding eigenvalue Λ(λ) satisfies the relation Thus we have Here {ω} represents the eigenvalues of operatorÔ, which have 3 N numbers, such as The eigenvalues {ω} of operatorÔ include two parts: The diagonal elements result in one single 2 (the center element) and several 1 (with degeneracy of 2 N +1 − 2); The non-diagonal elements give rise to roots of unity, ranging from square roots to N -th roots.
In general, the distribution of k-th roots of unity as eigenvalues of operatorÔ for a lattice of length N is given by: where µ(n) refers to the Möbius function in number theory. C(n, k) is the inline notation of the binomial coefficient n k . For k = 1 this formula gives the correct diagonal result if one excludes the eigenvalue 2. The 1 with degeneracy of 2 N +1 − 2 from the diagonal elements can be regarded as the 1-st roots of unity. In every specific k-th roots case, each kind of roots exp( 2πi k ), exp( 4πi k ),· · · , 1 has the same degeneracy. See Appendix A and B for detailed derivation.
According to the asymptotic behavior of the transfer matrix (Eq.4.5), we have that the leading terms of Λ(λ) have the following values, They are λ N multiplied by each eigenvalue of operatorÔ plus 1, respectively. They have a one-to-one correspondence with the eigenvalues of operatorÔ (Eq.5.4). Λ(λ), as a function of λ, is a polynomial of degree N . The functional relations (Eq.5.3) together with this asymptotic behavior (Eq.5.6), these N +1 conditions are sufficient to construct the T −Q relation [30] for Λ(λ).
Due to {ω}, other than XXX spin-1 2 chain case, there must exist coefficients before the two terms of the T − Q relation(Eq.5.9). Assuming ω is one eigenvalue of operatorÔ, then the corresponding coefficient of the leading term of Λ(λ) as λ → ∞ should be ω + 1, as shown in (Eq.5.6). Letã andd represent the coefficients before a(λ) and d(λ) terms, respectively. Based on the identities (Eq.5.3) and the asymptotic behavior (Eq.4.5) and (Eq.5.6), these coefficients should satisfy the following equations ã ·d = ω, a +d = ω + 1. (5.7) The equations have two sets of solutions: After carefully calculations, we find that only the first solution set is appropriate 8 . Basically speaking, this is because {ω} are eigenvalues of the operatorÔ, and the operatorÔ can be related to A(λ)+D 22 (λ), which correspond to a(λ) (Eq.2.19). In contrast, the operator D 11 (λ) corresponds to d(λ), and D 11 (λ) is proportional to identity operator (matrix) under the limit λ tends to ∞.

T − Q Ansatz and the spectrum
To explain the situation in the free Motzkin chain, we have presented a quick review of how to apply the functional BA approach to the XXX spin-1 2 chain in Section 3. We propose a generalized T − Q relation to describe the eigenvalues of the transfer matrix (Eq.2.17). This T − Q relation, together with the newly appearedÔ operator (Eq.4.2) in the current model, allows us to calculate all spectrum of the Hamiltonian.
Based on the properties (Eq.5.3), (Eq.5.6), and also the configuration of the coefficient ω (Eq.5.8), let us propose the following conjecture for Λ(λ), Here ω will go through each value in the set {ω}. The function Q(λ) is parameterized by N Bethe roots {λ j |j = 1, · · · , N } as follows: To ensure Λ(λ) is a polynomial, the residues of Λ(λ) at the poles λ j must vanish. So the N Bethe roots satisfy the Bethe equations Then the energy spectrum of the Hamiltonian (Eq.2.22) can be given in terms of the associated Bethe roots by

Conclusions
In this paper, we have studied the free Motzkin chain with periodic boundary conditions. Based on intrinsic properties of the R-matrix, we derive the functional relations (Eq.5.3) of the transfer matrix. These relations, together with the asymptotic behavior (Eq.5.6), allow us to construct a generalized T − Q relation (Eq.5.9) for the eigenvalues of the transfer matrix and the associated Bethe equations (Eq.5.11). The new parameter {ω} in the T − Q relation is related to the roots of unity. Its distribution can be described by the Möbius function of the number theory. Some patterns in number theory can be observed from the distribution (see equation (Eq.5.5), also the Table 1 in Appendix). For example, when the lattice site N is a prime number, the total numbers of the square roots and the highest N -th roots are both 2 N − 2. Since the degeneracies of all kinds of N -th roots are equal, the number 2 N − 2 is divisible by N , for any prime N . The exactly solutions allow one to study the thermodynamic properties. While now the Bethe equations have new parameters {ω}, they have different values in different sub-spaces. One can not directly apply the Yang-Yang thermodynamic to this model. Our method may apply to other Temperley-Lieb models with crossing unitarity breaking 9 .

Appendix A: Structure of operatorÔ and its eigenvalues
OperatorÔ is composed of the leading terms of operators A(λ) and D 22 (λ) as λ → ∞. To study this, besides the local operators I u in (Eq.2.7) and I d in (Eq.2.9),we further introduce the following notations, For total lattice sites number N , let us set Here for B 2 (λ), the subscript 2 refers to the second B type operator in the main text (Eq.2.16). While the subscript of B N is the number of total system sites, and these operators act nontrivially on all the N quantum spaces. Then for the system with N lattice sites, Because of the symmetry between A(λ) and D 22 (λ), we only need to analyze the structure of operator A(λ). According to the structure of the monodromy matrix and matrix multiplication, we have the following recursive relations for these operators: Subscripts of I u , I d , S + and S − refer to the quantum spaces on which these local operators act. Based on the above recursive relations, we have the following rules for the structure of the operator A N : • The operator A N is composed of local operators I u , I d , S − and S + .
• A N is the sum of products of N local operators acting on each of the quantum spaces(from 1 to N , respectively).
• A N is the sum of all possible operator products permutations(terms).
• S − and S + operators appear in pairs, and S − always comes before S + in one permutation(term). That means the quantum subspace number of S − is smaller than that of S + .
• I u can appear on its own. It can also come before S − or after S + . I d only appears between S − and S + .
• Multiple S − S + pairs appear in turns. They can not be nested within each other.
Appendix A1: ForÔ 3 ,Ô 6 andÔ 7 cases, the structure and the rules of the eigenvalue distribution In the following, we take N = 3 as an example. A 3 can be represented as Considering the symmetry between A 3 and D 3 , the operatorÔ 3 = A 3 + D 3 has the matrix form (A.6) In this matrix, A 3 occupies the left 14 columns, and D 3 occupies the right 14 columns. The number 2 = 1 + 1 in the center (the entry in the 14-th row and 14-th column, 2 × E MM , M = 3 3 +1 2 = 14) is the sum of value 1 from A 3 and another 1 from D 3 .
In this matrix form, we divide the first quantum space into 3 × 3 blocks by solid lines. The block in the center, |f 1 f 1 | = E 2,2 1 (we adopt the standard basis, subscript refers to quantum space) is the operatorÔ 2 . The boldface 3 × 3 diagonal matrix diag{1, 2, 1} in the middle is just the operatorÔ 1 . Then the 2 in the center can also be represented as 2 × E MM = 2 × E 2,2 1 E 2,2 2 E 2,2 3 in quantum subspaces. We can analyze the distribution of eigenvalues ofÔ 3 from operator A 3 . The term I u 1 I u 2 I u 3 corresponds to all the diagonal elements. Considering the symmetry between A 3 and D 3 , the operator O 3 has 2 4 − 1 = 15 diagonal elements, including one single 2, and 14 elements equal to 1.
Square roots of unity are caused by S − S + pairs. Each term of A 3 expression consists of 3 local operators. In each term with a S − S + pair, the center element E 2,2 (in its own subspace) of the rest operator I u (or I d ) will lead to a square root. In A 3 operator, these elements are and S − 1 E 2,2 2 S + 3 . Considering D 3 , there are 6 elements related to square roots inÔ 3 . This number can be calculated by combination 10 C(3, 2) × 2 = 6. The 6 elements are marked with red color in the above matrix. These elements will lead to 6 square roots of unity: −1, −1, −1, 1, 1 and 1.
The S − S + pairs are also the basic components of higher order roots of unity. Cubic roots come from terms with one S − S + pair. Unlike the case of square roots, in terms with a S − S + pair (3 local operators in each term), the rest one operator I u (or I d ) should take E 1,1 (or E 3,3 ). This E 1,1 (or E 3,3 ) makes the order plus one: 2 + 1 = 3. Then in A 3 operator, these elements are E 1,1 1 S − 2 S + 3 , S − 1 S + 2 E 1,1 3 and S − 1 E 3,3 2 S + 3 . Considering D 3 , there are 6 elements inÔ 3 . We calculate the number as C(3, 2) × C(1, 1) × 2 = 6. As marked in blue in theÔ 3 matrix form, the 6 elements will lead to 6 cubic roots: exp( 2πi 3 ), exp( 2πi 3 ), exp( 4πi 3 ), exp( 4πi 3 ), 1 and 1. Next let us take operator A 6 as an example to explain more rules of eigenvalues ofÔ. The operator Among them, the term I u 1 I u 2 I u 3 I u 4 I u 5 I u 6 corresponds to all the diagonal elements. Similarly, together with D 6 , there are 2 7 − 1 = 127 diagonal elements inÔ 6 , and one of them is 2. All the rest of the elements are 1.
The number of square roots of unity is (C(6, 2) + C(6, 4) + C(6, 6)) × 2 = 2 6 − 2 = 62, and half (62/2 = 31) of them are −1, the other half are 1. Single S − S + pair in one term leads to square roots, when all the rest operators take E 2,2 . Recurring S − S + pairs (2 or 3) in one term with the rest operators taking E 2,2 , will also give rise to square roots. This means we should also take terms with S − S + S − S + and S − S + S − S + S − S + into account.
The number of cubic roots is (C(6, 2) × C(4, 1) + C(3, 1)) × 2 = 126. The three kinds of cubic roots are exp( 2πi 3 ), exp( 4πi 3 ) and 1, and each has degeneracy 126 3 . In the expression of operator A 6 , let us consider terms with one S − S + pair and four I u (or I d ) operators. When one I u (or I d ) operator takes E 1,1 (or E 3,3 ) in its subspace, and the other three I u take E 2,2 , each such choice in one term will lead to a cubic root. And recurring I u S − S + (I u takes E 11 ), S − I d S + (I d takes E 33 ) and S − S + I u (I u takes E 11 ) patterns in one term also lead to cubic roots.
The number of 4-th roots is (C(6, 2) × C(4, 2)) × 2 = 180. The four kinds of 4-th roots are i, −1, −i and 1, and each kind has degeneracy 180 4 . 4-th roots come from terms with a S − S + pair and two I u (or I d ) taking E 1,1 (or E 3,3 ) in its subspace, and the other two I u (or I d ) taking E 2,2 . Two S − S + pairs in one term with the rest two operators taking E 2,2 , will reduce to square roots, as we mentioned above.
The second consequence is that the operator maps basis states to basis states. Apart from the state |f ⊗ |f ⊗ ... ⊗ |f ⊗ |f which is mapped to itself by both A N and D N and so is an eigenstate with eigenvalue 2, the operator acts as a permutation matrix on the basis states. It acts as a left cyclic shift of all sites not equal to |f . For example, on a chain of length 3 it maps |u ⊗ |d ⊗ |u to |d ⊗ |u ⊗ |u . On a chain of length 6, the state |u ⊗ |u ⊗ |d ⊗ |f ⊗ |u ⊗ |d can be mapped to |u ⊗ |d ⊗ |u ⊗ |f ⊗ |d ⊗ |u .
As a result, finding the spectrum ofÔ is equivalent to finding the cycle decomposition of the corresponding permutation. This is also equivalent to the combination rules in the last section. We can write down the distribution of eigenvalues for arbitraryÔ N .
The number of k-th roots for a lattice of length N is given by: where µ(n) refers to the Möbius function. For k = 1 this formula gives the correct diagonal result if one excludes the eigenvalue 2.
Some results in number theory can be observed from the distribution. For example, when N is a prime number, the total number of square roots is 2 N − 2. It is equal to the total number of the highest N -th roots in this case. As we mentioned above, the degeneracies of all kinds of N -th roots are equal. This will lead to a result in number theory: 2 N − 2 is divisible by N , for N is prime.
This was computed with the following implementation in Mathematica: