Composing arbitrarily many $SU(N)$ fundamentals

We compute the multiplicity of the irreducible representations in the decomposition of the tensor product of an arbitrary number $n$ of fundamental representations of $SU(N)$, and we identify a duality in the representation content of this decomposition. Our method utilizes the mapping of the representations of $SU(N)$ to the states of free fermions on the circle, and can be viewed as a random walk on a multidimensional lattice. We also derive the large-$n$ limit and the response of the system to an external non-abelian magnetic field. These results can be used to study the phase properties of non-abelian ferromagnets and to take various scaling limits.


Introduction
Symmetries in physics and in mathematics are often encoded in unitary groups. In quantum mechanics, in particular, many physical systems realize representations of unitary groups. The canonical example is spin, an irreducible representation of the group SU(2) (a double cover of the group of spatial rotations SO(3)), with isospin (also SU(2)) and elementary particle "flavor" symmetry (SU(3)) constituting more exotic examples. Systems consisting of several components, each carrying an irreducible representation of a symmetry group, transform under the same symmetry, in the di-rect sum of the representations of the components. The decomposition of the states of the full system into irreducible representations of the symmetry group becomes, then, a problem of physical relevance.
The group theoretical techniques of decomposing a direct sum of two irreducible representations (irreps) are well established [1]. In particular, the irrep content of the decomposition of the tensor product of irreps of SU(N) has been studied in the mathematics literature [2,3]. The result for the tensor product of fundamental irreps, in particular, can be derived from Schur-Weyl duality [4]. 1 The existing derivations, however, use a rather heavily mathematical language and the results are in a somewhat implicit form. A derivation in a physically motivated language and in a way accessible to physicists was desirable. In [5] we undertook this task for the simplest case of SU (2). Specifically, we computed the multiplicity of spin j irreps arising from composing n spins s (related results were also derived in [6,8] and used in [7], and were further elaborated in [9]). We also derived explicit expressions for the multiplicities arising from spins corresponding to indistinguishable particles by distilling the totally symmetric (for bosons) or antisymmetric (for fermions) parts of the decomposition. We could recast the composition problem as random walks on a one-dimensional lattice, or equivalently generalized Dyck or Lukasiewicz paths [10] with appropriate boundary conditions. Even for the SU(2) case the results are interesting, exhibiting nontrivial scaling properties in the large-n and large-s limits, especially for bosonic or fermionic spins. In [5] we provided the asymptotic expressions in these limits, and implicit. In the present work we simplify somewhat the problem by restricting to the direct sum of n fundamental representations of SU(N). As we shall see, this problem still has a rich structure.
Our motivation for carrying out this analysis is manifold. Firstly, we wished to present a complete and explicit solution to the problem. Moreover, we wanted to present an intuitive and pedagogical approach, building on physics connections and intuition. In this vein, we used the fact that this problem can be related to a system of N free quantum mechanical fermions on a circle. Finally, systems of many SU(N) degrees of freedom arise in condensed matter and ultracold systems [11,12], spin chains [13,14], particle physics [15,16], matrix models [17], etc., and we wanted to pave the road for possible physical applications.
In this paper we present a full solution to the problem of decomposing the direct sum of an arbitrary number of SU(N) irreps and apply it to the case of fundamental irreps, deriving explicit formulae for the representation content of the decomposition.
The example of SU(3) is analyzed in some detail. We point out a random walk connection, this time on a N − 1-dimensional lattice, and identify a duality relation in the distribution of irreps. We also consider the limit of a large number n of fundamentals, as a prelude to more general and intricate limiting cases. Finally, we calculate the response of the system to coupling it with a nonabelian magnetic field, which is relevant in the study of the thermodynamics of interacting SU(N) systems. The structure and phase transitions of the thermodynamic limit of such systems are treated in [20], and nontrivial scaling limits involving both N and n taken large will be investigated in forthcoming publications.
We present our results in the upcoming sections, hinting at additional physical applications and possible extensions of our analysis in the conclusions.

Fermion representation of SU(N) group theory
We review here the description of irreducible representations (irreps) of SU(N) as Nfermion energy eigenstates on the circle [18] and the corresponding composition rules in the fermion picture. This representation affords a conceptual standpoint that is well-suited to our considerations and will be useful throughout the rest of the paper.
It also leads naturally to a duality relation between groups and representations.

Quantum mechanics on the U(N) manifold
The correspondence of irreps of SU(N) and N-fermion energy eigenstates can be most naturally established by considering the action of a particle moving freely on the group manifold U(N) [19]. The Lagrangian is the kinetic energy of the particle where overdot signifies time derivative. This is essentially a unitary matrix model.
One approach to find the energy eigenstates of the corresponding Hamiltonian is to note that, upon defining momenta and quantizing, the Hermitian matrix operators Left and right multiplications of U by unitary matrices are generated by L and R, respectively, which conforms with their commutation relations. L and R have equal and nations ∑ a r aa (U) = Tr r (U) = χ r (U), the latter being the character of the irrep. We recover the same spectrum, with states corresponding to irreps r of SU(N), but now with degeneracy one due to the trace.
The condition L + R = 0, or, classically, the vanishing of the matrix commutator [U,U] = 0, implies (using the equation of motion for U) that U(t) and U(0) commute as matrices. Hence, U(t) can be diagonalized by a time-independent unitary conjugation. Implementing this at the level of the Lagrangian by writing and similarly for other sets of N quantities), ∆(z) is given by The prefactor (z 1 · · · z N ) − N−1 2 is a pure phase introduced to make ∆(z) invariant under uniform shifts of the x i 's.
Upon incorporating one factor ∆(z) into the wavefunction and free fermion energy eigenstates. In particular, the fermion energy eigenstate corresponding to an irrep r and energy E = C 2 (r) is with χ r (z) = Tr r(U) the character of irrep r in terms of the eigenvalues z j of U.

Correspondence with Young tableaux
It will be useful to establish the correspondence of fermion states and Young tableaux of SU(N). The single-particle spectrum on the circle consists of discrete momentum eigenstates with eigenvalue k = 0, ±1, ±2, . . . and energy E k = k 2 /2. An N-fermion energy eigenstate corresponds to filling n of the single-particle states with fermions.
The (unnormalized) wavefunction corresponding to this state is given by the Slater The total momentum of the fermions k = k 1 + · · · + k n corresponds to the states picking up a phase e ick upon the shift x j → x j + c, that is, upon U → e ic U. It thus represents the U(1) charge of the state. We may shift all momenta by a constant, changing k and the U(1) charge without affecting the SU(N) part of the states, which can then be labeled by the Alternatively, we can neutralize the U(1) charge by introducing the , similarly to the prefactor introduced in ∆(z). With this additional prefactor, (2.7) maps to (2.4) for the singlet representation for which The final step involves expressing k j − k N in terms of new "bosonic" variables ℓ j as The non-negative, ordered integers ℓ j represent the lengh of rows j = 1, 2, . . . , N − 1 of the Young tableau of the irrep corresponding to the fermionic state. The condition ℓ N = 0 is the statement that Young tableau columns of length N give a singlet, contributing at most to the U(1) charge, and are eliminated. Fig. 1 depicts an example of a fermionic state equivalent to a Young tableau with five rows. As stated, the total momentum and energy of the fermion state map to the U(1) charge and the quadratic Casimir of the representation.
The transition from k j to ℓ j is, in fact, bosonization, in the sense that two or more ℓ i 's can be equal. For reference, we also record the different parametrization The non-negative integers m j are unrestricted and thus provide a "field theoretical" bosonization of the system. They can be thought of as corresponding to excitations in the first N − 1 positive modes j of a second quantized boson field φ(x).

Composition of representations
We will now present the method of composing irreps and derive explicit expressions for the case of composition of fundamentals. We will also work out explicitly some examples of low rank groups.

General method
The composition of representations, and their decomposition into irreps, becomes particularly convenient in the fermion picture.
Consider the direct product r 1 × r 2 of two (possibly reducible) representations r 1 and r 2 . The basic relation implies, through (2.5 and 2.6), that their corresponding fermion states are related as The advantage of the fermionic representation is that from the wavefunction ψ(z) we can simply read off the irreducible components by expanding in monomials z k 1 1 · · · z k N N and mapping them to Young tableaux. This will be demonstrated in detail in the upcoming sections for the composition of a large number of irreps. We record here the characters for the fundamental f , doubly symmetric (a single row with two boxes) s, doubly antisymmetric (two rows with one box each) a, and adjoint ad irreps To obtain the above note that the singlet corresponds to For the fundamental representation we increase k 1 by one, leaving the other k i 's intact.
For the symmetric representation we change k 1 by two and for the antisymmetric representation we change k 1 as well as k 2 by unity. Then, from (2.6), (2.7) and (2.4) we obtain the first three lines above. The conjugate to the fundamental representation χ f is just χ f with the z i replaced by their inverses. By decomposing χ fχ f = χ ad + 1 we obtain χ ad as above. These will constitute basic building blocks when we examine the product of a large number of such irreps.

Composition of fundamentals
The composition of two irreps r ⊗ f , one of which, specifically f , is the fundamental, is particularly intuitive in the fermionic picture. It corresponds to exciting one of the fermions in the state representing the original irrep by one unit of momentum. This gives, in principle, N resulting states, according to which of the N fermions is excited.
However, excitations that would make two fermions occupy the same momentum level are impossible.
It is easy to see that this reproduces the usual rules for Young tableau composition, adding a single box. Adding one unit to k i amounts to increasing in which case this is not possible. Further, exciting the last fermion k N by 1 (when this is possible) amounts to subtracting 1 from all ℓ j , that is, removing a column of length N from the Young tableau.
The fermion wavefunction of the state corresponding to r ⊗ f is, using (3.2) and (3.3), simply This can be used to obtain the fermionic state corresponding to the composition of several (n in number) fundamental irreps f . The original, singlet state is simply ∆(z) and an iteration of the above formula yields The irreducible components of the resulting representation can be "read off" from the terms of ψ N,n (z) expanded as a polnomial in z j : the coefficient of the terms z k 1 1 · · · z k N N such that k 1 > k 2 > · · · > k N , represents the multiplicity of the irrep corresponding to Young tableau lengths ℓ j = k j − k N + j − N, according to (2.8). The remaining terms, making ψ N,n properly fermionic, can be ignored.

Elementary example: Fundamental SU(2) irreps
As a simple example, we can apply the above method to the case of SU(2) and make contact with previous results. In this case (3.5) becomes simply with k 1 > k 2 can be isolated by writing z 1 = zx, z 2 = zx −1 , which gives z k 1 1 z k 2 2 = z k 1 +k 2 x k 1 −k 2 . So powers of z count total U(1) charge, and focusing on terms x k for k > 0 with j being the spin of the irrep. We obtain with the sum running in integer decrements from n/2. The coefficient of the term is the multiplicity of spin j in the decomposition, thus reproducing the result of [5]. 2 A direct interpretation of (3.7) is obtained by focusing on the factor (x + x −1 ) n . Its expansion in x k produces the number of states D n,k at spin S z = k/2. Irreps with total spin j k/2 contribute one state at S z = k/2, so we can isolate the number of irreps at j = k/2 by subtracting the contribution from j k/2 + 1, that is, D n,k+2 . Given that ℓ = 2j, we have which is precisely what the prefactor x − x −1 is achieving, upon identifying k = ℓ + 1.
The large-n, large-j limit is interesting, and it is instructive to compute it in two alternative ways. j e −2j 2 /n , (3.11) which coincides, for j ≫ 1, with the result of [8,5]. Note that the distribution (3.11), for fixed n, has a maximum at j = √ n/2 with a value of order O 2 n /n .
An alternative method is based on random walks, leading to a "diffusion process," which does not rely on the exact formula for d n,j and is generalizable to higher groups.
We focus on the factor (x + x −1 ) n in (3.7) leading to D n,k . Each additional composition with a spin 1 2 shifts the value of k by either +1 or −1 and similarly the corresponding values for ℓ. Explicitly, D n+1,k = D n,k+1 + D n,k−1 . (3.12) In our case D n,k = n n/2 + k/2 which indeed satisfies (3.12).
We will view (3.12) as a discrete evolution process and, for large n and k, we will 2 One may readily check that indeed n/2 ∑ j=j min (2j + 1)d n,j = 2 n , where j min equals 0 or 1/2 depending on whether or not n is even or odd. take the continuum limit. To do this, we set D n,k to a continuous function D(n, k) end extract an explicit scaling factor 2 n to account for the doubling of states per extra added spin Then, a Taylor expansion of F(n, k) in (3.12) to leading orders in n, k implies that this obeys the heat equation Note that the extra factor 2 n in (3.13) ensures that there is a smooth continuum limit of (3.12), and F(n, k) becomes a continuous distribution. The solution of the above equation satisfying the initial condition D 0,k = δ 0,k , that is, F(0, k) = δ(k), is given by for n 0 and k ∈ (−∞, ∞), and thus Then (3.9) becomes in this limit Finally, changing variable to j = ℓ/2 and taking into account that in the continuum limit d n,j becomes a distribution, so that d(n, j)dj = d(n, ℓ)dℓ = 2d(n, ℓ = 2j)dj, giving an extra factor of 2, we reproduce (3.11). This implies the proper normalization where now the integration is over physical (positive) values of j, with the multiplicity 2j + 1 ≃ 2j (valid for large j) included.

Fundamental SU(3) irreps
A more nontrivial case is that of the composition of n fundamental SU (3) irreps. For

SU(3), the fermion state (3.5) becomes
(as in the SU(2) case, we omitted the prefactor in the Vandermonde determinant (2.4).) To isolate terms z k 1 1 z k 2 2 z k 3 3 with k 1 > k 2 > k 3 we adopt the parametrization for which and focus on terms in the wavefunction with positive, ordered powers of x and y.
Clearly z counts the U(1) charge, so we can put z = 1. We also divide the wavefunction by x 2 y to make its terms proportional to x ℓ 1 y ℓ 2 . The reduced form of the wavefunction becomes The generating function for the irrep content of the product of n fundamentals of

SU(3) can then be expressed as
with d n;ℓ 1 ,ℓ 2 the multiplicity of the irrep {ℓ 1 , ℓ 2 } in the product, and | + denoting that only positive, ordered powers x ℓ 1 y ℓ 2 are kept. As a demonstration, we list the nonvanishing d n;ℓ 1 ,ℓ 2 for the first few values of n which reproduce the correct decomposition of (⊗ f ) n . The coefficients d n;ℓ 1 ,ℓ 2 for general n can be calculated explicitly, but we defer their derivation to the general SU(N) case, where we will present a complete treatment.

Composition of many fundamental SU(N) irreps
The fermion method described in the previous sections can be used to derive explicit combinatorial expressions for the components of the decomposition of an arbitrary number n of fundamental U(N) irreps. We will work directly with the fermion momenta k i , leaving the transition to Young tableau row lengths ℓ i in (2.8) for the end.

Combinatorial expressions
The fermion state for the product of n SU(N) irreps is given by (3.5) and gives the generating function of irreps in the momentum parametrization k (recall that we use boldface letters for sets of N quantities such as k 1 , . . . , k N ) The ψ N,n is antisymmetric in the z i , and as a consequence the d n;k appearing above are fully antisymmetric in the k i . When the k i are in decreasing order, d n;k gives the multiplicity of the irrep labeled by k 1 > · · · > k N .
To derive an explicit combinatorial expression for the multiplicity we first focus on the coefficients produced by the term (z 1 + · · · + z N ) n , denoted by D n;k . We have that is, where we used the standard expression for the multinomial coefficient. The Vandermonde factor in (4.1) (again, omitting its phase prefactor) contains N! terms, each augmenting the powers of z i in each term in the expansion. This corresponds to shifting the indices k 1 , . . . , k N in the D n;k 1 ,...,k N summation by the opposite of the corresponding powers; that is, with S i being discrete shift operators acting as Observing that we see that for any polynomial function P(k 1 , . . . , k N ) the Vandermonde of shift operators ∏ j>i (S i − S j ) acting on the last term in (4.6) will produce a polynomial of degree N(N − 1)/2 in the numerator with leading coefficient 1. Since it must also be fully antisymmetric in the k i , the only such polynomial is the Vandermonde of k i . We obtain The above is written in the parametrization with a fixed total momentum. To restore invariance under total momentum shifts k i → k i + c, we re-express the k i as which reduce to k i when the constraint k = n + N(N − 1)/2 is satisfied, but otherwise are invariant under the common shift k i → k i + c, and thus k → k + cN. We obtain The transcription in terms of Young tableau row lengths Formulae (4.10) and (4.11) hold whenever the N-ality condition ℓ = n (mod N), or The above multiplicities reproduce, as they should, the total number of states as where dim(k) is the dimension of irrep {k 1 , . . . , k n }, given by where ∆(k) is the Vandermonde determinant in (2.4) but without the prefactor. To show this, we act on both sides of (4.1) with the Vandermonde of derivative operators . The term ∂ n 1 1 · · · ∂ n N N in ∆(∂), with n 1 , . . . , n N a permutation of 0, 1, . . . , N−1, gives a nonzero result only upon acting on the corresponding term z n 1 1 · · · z n N N in ∆(z), equal to 0! 1! · · · (N − 1)!. Accounting for all N! such terms we have and thus On the right hand side, we can antisymmetrize ∏ i z k i i into the full Slater determinant (2.7), taking advantage of the antisymmetry of d n;k 1 ,...,k N , and restrict the summation to k 1 >· · ·>k N . Acting with ∆(∂) produces a symmetric polynomial in z i with coefficients polynomials in k i , and upon setting z i = 1 it gives an antisymmetric polynomial in k i of degree N(N−1)/2, which is necessarily proportional to the Vandermonde ∆(k).
By acting on the ground state k i = N − i, in which case the Slater state is simply the Vandermonde ∆(z), and using (4.15), the overall coefficient is found to be N!, for a result of N! ∆(k). Equating left and right hand side results we obtain (4.16) which, upon using (4.13), reproduces (4.12).
We note that the above procedure can be interpreted as a random walk on the (N−1)dimensional lattice spanned by x i = k i − k i+1 , i = 1, . . . , N − 1 with starting point Each additional term in (4.2) increasing n to n + 1 can be interpreted as one of N possible steps on the lattice taking is implementing this image method: the term D n;k 1 ,...,k N corresponds to unrestricted walks starting from x 1 = · · · = x N = 0, and the N! terms in the expansion of the product of shift operators generates the images. This is a useful picture, especially in the large-n limit where the random walk essentially becomes a Brownian motion obeying a diffusion equation, but we will not elaborate it further here.
We conclude this section by pointing out that the number of SU(N) irreps d n;ℓ can be interpreted as an N-dimensional Catalan pyramid, just as the corresponding number of SU(2) irreps maps to the standard Catalan triangle numbers. We can define the N-dimensional Catalan pyramid C(ℓ 1 , . . . , ℓ N ) as the number of words with ℓ 1 letters X 1 , ℓ 2 letters X 2 , . . . , ℓ N letters X N such that no initial word has more X i than X i+1 for all i = 1, 2, . . . , N − 1 (N = 2 gives the usual Catalan triangle). This clearly reproduces the rules of composing n = ℓ 1 +ℓ 2 +· · ·+ℓ N SU(N) fundamental irreps, X i representing boxes in row i, and maps to our d n;ℓ . We can also define truncated Catalan pyramids above Catalan structures also clearly admit an interpretation as the number of random walks of the type described in the previous paragraph on the (N − 1)-dimensional integer lattice wedge x 1 x 2 · · · x N−1 0.

Momentum density and a group duality
We conclude by giving a "second quantized" expression for the d n,k 1 ,...,k N that is useful in the large-N, n limit. Thinking of the k i as a distribution of fermions on the positive momentum lattice s = 0, 1, . . . , we define the discrete momentum density of fermions ρ s equal to 1 on points s of the momentum lattice where there is a fermion and zero elsewhere; that is, In the above, M is a cutoff momentum that can be chosen arbitrarily as long as it is bigger than all the k i . Then (4.8) can be written as The integer M could in principle be taken to infinity. However, keeping it finite serves to demonstrate an interesting particle-hole duality of the formulae. Definẽ That is, in the decomposition of the tensor product of n fundamentals of SU(N), the multiplicity of any given irrep is the same as the one for its dual irrep in the product of n fundamentals of SU(M − N + 1). Note that this relation holds for any M such This duality between SU(N) and SU(M − N + 1) can be turned into a self-duality if This will be guaranteed, e.g., in the case n < N. Then (4.22) states that dual irreps in the decomposition come with equal multiplicities.
The above group duality manifests as a particle-hole duality in the fermion density picture, but it can also be understood in terms of Young tableau composition rules: when adding one additional box, the rules are symmetric with respect to rows or columns, the only difference arising when a column of length N is formed, which is then eliminated. Therefore, starting with the singlet and adding single boxes (fundamentals), the obtained set of irreps will be symmetric under row-column exchange as long as n < N. This also shows that a similar duality will emerge in the composition of any number of self-dual irreps, such as, e.g., ℓ 1 = 2, ℓ 2 = 1.

Large-n limit
The limit n ≫ 1 is similar to the large-n limit in the SU(2) case. We can use either a Stirling approximation in the combinatorial formula for d n;ℓ 1 ,...,ℓ N−1 or d n;k 1 ,...,k N , or a random walk approach as in the SU(2) case. Either method leads to the desired result.
Here we detail the derivation using the Stirling approximation n! ≃ √ 2πn n+1/2 e −n , n ≫ 1 , (4.23) in the exact expression (4.10) in which the total momentum is not fixed and shift invariance of the k i 's is manifest. Assuming k i ≫ 1, we obtain the intermediate expres-sion towards a continuous distribution d n;k 1 ,...,k N → d(n; k) (4.24) In the n ≫ 1 limit, a saddle point expansion in the product above reveals that Using the identity , and taking the n ≫ 1 limit we obtain and (4.25) becomes an expression still invariant under uniform shifts of the k i 's. We see that in the large-n limit the distribution of row lengths becomes Gaussian with a polynomial prefactor, and typical row lengths of the obtained representations scale as √ n.
The above result can be expressed in terms of the dimension and quadratic Casimir of each irrep r = (k 1 , . . . , k N ), given by

Coupling to magnetic fields
As a physical application of the results of this paper, we examine the partition function of a system of many fundamental SU(N) "spins" (or better, flavors) coupled to an external non-abelian magnetic field. This will require the evaluation of traces of exponentials of the total flavor operators in arbitrary irreps. If the full Hamiltonian is the one described by (5.1), the flavors are not interacting and the full partition function is the n th power of the single-flavor partition function  The matrix U can be written as x j H j V −1 , (5.8) where θ a are group parameters, V is a diagonalizing matrix, and H j are the Cartan generators, including the U(1) part, in the diagonal basis with matrix elemens (H j ) ik = δ ij δ ik . (5.9) Therefore, the energy eigenstates provide the traces (characters) x j H j . (5.10) On the other hand, as explained in section 2, the energy eigenstates are given as The irrep r with Young tableau lengths ℓ j maps to the fermionic momenta k j as expressed in (2.8).
Note that, (5.12) contains one extra magnetic variable, since there is also a term coupling to the U(1) charge of U (the identity matrix). This can easily be eliminated by shifting all k j by a constant such that k 1 + · · · + k N = 0 which neutralizes the U(1) charge. Alternatively, we can simply choose magnetic parameters B j such that B 1 + · · · + B j = 0, ensuring that the coupling to the U(1) part vanishes. Further, the magnetic variables B j appearing in the trace (5.12) couple to H j in the specific basis (5.9). Moving to any other basis, and in particular to one where the H j are traceless, would simply amount to replacing B j by appropriate linear combinations.
The above traces constitute the basic components needed to determine the statistical mechanics of any collection of flavors with Hamiltonian of the type considered in this section.

Concluding remarks
The SU(N) decomposition problem is solvable in explicit form, at least for the direct product of many "flavors" (fundamental irreps). The results in this paper open the road to various applications and generalizations.
Perhaps the most physically relevant question is the thermodynamics and phase transitions of a large collection of SU(N) flavors in the presence of an external nonabelian magnetic field. This is in principle accessible from the results of the present work. For interactions of ferromagnetic type, the system can be shown to exhibit ferromagneticlike phase transition, but its phase structure and properties turn out to be quite complicated and nontrivial, exhibiting qualitatively different behavior compared to the standard SU(2) ferromagnet. This system is treated in the follow-up paper [20].
Another interesting issue, relevant to various theoretical physics contexts, is the scaling properties of the system as the number of flavors n and the size of the group N both go to infinity. The simplifications and special properties of the limit of large N are known, and have been extensively used in particle physics [15][16][17]. In our case, the existence of the additional parameter n enriches the structure, and it turns out that different limits can be reached according to the relation between n and N as they both go to infinity. Again, in these limits phase transitions appear depending on the relation of the two parameters. This will also be investigated in a future publication.
Finally, the decomposition of the direct product of irreps other than the fundamental can be considered. The relevant formulae quickly grow complicated and unwieldy (see (3.3)), but the thermodynamic properties and the various large-n, N limits could still be analytically accessible, and are topics for future investigation.